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

[Olivier Nuñez]

Bueno, se parece a un problema de separación.

>   dat <- data.frame(x1 = rep(c(0,1,1,0), each = n),
+                     x2= rep(c(0,1), times = n))
> 
>   odds <- exp( -1 - 4 * dat$x1 + 2*dat$x2 - 1 * dat$x1 * dat$x2 )
> 
>   pr <- odds/(1+odds)
>     dat$y <- rbinom(n,1,pr)
> table(dat)
, , y = 0

   x2
x1    0   1
  0 730 278
  1 730 278

, , y = 1

   x2
x1    0   1
  0 270 722
  1 270 722



En otras palabras, conociendo "y" y "x2" es "muy difícil" averiguar si "x1" tomo el valor 0 o 1.
Por lo tanto, resulta difícil evaluar tanto el efecto principal de esta ultima variable como su interacción con x1. 
Te aconsejo optar por un diseño más adecuado y utilizar un regresión logística penalizada (ejemplo: función bayesglm del paquete arm).
Un saludo. Olivier

----- Mensaje original -----
De: "Carlos J. Gil Bellosta" <cgb en datanalytics.com>
Para: "Emilio Torres Manzanera" <torres en uniovi.es>
CC: "r-help-es" <r-help-es en r-project.org>
Enviados: Jueves, 22 de Enero 2015 13:46:44
Asunto: Re: [R-es] Simulación de modelo logit con interacción

Hola, ¿qué tal?

Compara los dos gráficos siguientes:

set.seed(1234)

n <- 100000

logisticsimulation <- function(n){
  dat <- data.frame(x1 = rep(c(0,1), each = n),
                    x2 = rep(c(0,1), times = n))
  odds <- exp(-1 - 4 * dat$x1 + 2*dat$x2 - 1 * dat$x1 * dat$x2 )
  dat$x1 <- factor(dat$x1)
  dat$x2 <- factor(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)

res <- as.data.frame(res)

plot(res)


logisticsimulation <- function(n){
  dat <- data.frame(x1 = rep(c(0,1), each = n),
                    x2 = rep(c(0,1), times = n))
  odds <- exp(-1 - 4 * dat$x1 + 7*dat$x2 - 1 * dat$x1 * dat$x2 )
  dat$x1 <- factor(dat$x1)
  dat$x2 <- factor(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)

res <- as.data.frame(res)

plot(res)


El primero es tu caso. El segundo, uno menos extremo. Verás que tienes
problemas de estabilidad en el coeficiente de x2 porque el coeficiente
es demasiado grande. Tal vez necesites un n mayor para estimar tus
coeficientes bien si tienen ese rango. Un x2 distinto de cero es toda
una rareza.

Salud,

Carlos J. Gil Bellosta
http://www.datanalytics.com

El día 22 de enero de 2015, 12:28, Emilio Torres Manzanera
<torres en uniovi.es> escribió:
> 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

_______________________________________________
R-help-es mailing list
R-help-es en r-project.org
https://stat.ethz.ch/mailman/listinfo/r-help-es