diff --git a/scripts/ce089-06.R b/scripts/ce089-06.R index 69bfd1f039b0005982caa7bacf49e9a79b488ad0..3107a37d2f24a2cc2961f73e362085ae116f75d1 100644 --- a/scripts/ce089-06.R +++ b/scripts/ce089-06.R @@ -1,4 +1,5 @@ #----------------------------------------------------------------------- +# 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. #-----------------------------------------------------------------------