Hola, Tengo una base de datos con estructura jerárquica en la que quiero clasificar observaciones en distintas categorías. En el caso más simple, tengo una variable con dos categorías (variable Y1) y dentro de cada una de ellas hay otras dos categorías (variable Y2). Además tengo una variable explicativa cuantitativa discreta X. El banco de datos sería de este tipo: X Y1 Y2 5 0 1 9 0 0 2 0 1 8 0 0 9 0 0 6 0 0 4 1 0 5 1 1 5 1 1 6 1 1 0 1 1 6 1 0 El enfoque estadístico que quiero plantear es un modelo lineal generalizado mixto con efectos aletorios para tener en cuenta esta estructura jerárquica o multinivel. Para analizar si este enfoque es adecuado he simulado unos datos. En primer lugar, trato de clasificar las observaciones en las categorías de Y1 mediante un glm y obtengo estimadores de los coeficientes que son muy similares a los valores simulados, como era esperable. Sin embargo, cuando trato de utilizar un glmer con un intercepto aleatorio y pendientes fijas para tratar de distinguir entre categorías, obtengo estimadores muy alejados de los valores simulados. He leído mucha de la bibliografía específica (por ejemplo, los capítulos relacionados con los glmer de los libros de Bates o Goldstein) y he buscado por internet pero no logro encontrar el error en mi procedimiento. Adjunto el código de la simulación que llevo a cabo en R por si alguien me pudiera hacer alguna indicación. Muchas gracias de antemano. Un saludo, Guillermo ------------ próxima parte ------------ A non-text attachment was scrubbed... Name: simulation_glmer.R Type: application/octet-stream Size: 2498 bytes Desc: URL: <https://stat.ethz.ch/pipermail/r-help-es/attachments/20150610/4f2e37e6/attachment.obj>
Guillermo,
me temo que en tu simulación, el enfoque multinivel carezca de sentido.
Ten en cuenta que en este tipo de modelo la agrupación de los datos, es decir el
segundo nivel en la jerarquía,
no puede ser en sólo 2 categorías o grupos, sino en un numero considerable de
grupos que justifique el análisis de la variabilidad entre dichos grupos.
Prueba más bien la siguiente simulación:
set.seed(100)
beta_0 <- -0.9
beta_1 <- 0.02
sigma <- 1
q=200; n=50
re <- rnorm(q)*sigma # efectos aleatorios asociado al grupo
u <- rep(re,each=n)
x <- floor(runif(n*q,min=0,max=10))
p <- exp( beta_0 + beta_1*x + u) / (1+ exp( beta_0 + beta_1*x + u))
y <-rbinom(n*q,1,p)
datos=data.frame(y,u,grupo=rep(1:q,each=n))
fit=glmer(y~x+(1|grupo),data=datos,family=binomial)
summary(fit)
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: y ~ x + (1 | grupo)
Data: datos
AIC BIC logLik deviance df.resid
11603.3 11625.0 -5798.7 11597.3 9997
Scaled residuals:
Min 1Q Median 3Q Max
-2.8102 -0.6685 -0.4591 0.8800 4.0097
Random effects:
Groups Name Variance Std.Dev.
grupo (Intercept) 0.9758 0.9878
Number of obs: 10000, groups: grupo, 200
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.911938 0.082754 -11.020 < 2e-16 ***
x 0.027207 0.008101 3.359 0.000783 ***
---
Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
----- Mensaje original -----
De: "Guillermo Vinue" <Guillermo.Vinue en uv.es>
Para: r-help-es en r-project.org
Enviados: Miércoles, 10 de Junio 2015 13:10:26
Asunto: [R-es] Duda glmer
Hola,
Tengo una base de datos con estructura jerárquica en la que quiero
clasificar observaciones en distintas categorías.
En el caso más simple, tengo una variable con dos categorías (variable
Y1) y dentro de cada una de ellas hay otras dos categorías (variable
Y2). Además tengo una variable explicativa cuantitativa discreta X.
El banco de datos sería de este tipo:
X Y1 Y2
5 0 1
9 0 0
2 0 1
8 0 0
9 0 0
6 0 0
4 1 0
5 1 1
5 1 1
6 1 1
0 1 1
6 1 0
El enfoque estadístico que quiero plantear es un modelo lineal
generalizado mixto con efectos aletorios para tener en cuenta esta
estructura jerárquica o multinivel.
Para analizar si este enfoque es adecuado he simulado unos datos. En
primer lugar, trato de clasificar las observaciones en las categorías de
Y1 mediante un glm y obtengo estimadores de los coeficientes que son muy
similares a los valores simulados, como era esperable.
Sin embargo, cuando trato de utilizar un glmer con un intercepto
aleatorio y pendientes fijas para tratar de distinguir entre categorías,
obtengo estimadores muy alejados de los valores simulados.
He leído mucha de la bibliografía específica (por ejemplo, los capítulos
relacionados con los glmer de los libros de Bates o Goldstein) y he
buscado por internet pero no logro encontrar el error en mi procedimiento.
Adjunto el código de la simulación que llevo a cabo en R por si alguien
me pudiera hacer alguna indicación.
Muchas gracias de antemano.
Un saludo,
Guillermo
_______________________________________________
R-help-es mailing list
R-help-es en r-project.org
https://stat.ethz.ch/mailman/listinfo/r-help-es
Muchas gracias por tu respuesta y el código Olivier, me es de mucha ayuda. Un saludo, Guillermo> Guillermo, > > me temo que en tu simulación, el enfoque multinivel carezca de sentido. > Ten en cuenta que en este tipo de modelo la agrupación de los datos,es decir el segundo nivel en la jerarquía,> no puede ser en sólo 2 categorías o grupos, sino en un numeroconsiderable de grupos que justifique el análisis de la variabilidad entre dichos grupos.> Prueba más bien la siguiente simulación: > > > set.seed(100) > beta_0 <- -0.9 > beta_1 <- 0.02 > sigma <- 1 > q=200; n=50 > re <- rnorm(q)*sigma # efectos aleatorios asociado al grupo > u <- rep(re,each=n) > x <- floor(runif(n*q,min=0,max=10)) > p <- exp( beta_0 + beta_1*x + u) / (1+ exp( beta_0 + beta_1*x + u)) > y <-rbinom(n*q,1,p) > datos=data.frame(y,u,grupo=rep(1:q,each=n)) > > fit=glmer(y~x+(1|grupo),data=datos,family=binomial) > summary(fit) > > > Generalized linear mixed model fit by maximum likelihood (Laplace > Approximation) [glmerMod] > Family: binomial ( logit ) > Formula: y ~ x + (1 | grupo) > Data: datos > > AIC BIC logLik deviance df.resid > 11603.3 11625.0 -5798.7 11597.3 9997 > > Scaled residuals: > Min 1Q Median 3Q Max > -2.8102 -0.6685 -0.4591 0.8800 4.0097 > > Random effects: > Groups Name Variance Std.Dev. > grupo (Intercept) 0.9758 0.9878 > Number of obs: 10000, groups: grupo, 200 > > Fixed effects: > Estimate Std. Error z value Pr(>|z|) > (Intercept) -0.911938 0.082754 -11.020 < 2e-16 *** > x 0.027207 0.008101 3.359 0.000783 *** > --- > Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1 > > > > ----- Mensaje original ----- > De: "Guillermo Vinue" <Guillermo.Vinue en uv.es> > Para: r-help-es en r-project.org > Enviados: Miércoles, 10 de Junio 2015 13:10:26 > Asunto: [R-es] Duda glmer > > Hola, > > Tengo una base de datos con estructura jerárquica en la que quiero > clasificar observaciones en distintas categorías. > > En el caso más simple, tengo una variable con dos categorías (variable > Y1) y dentro de cada una de ellas hay otras dos categorías (variable > Y2). Además tengo una variable explicativa cuantitativa discreta X. > > El banco de datos sería de este tipo: > > X Y1 Y2 > 5 0 1 > 9 0 0 > 2 0 1 > 8 0 0 > 9 0 0 > 6 0 0 > 4 1 0 > 5 1 1 > 5 1 1 > 6 1 1 > 0 1 1 > 6 1 0 > > El enfoque estadístico que quiero plantear es un modelo lineal > generalizado mixto con efectos aletorios para tener en cuenta esta > estructura jerárquica o multinivel. > > Para analizar si este enfoque es adecuado he simulado unos datos. En > primer lugar, trato de clasificar las observaciones en las categorías de > Y1 mediante un glm y obtengo estimadores de los coeficientes que son muy > similares a los valores simulados, como era esperable. > > Sin embargo, cuando trato de utilizar un glmer con un intercepto > aleatorio y pendientes fijas para tratar de distinguir entre categorías, > obtengo estimadores muy alejados de los valores simulados. > > He leído mucha de la bibliografía específica (por ejemplo, los capítulos > relacionados con los glmer de los libros de Bates o Goldstein) y he > buscado por internet pero no logro encontrar el error en mi procedimiento . > > Adjunto el código de la simulación que llevo a cabo en R por si alguien > me pudiera hacer alguna indicación. > > Muchas gracias de antemano. > > Un saludo, > > Guillermo > > _______________________________________________ > R-help-es mailing list > R-help-es en r-project.org > https://stat.ethz.ch/mailman/listinfo/r-help-es > >