Efectivamente, la normalidad tarda en llegar (un problema que merece ser
investigado)
En cualquier caso, parece que ambos diseños dan varianzas asintóticas similares
de los estimadores:
> n=1000
>
> dat <- data.frame(x1=sample(0:1, n,replace=TRUE),x2=sample(0:1,
> n,replace=TRUE))
> dat$intercept=1
> odds <- exp(-1 - 4 * dat$x1 + 7*dat$x2 - 1 *dat$x1* dat$x2 )
> pr <- odds/(1+odds)
> D=diag(pr*(1-pr))
> X=as.matrix(dat)
> I=t(X)%*%D%*%X
> solve(I)
x1 x2 intercept
x1 0.4716741 -0.451997095 -0.014952896
x2 -0.4519971 0.470731667 -0.005214237
intercept -0.0149529 -0.005214237 0.020017370
>
> dat$x1= rep(c(0,1), each = n/2)
> dat$x2 = rep(c(0,1), times = n/2)
> odds <- exp(-1 - 4 * dat$x1 + 7*dat$x2 - 1 * dat$x1 * dat$x2 )
> pr <- odds/(1+odds)
>
> D=diag(pr*(1-pr))
> X=as.matrix(dat)
> I=t(X)%*%D%*%X
> solve(I)
x1 x2 intercept
x1 0.45113285 -0.430788206 -0.014755274
x2 -0.43078821 0.451132851 -0.005589371
intercept -0.01475527 -0.005589371 0.020161833
>
----- Mensaje original -----
De: "Carlos J. Gil Bellosta" <c...@datanalytics.com>
Para: "Olivier Nuñez" <onu...@unex.es>
CC: "Emilio Torres Manzanera" <tor...@uniovi.es>, "r-help-es"
<r-help-es@r-project.org>
Enviados: Viernes, 23 de Enero 2015 11:14:44
Asunto: Re: [R-es] Simulación de modelo logit con interacción
Hola, ¿qué tal?
Cierto, cierto, habÃa un error en el código que publiqué. Pero el
diagnóstico es parecido. Cuando los datos se generan con el
coeficiente de x2 igual a 7, los coeficientes estimados tienen una
distribución extraña, bimodal (aparentemente), en lugar de
_normalmente_ distribuida alrededor del 7 como se espera. Supongo que
depende del número de casos en que x2 = 1 e y = 0.
Un saludo,
Carlos J. Gil Bellosta
http://www.datanalytics.com
El dÃa 23 de enero de 2015, 9:54, Olivier Nuñez <onu...@unex.es> escribió:
> Emilio,
>
> espero no haberte generado mucha confusión con mi anterior respuesta.
> El problema no es de separación sino más bien de tamaño muestral.
> Al coger el código de Carlos, obtenÃa que "y" y "x1" eran sistemáticamente
> independiente (la tabla table(dat$y,dat$x1) tiene columnas proporcionales).
> En el diseño de Carlos, el código correcto para simular los datos ha de ser
> dat$y= rbinom(2*n,1,pr).
> Con tu diseño, si elijo un tamaño muestral suficientemente grande, obtengo
> estimaciones razonables de los parámetros de tu modelo:
>
>> res <- logisticsimulation(10000)
>> apply(res,2,median)
> (Intercept) x1 x2 x1:x2
> -0.9955877 -3.9938632 6.9967218 -0.9913839
>
> Un saludo. Olivier
>
>
> ----- Mensaje original -----
> De: "Emilio Torres Manzanera" <tor...@uniovi.es>
> Para: r-help-es@r-project.org
> Enviados: Jueves, 22 de Enero 2015 12:28:14
> Asunto: [R-es] Simulación de modelo logit con interacción
>
> Hola,
> Deseo simular un modelo logit con interacción, estimar sus coeficientes y
> comprobar si son o no parecidos al modelo teórico. Con este ejemplo obtengo
> que los coeficientes estimados no se asemejan mucho a los originales. ¿Se le
> ocurre a alguien cuál es el motivo de esta discrepancia? ¿y cómo solucionarlo?
> Muchas gracias
> Emilio
>
> logisticsimulation <- function(n){
> dat <- data.frame(x1=sample(0:1, n,replace=TRUE),
> x2=sample(0:1, n,replace=TRUE))
> odds <- exp(-1 - 4 * dat$x1 + 7*dat$x2 - 1 *dat$x1* dat$x2 )
> pr <- odds/(1+odds)
> res <- replicate(100, {
> dat$y <- rbinom(n,1,pr)
> coef(glm(y ~ x1*x2, data = dat, family = binomial()))
> })
> t(res)
> }
>
> res <- logisticsimulation(100)
> apply(res,2,median)
> ## (Intercept) x1 x2 x1:x2
> ## -1.0986123 -18.4674562 20.4823593 -0.0512933
>
> DeberÃan salir -1, -4, 7, 1
>
> _______________________________________________
> R-help-es mailing list
> R-help-es@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es
>
> _______________________________________________
> R-help-es mailing list
> R-help-es@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es
_______________________________________________
R-help-es mailing list
R-help-es@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-help-es
