20 Paquete MASS
El paquete MASS de Ripley (2023) se utiliza para estimar modelos glmm por medio de Penalized Quasi-Likelihood (PQL). Al visitar este enlace se encontrará la página de apoyo del paquete, allí se puede consultar el manual de referencia y las viñetas.
20.1 Función glmmPQL
La función glmmPQL es la principal función del paquete glmmMAAS. Esta función sirve para ajustar un glmm y su estructura es la siguiente:
glmmPQL(fixed, random, family, data, correlation, weights,
control, niter = 10, verbose = TRUE, ...)Los principales argumentos de la función son:
formula: es una fórmula similar a la usada en el modelo lineal clásico. Un ejemplo de fórmula seríay ~ 1 + x1 + x2 + (1 + x2 | grupo)con la cual se indican los efectos fijos y los efectos aleatorios del modelo. Más abajo hay una tabla con más detalles sobre la fórmula.random: es una fórmula solo con lado derecho. Si queremos indicar intercepto aleatorio se escribe~ 1 | groupy si se desea intercepto y pendiente aleatoria se escribe~ 1 + x2 | group.data: marco de datos donde están las variables.family: argumento para indicar la distribución de la variable respuesta.
Ejemplo: modelo Bernoulli (binomial) con intercepto aleatorio
En este ejemplo vamos a simular observaciones \(n_i=50\) observaciones para \(G=10\) grupos (en total 500 obs) que tengan la estructura mostrada abajo. El objetivo del ejemplo es ilustrar el uso de la función glmmPQL para estimar los parámetros del modelo.
El vector de parámetros de este modelo es \(\boldsymbol{\Theta}=(\beta_0=4, \beta_1=-8, \sigma_{b0}=2)^\top\).
El código para simular las 500 observaciones se muestra a continuación. Observe que se fijó la semilla en dos ocasiones para que el lector pueda replicar el ejemplo y obtener los mismos resultados.
logit_inv <- function(x) exp(x) / (exp(x) + 1) # Función enlace inversa
ni <- 50
G <- 10
nobs <- ni * G
grupo <- factor(rep(x=1:G, each=ni))
obs <- rep(x=1:ni, times=G)
set.seed(12345)
x <- runif(n=nobs, min=0, max=1)
set.seed(12345)
b0 <- rnorm(n=G, mean=0, sd=sqrt(4)) # Intercepto aleatorio
b0 <- rep(x=b0, each=ni) # El mismo intercepto aleatorio pero repetido
prob <- logit_inv(4 - 8 * x + b0)
set.seed(12345)
y <- rbinom(n=nobs, size=1, prob=prob)
datos <- data.frame(obs, grupo, b0, x, prob, y)Para estimar los parámetros del modelo se usa la función glmmPQL de la siguiente forma.
library(nlme)
library(MASS)
fit1 <- glmmPQL(y ~ x, random = ~ 1 | grupo, niter=20,
family = binomial(link = "logit"), data = datos)La función summary se puede usar sobre el objeto fit1 para obtener una tabla de resumen, a continuación se ilustra el uso y la salida de summary.
## Linear mixed-effects model fit by maximum likelihood
## Data: datos
## AIC BIC logLik
## NA NA NA
##
## Random effects:
## Formula: ~1 | grupo
## (Intercept) Residual
## StdDev: 3.089343 1.195473
##
## Variance function:
## Structure: fixed weights
## Formula: ~invwt
## Fixed effects: y ~ x
## Value Std.Error DF t-value p-value
## (Intercept) 4.864553 1.1476975 489 4.238533 0
## x -7.811406 0.9982813 489 -7.824854 0
## Correlation:
## (Intr)
## x -0.493
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.88958933 -0.20167529 0.02509534 0.24445852 8.17974927
##
## Number of Observations: 500
## Number of Groups: 10Según el resultado anterior \(\hat{\boldsymbol{\Theta}}=(\hat{\beta}_0=4.86, \hat{\beta}_1=-7.81, \hat{\sigma}_{b0}=3.09)^\top\) mientras que el vector real de parámetros es \(\boldsymbol{\Theta}=(\beta_0=4, \beta_1=-8, \sigma_{b0}=2)^\top\).
20.2 What it means for AIC NA in glmmPQL (MASS) summary output?
Esta pregunta fue hecha en StackOverFlow y puede ser consultada en este enlace.
20.3 Two questions on the results from glmmPQL(MASS)
Esta pregunta fue hecha en la comunidad de usuarios de R y puede ser consultada en este enlace.
20.4 How do I interpret the variance of random effect in a generalized linear mixed model?
Esta pregunta fue hecha en CrossValidated y puede ser consultada en este enlace.