Complementa código para taxa de cobertura de IC bootstrap.

scripts/ce089-06.R

 #-----------------------------------------------------------------------
+# Aplicação de bootstrap. Estudo de caso com modelos não lineares.
 
 library(alr3)
 str(turk0)
@@ -28,6 +29,7 @@ n1 <- nls(Gain ~ Int + Ass * (1 - 2^(-A/Mei)),
 cbind(coef(n0), coef(n1))
 
 #-----------------------------------------------------------------------
+# Criando funções que fazem reamostragem dos dados e ajustam o modelo.
 
 n <- nrow(turk0)
 i <- 1:n
@@ -47,6 +49,7 @@ myboot1<- function() {
 }
 
 #-----------------------------------------------------------------------
+# Replicando.
 
 b0 <- replicate(999, myboot0())
 b1 <- replicate(999, myboot1())
@@ -157,7 +160,10 @@ library(latticeExtra)
 
 help(boot, help_type = "html")
 
-str(turk0)
+# dput(as.list(coef(n1)))
+start <- list(Int = 622.958054146589,
+              Ass = 178.251908347242,
+              Mei = 0.097321927254495)
 
 # Criar uma função com dois argumentos: o data.frame original e um vetor
 # que vai representar o índice das linhas usado para reamostrar dos
@@ -165,7 +171,7 @@ str(turk0)
 fitmodel <- function(dataset, index) {
     n0 <- nls(Gain ~ Int + Ass * (1 - 2^(-A/Mei)),
               data = dataset[index, ],
-              start = list(Int = 625, Ass = 180, Mei = 0.1))
+              start = start)
     c(coef(n0), s2 = deviance(n0)/df.residual(n0))
 }
 
@@ -234,7 +240,8 @@ bci
 
 str(bci)
 
-help(boot.ci, help_type = "html")
+#-----------------------------------------------------------------------
+# Abrir as funções para estudar o código fonte.
 
 # Abrir o código fonte para entender.
 boot.ci
@@ -246,6 +253,9 @@ boot:::stud.ci
 boot:::perc.ci
 boot:::bca.ci
 
+#-----------------------------------------------------------------------
+# Fazendo o debugging das funções.
+
 debug(norm.ci)
 norm.ci(b0)
 undebug(norm.ci)
@@ -254,16 +264,79 @@ debug(confint)
 confint(b0, type = "basic")
 undebug(confint)
 
+#-----------------------------------------------------------------------
 #-----------------------------------------------------------------------
 # Estudo de simulação sobre o desempenho dos tipos de intervalo.
 
-# Gerar observações assumindo um modelo linear com resposta Poisson.
+#--------------------------------------------
+# 1. Gerar observações assumindo um modelo linear com resposta Poisson.
+
+simul_model <- function(n = 50,
+                        x = rep(seq(0, 3, length.out = 10), times = 5),
+                        beta = c(-0.5, 1)) {
+    eta <- beta[1] + beta[2] * x
+    mu <- exp(eta)
+    y <- rpois(n = n, lambda = mu)
+    return(data.frame(x = x, y = y))
+}
+
+# Vendo o resultado do modelo.
+with(simul_model(),
+     expr = {
+         plot(y ~ x)
+     })
+
+#--------------------------------------------
+# 2. Para cada novo conjunto de dados, obter os IC com cada método
+# usando um R = 999.
+
+# Função que ajusta e retorna as estatísticas de interesse.
+fitmodel <- function(dataset, index) {
+    g0 <- glm(y ~ x,
+              data = dataset[index, ],
+              family = poisson(link = log))
+    return(c(coef(g0), deviance(g0)))
+}
+
+# Valores usados para os parâmetros.
+beta <- c(-0.5, 1)
 
-# Para cada novo conjunto de dados, obter os IC com cada método usando
-# um R = 999.
+# Corre as simulações bootstrap (reamostragens).
+btst <- boot(data = simul_model(beta = beta),
+             statistic = fitmodel,
+             R = 999)
 
-# Repetir isso 100 vezes.
+# Obtém os intervalos de confiança para uma das estatísticas.
+bci <- boot.ci(btst,
+               type = "all",
+               index = c(2, 2))
+str(bci)
 
-# Verificar a cobertura dos intervalos de confiança.
+# Extração dos limites do IC.
+ics <- sapply(bci[-(1:3)],
+              FUN = function(x) {
+                  # Pega as duas últimas colunas.
+                  x[, ncol(x) - 1:0]
+              })
+ics
+
+# Verificar se contém o não o valor do parâmetro.
+inside <- apply(ics,
+                MARGIN = 2,
+                FUN = function(lim) {
+                    prod(lim - beta[2]) < 0
+                })
+inside
+
+# Comprimento de cada IC.
+apply(ics,
+      MARGIN = 2,
+      FUN = diff)
+
+#--------------------------------------------
+# 3. Repetir isso 100 vezes.
+
+#--------------------------------------------
+# 4. Verificar a taxa de cobertura dos intervalos de confiança.
 
 #-----------------------------------------------------------------------