¿Hay una manera de obtener "efectos marginales" de un objeto ''glmer''
lme4 marginal-effects (3)
Estoy estimando el modelo logit de efectos aleatorios utilizando glmer y me gustaría informar Efectos marginales para las variables independientes. Para los modelos glm , el paquete mfx ayuda a calcular los efectos marginales. ¿Hay algún paquete o función para objetos glmer ?
Gracias por tu ayuda.
Un ejemplo reproducible se da a continuación.
mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")
mydata$rank <- factor(mydata$rank) #creating ranks
id <- rep(1:ceiling(nrow(mydata)/2), times=c(2)) #creating ID variable
mydata <- cbind(mydata,data.frame(id,stringsAsFactors=FALSE))
set.seed(12345)
mydata$ran <- runif(nrow(mydata),0,1) #creating a random variable
library(lme4)
cfelr <- glmer(admit ~ (1 | id) + rank + gpa + ran + gre, data=mydata ,family = binomial)
summary(cfelr)
Esta es una respuesta mucho menos técnica, pero quizás proporciona un recurso útil. Soy un fan del paquete sjPlot que proporciona gráficos de efectos marginales de objetos glmer, como por ejemplo:
library(sjPlot)
sjp.glmer(cfelr, type = "eff")
El paquete ofrece muchas opciones para explorar también los efectos fijos y aleatorios de un modelo glmer. https://github.com/strengejacke/sjPlot
Saludos ben
Mi solución no responde a la pregunta,
"¿Hay una manera de obtener" efectos marginales "de un objeto glmer ",
sino mas bien
"¿Hay alguna forma de obtener coeficientes de regresión logística marginal de una regresión logística condicional con una intersección aleatoria?"
Solo ofrezco este informe porque el ejemplo reproducible proporcionado fue una regresión logística condicional con una intercepción aleatoria y tengo la intención de ser útil. Por favor, no votes abajo; Voy a anotar si esta respuesta se considera demasiado fuera de tema.
El código R se basa en el trabajo de Patrick Heagerty (haga clic en "Ver sin procesar" para ver el pdf) , e incluyo un ejemplo reproducible a continuación de mi versión github de su paquete lnMLE (disculpe las advertencias en la instalación. El paquete de Patrick no CRAN). Estoy omitiendo la salida para todos excepto la última línea, compare , que muestra los coeficientes de efectos fijos de lado a lado.
library(devtools)
install_github("lnMLE_1.0-2", "swihart")
library(lnMLE)
## run the example from the logit.normal.mle help page
## see also the accompanying document (click ''View Raw'' on page below:)
## https://github.com/swihart/lnMLE_1.0-2/blob/master/inst/doc/lnMLEhelp.pdf
data(eye_race)
attach(eye_race)
marg_model <- logit.normal.mle(meanmodel = value ~ black,
logSigma= ~1,
id=eye_race$id,
model="marginal",
data=eye_race,
tol=1e-5,
maxits=100,
r=50)
marg_model
cond_model <- logit.normal.mle(meanmodel = value ~ black,
logSigma= ~1,
id=eye_race$id,
model="conditional",
data=eye_race,
tol=1e-5,
maxits=100,
r=50)
cond_model
compare<-round(cbind(marg_model$beta, cond_model$beta),2)
colnames(compare)<-c("Marginal", "Conditional")
compare
La salida de la última línea:
comparar
Marginal Conditional
(Intercept) -2.43 -4.94
black 0.08 0.15
Intenté el ejemplo reproducible dado, pero tuve problemas con las implementaciones glmer y lnMLE; de nuevo, solo glmer() la salida correspondiente a los resultados de comparación y las advertencias de la llamada glmer() :
##original question / answer... glmer() function gave a warning and the lnMLE did not fit well...
mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")
mydata$rank <- factor(mydata$rank) #creating ranks
id <- rep(1:ceiling(nrow(mydata)/2), times=c(2)) #creating ID variable
mydata <- cbind(mydata,data.frame(id,stringsAsFactors=FALSE))
set.seed(12345)
mydata$ran <- runif(nrow(mydata),0,1) #creating a random variable
library(lme4)
cfelr <- glmer(admit ~ (1 | id) + rank + gpa + ran + gre,
data=mydata,
family = binomial)
El cual dio
Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.00161047 (tol = 0.001, component 2)
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
- Rescale variables?
pero tontamente seguí sin logit.normal.mle a logit.normal.mle , tratando de aplicar el logit.normal.mle al ejemplo dado. Sin embargo, el modelo condicional no converge ni produce estimaciones de error estándar,
summary(cfelr)
library(devtools)
install_github("lnMLE_1.0-2", "swihart")
library(lnMLE)
mydata$rank2 = mydata$rank==2
mydata$rank3 = mydata$rank==3
mydata$rank4 = mydata$rank==4
cfelr_cond = logit.normal.mle(meanmodel = admit ~ rank2+rank3+rank4+gpa+ran+gre,
logSigma = ~1 ,
id=id,
model="conditional",
data=mydata,
r=50,
tol=1e-6,
maxits=500)
cfelr_cond
cfelr_marg = logit.normal.mle(meanmodel = admit ~ rank2+rank3+rank4+gpa+ran+gre,
logSigma = ~1 ,
id=id,
model="marginal",
data=mydata,
r=50,
tol=1e-6,
maxits=500)
cfelr_marg
compare_glmer<-round(cbind(cfelr_marg$beta, cfelr_cond$beta,summary(cfelr)$coeff[,"Estimate"]),3)
colnames(compare_glmer)<-c("Marginal", "Conditional","glmer() Conditional")
compare_glmer
La última línea de la cual revela que el modelo condicional de cfelr_cond no evaluó un modelo condicional, solo devolvió los coeficientes marginales y no los errores estándar.
> compare_glmer
Marginal Conditional glmer() Conditional
(Intercept) -4.407 -4.407 -4.425
rank2 -0.667 -0.667 -0.680
rank3 -1.832 -1.833 -1.418
rank4 -1.930 -1.930 -1.585
gpa 0.547 0.548 0.869
ran 0.860 0.860 0.413
gre 0.004 0.004 0.002
Espero poder solucionar estos problemas. Cualquier ayuda / comentario apreciado. Daré actualizaciones de estado cuando pueda.
Podría usar el ggeffects-package (ejemplos en las package-vignettes ). Entonces, para tu código, esto podría verse así:
library(ggeffects)
# dat is a data frame with marginal effects
dat <- ggpredict(cfelr, term = "rank")
plot(dat)
o usa, como describió Benjamin, el Podía usar sjPlot-package , usando la función plot_model() con plot-type "pred" (esto simplemente envuelve el paquete ggeffects para diagramas de efectos marginales):
library(sjPlot)
plot_model(cfelr, type = "pred", term = "rank")