Emilio Torres Manzanera
2015-Jan-22 11:28 UTC
[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
Hazlo con n = 10000 Estás modelizando una variable binaria con otras 2 binarias: "muy chungo" Un saludo. Isidro Hidalgo Arellano Observatorio Regional de Empleo Consejería de Empleo y Economía http://www.jccm.es> -----Mensaje original----- > De: R-help-es [mailto:r-help-es-bounces en r-project.org] En nombre de > Emilio Torres Manzanera > Enviado el: jueves, 22 de enero de 2015 12:28 > Para: r-help-es en r-project.org > 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
Carlos J. Gil Bellosta
2015-Jan-22 12:46 UTC
[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
<torres en 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 en r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es
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" <cgb en datanalytics.com> Para: "Emilio Torres Manzanera" <torres en uniovi.es> CC: "r-help-es" <r-help-es en 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 <torres en 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 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, 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
Carlos J. Gil Bellosta
2015-Jan-23 10:14 UTC
[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