[R-es] Bootstrap
Jose Luis Cañadas
canadasreche en gmail.com
Jue Dic 27 16:55:59 CET 2012
Hola.
Por lo que entiendo lo que quieres hacer es calcular los coeficientes
mediante bootstrap de un glm. Supongo que para calcular Intervalos de
confianza.
No conozco la librería bootstrap , pero he utilizado la librería boot y
la librería car que te pueden ayudar.
Con la librería car viene la función bootCase que es bastante sencilla
de utilizar.
library(car)
tu.modelo <- glm(y~x1+x2+x3, data=datos, family=binomial)
betahat <- coef(tu.modelo)
betahat.boot <- bootCase(tu.modelo)
summary(betahat.boot)
Y para calcular los I.C
apply(betahat.boot, 2, function(x) quantile(x, c(0.25,0.975) ))
Con la librería boot sería
library(boot)
coeficientes <- function(formula,data,indices){
d <- data[indices,]
fit <- glm(formula,data=d, family=binomial)
return(coef(fit))
}
res.boot <- boot(data=xdata ,statistic=coeficientes, R=2000, formula =
y~x1+x2+x3)
boot.ci(res.boot,index=1, type="basic")#IC para el primer parámetro
boot.ci(res.boot,index=2, type="basic")#IC para el segundo parámetro
Espero que te sirva.
Saludos.
El 27/12/12 16:18, dolors giralt casellas escribió:
> Hola, buenas tardes
>
> estoy intentando hacer un bootstrap de un modelo, pero me da el siguiente
> error:
>
> "Error in FUN(newX[, i], ...) :
> unused argument(s) (list(age = c(33, 47, 49, 56, 60, 64, 64, 66, 68, 69,
> 71, 71, 72, 73, 74, 75, 75, 76, 78, 81, 83, 83, 36, 43, 46, 47, 49, 49, 51,
> 51, 52, 52, 53, 54, 54, 54, 55, 56, 56, 57, 57, 58, 58, 58, 58, 59, 59, 60,
> 61, 62, 63, 64, 65, 65, 66, 66, 67, 68, 69, 69, 71, 71, 71, 71, 71, 72, 72,
> 72, 72, 72, 73, 73, 73, 74, 74, 74, 74, 74, 75, 75, 76, 78, 79, 81, 81, 84,
> 84, 85, 88, 91,"
>
> tengo un data frame con 4 variables y uso las 4 para generar el modelo,
> pero no entiendo porque me da este error (la variable en question es la
> xdata[,1])
>
> muchas gracias!
>
>
> library (bootstrap)
>
> n <- 1794
>
> theta <- function(xdata)
> {
> coef(glm(xdata[,2]~xdata[,1]+xdata[,3]+xdata[,4],xdata ,
> family=binomial))[4]
> }
>
> results <- bootstrap(1:n,1000,theta,xdata)
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-help-es mailing list
> R-help-es en r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es
Más información sobre la lista de distribución R-help-es