Bueno, se parece a un problema de separación. > dat <- data.frame(x1 = rep(c(0,1,1,0), each = n), + x2= rep(c(0,1), times = n)) > > odds <- exp( -1 - 4 * dat$x1 + 2*dat$x2 - 1 * dat$x1 * dat$x2 ) > > pr <- odds/(1+odds) > dat$y <- rbinom(n,1,pr) > table(dat) , , y = 0
x2 x1 0 1 0 730 278 1 730 278 , , y = 1 x2 x1 0 1 0 270 722 1 270 722 En otras palabras, conociendo "y" y "x2" es "muy difícil" averiguar si "x1" tomo el valor 0 o 1. Por lo tanto, resulta difícil evaluar tanto el efecto principal de esta ultima variable como su interacción con x1. Te aconsejo optar por un diseño más adecuado y utilizar un regresión logística penalizada (ejemplo: función bayesglm del paquete arm). Un saludo. Olivier ----- Mensaje original ----- De: "Carlos J. Gil Bellosta" <c...@datanalytics.com> Para: "Emilio Torres Manzanera" <tor...@uniovi.es> CC: "r-help-es" <r-help-es@r-project.org> Enviados: Jueves, 22 de Enero 2015 13:46:44 Asunto: Re: [R-es] Simulación de modelo logit con interacción 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 <tor...@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@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 _______________________________________________ R-help-es mailing list R-help-es@r-project.org https://stat.ethz.ch/mailman/listinfo/r-help-es