[R-es] Simulación de modelo logit con interacción

Olivier Nuñez onunez en unex.es
Vie Ene 23 09:54:15 CET 2015


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" <torres en uniovi.es>
Para: r-help-es en 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 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