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

[Carlos J. Gil Bellosta]

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