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

[eric]

Pero que gran comunidad ... si el mundo funcionara con esta generosidad, el 
estado de cosas sería diferente .... gracias por existir.

Slds, eric.




On 23/01/15 16:24, Emilio Torres Manzanera wrote:
> ¡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 <onu...@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" <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
>
>

-- 
Forest Engineer
Master in Environmental and Natural Resource Economics
Ph.D. student in Sciences of Natural Resources at La Frontera University
Member in AguaDeTemu2030, citizen movement for Temuco with green city standards 
for living

Nota: Las tildes se han omitido para asegurar compatibilidad con algunos 
lectores de correo.

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