[R-es] Simulación de modelo logit con interacción
Jorge I Velez
jorgeivanvelez en gmail.com
Sab Ene 24 01:38:38 CET 2015
Gracias Emilio por el codigo en R.
Hace poco estuve revisando una situacion similar a la que ilustras en tu
mensaje (i.e., seleccion de modelos e inferencia). Dale una mirada a estos
dos articulos:
http://www-stat.wharton.upenn.edu/~berkr/PoSI-submit.pdf
http://statweb.stanford.edu/~ckirby/brad/papers/2013ModelSelection.pdf
Saludos cordiales,
Jorge.-
2015-01-24 6:24 GMT+11:00 Emilio Torres Manzanera <torres en uniovi.es>:
> ¡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 mailing list
> R-help-es en r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es
>
[[alternative HTML version deleted]]
Más información sobre la lista de distribución R-help-es