thr3ads.net - R help es - [R-es] Simulación de modelo logit con interacción [Jan 2015]

If this information is useful, please help other people find it:
Share via:

Carlos J. Gil Bellosta

2015-Jan-23 10:14 UTC

[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 <onunez en 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" <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
>
> _______________________________________________
> R-help-es mailing list
> R-help-es en r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es

Olivier Nuñez

2015-Jan-23 11:23 UTC

head link

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

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" <cgb en datanalytics.com>
Para: "Olivier Nuñez" <onunez en unex.es>
CC: "Emilio Torres Manzanera" <torres en uniovi.es>,
"r-help-es" <r-help-es en 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 <onunez en 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" <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
>
> _______________________________________________
> R-help-es mailing list
> R-help-es en r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es

Emilio Torres Manzanera

2015-Jan-23 19:24 UTC

head link

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

¡Mil gracias!
Parece que el tamaño ¡sí que importa! aunque solo sean 4 coeficientes. Para
estimar los 4 coeficientes con cierta precisión, se necesitan unos 2000 datos.
La próxima vez que vea un estudio sobre si trabaja (sí/no) en función del sexo
(hombre/mujer) y estado civil (casado/soltero), con un tamaño de muestra de
n=400, no sé qué le diré al autor...

Respecto a los coeficientes, la verdad es que eran feos (eso me pasa por copiar
y pegar de otros ejemplos). Pongo abajo el código definitivo que he utilizado
para que mi pobre ordenador pudiera correr la simulación en un tiempo razonable.

¡Gracias de nuevo!
Emilio

## ===========================================================library(dplyr)
## funcion auxiliar
hacertablafrecuencias <- function(x, vars=NULL, freq=NULL, sort=FALSE){
  evaldp <- function (.data, .fun.name, ..., .envir = parent.frame())
    {
      args <- list(...)
      args <- partial_eval(args, env = .envir)
      args <- unlist(args)
      if (is.function(.fun.name))
        .fun.name <- as.character(substitute(.fun.name))
      code <- paste0(.fun.name, "(.data,", paste0(args, collapse =
","),
                     ")")
      eval(parse(text = code, srcfile = NULL))
    }
  if(is.null(vars)) vars <- names(x)
  nms <- c(vars, "freq")
  nmsfreq <- make.names(nms, unique = TRUE)[length(nms)]
  if (is.null(freq)) {
    action <- paste( nmsfreq,"=n()" ) ##quote(n())
  } else {
    if(freq=="freq") nmsfreq <- "freq"
    vars <- vars[ vars!= "freq"]
    action <- paste(nmsfreq, "=sum(",freq,")")
  }
  if(nmsfreq!="freq") warning(paste("freq column:",nmsfreq))
  out <- group_by_(x, .dots = vars) %>%
    evaldp(summarise, action) %>%
      ungroup()
  if (!sort) {
    out
  } else {
    arrange(out, desc(freq))
  }
}

 ## modelo : -1 - 4 * dat$x1 + 7 * dat$x2 - 1 *dat$x1* dat$x2 

logisticsimulation <- function(n){
  dat <- data.frame(x1=rep(c(0:1), each=n),
                    x2=rep(c(0:1), time=n))
  odds <- exp(-1 - 4 * dat$x1 + 7*dat$x2 - 1 *dat$x1* dat$x2 )
  pr <- odds/(1+odds)
  res <- replicate(100, {
    dat$y <- rbinom(2*n,1,pr)
    ttdat <- hacertablafrecuencias(dat)
    coef(glm(y ~ x1*x2, data = ttdat, family = binomial(),weights=ttdat$freq))
  })
  t(res)
}

res <- logisticsimulation(400) # son en realidad 800 datos. Ajusta de pena
apply(res,2,median)
## (Intercept)          x1          x2       x1:x2 
##  -0.9946226  -4.2473363  19.8453136  -0.5429929

res <- logisticsimulation(1000) # son en realidad 2000 datos. Buen ajuste
apply(res,2,median)
## (Intercept)          x1          x2       x1:x2 
##   -1.004794   -4.120370    7.243097   -1.267561






On Friday, January 23 2015, 12:23:25, Olivier Nuñez <onunez en unex.es>
wrote:
> 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" <cgb en datanalytics.com>
> Para: "Olivier Nuñez" <onunez en unex.es>
> CC: "Emilio Torres Manzanera" <torres en uniovi.es>,
"r-help-es" <r-help-es en 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 <onunez en 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" <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
>>
>> _______________________________________________
>> R-help-es mailing list
>> R-help-es en r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-help-es

-- 
================================================Emilio Torres Manzanera
Fac. de Comercio - Universidad de Oviedo
c/ Luis Moya 261, E-33203 Gijón (Spain)
Tel. 985 182 197 email: torres en uniovi.es

R help es - Jan 2015 - Simulación de modelo logit con interacción

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

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

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