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