Probabilidad y estimaciones de efectos mixtos Regresión logística

Primero simulemos algunos datos para una regresión logística con partes fijas y aleatorias:

set.seed(1)
n    <- 100
x    <- runif(n)
z    <- sample(c(0,1), n, replace=TRUE)
b    <- rnorm(2)
beta <- c(0.4, 0.8)
X    <- model.matrix(~x)
Z    <- cbind(z, 1-z)
eta  <- X%*%beta + Z%*%b
pr   <- 1/(1+exp(-eta))
y    <- rbinom(n, 1, pr)

Si solo quisiéramos ajustar una regresión logística sin partes aleatorias, podríamos usar la glmfunción:

glm(y~x, family="binomial")

glm(y~x, family="binomial")$coefficients
# (Intercept)           x 
#  -0.2992785   2.1429825 

O construir nuestra propia función del log-verosimilitud

donde $\Lambda(\eta_i)=\frac{1}{1+\exp(-\eta_i)}$ y $\eta_i=\boldsymbol{X_i' \beta}$ y use optim()para estimar los parámetros que lo maximizan, como en el siguiente código de ejemplo:

ll.no.random <- function(theta,X,y){
  beta <- theta[1:ncol(X)]
  eta  <- X%*%beta
  p    <- 1/(1+exp(-eta))
  ll   <- sum( y*log(p) + (1-y)*log(1-p) )
  -ll
}

optim(c(0,1), ll.no.random, X=X, y=y)

optim(c(0,1), ll.no.random, X=X, y=y)$par
# -0.2992456  2.1427484

que, por supuesto, proporciona las mismas estimaciones y maximiza la probabilidad logarítmica para el mismo valor. Para efectos mixtos, nos gustaría algo como

library(lme4)
glmer(y~x + (1|z), family="binomial")

Pero, ¿cómo podemos hacer lo mismo con nuestra propia función? Dado que la probabilidad es

y la integral no tiene una expresión de forma cerrada, necesitamos usar la integración numérica como la Cuadratura Gaussiana. Podemos usar el paquete statmodpara obtener algunas cuadraturas, digamos 10

library(statmod)    
gq <- gauss.quad(10)
w  <- gq$weights
g  <- gq$nodes

ACTUALIZACIÓN: Al usar estas ubicaciones en cuadratura y los pesos para ( aquí), podemos aproximar la integral sobre mediante una suma de los términos con sustituido por y todo término multiplicado por los pesos respectivos . Por lo tanto, nuestra función de probabilidad debería ser ahora$g_r$$w_r$$r=1,...,R$$R=10$$b_j$$R$$g_r$$b_j$$w_r$

Además, debemos tener en cuenta la varianza de la parte aleatoria, leí que esto se puede lograr reemplazando en nuestra función con where , por lo que en la función de probabilidad anterior reemplazamos 's con ' s y no 's.$b_j \sim N(0,\sigma_b^2)$$\eta$$\sigma_j \theta_j$$\theta_j \sim N(0,1)$$\theta$$g$$\beta$

Un problema computacional que no entiendo es cómo sustituir los términos, ya que los vectores no tendrán la misma longitud. Pero probablemente no entiendo eso, porque me falta algo crucial aquí o no entendí mal cómo funciona este método.

No vi cómo "los vectores no serán de la misma longitud", aclare su pregunta.

En primer lugar, para la integral con una dimensión inferior a 4, los métodos numéricos directos como la cuadratura son más eficientes que MCMC. Estudié estas preguntas por un tiempo y me complacería discutir este problema con usted.

Para la regresión logística de efectos mixtos, el único Rcódigo explícito que he encontrado es del libro del Prof. Demidenko, Modelos mixtos: teoría y aplicaciones , puede descargar el código a través de la columna "SOFTWARE Y DATOS" en la página web. El logMLEgh()se puede encontrar en \mixed_models_data.zip\MixedModels\Chapter07. No usó el statmodpaquete para obtener las cuadraturas, pero escribió su propia función gauher(). Hay algunos errores menores en el código y los he discutido con el autor, pero todavía es muy útil comenzar desde su código y libro. Puedo proporcionar la versión corregida si es necesario.

Otro problema es que, si desea obtener estimaciones precisas, optim()no es suficiente, es posible que deba usar métodos como la puntuación de Fisher, como en glm().