[R-es] Simulación de modelo logit con interacción
Emilio Torres Manzanera
torres en uniovi.es
Vie Ene 23 20:24:09 CET 2015
¡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
Más información sobre la lista de distribución R-help-es