Probabilidade e estimativas para efeitos mistos Regressão logística

Primeiro vamos simular alguns dados para uma regressão logística com partes fixas e aleatórias:

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)

Se apenas desejássemos ajustar uma regressão logística sem partes aleatórias, poderíamos usar a glmfunção:

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

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

Ou construindo nossa própria função da probabilidade logarítmica

onde $\Lambda(\eta_i)=\frac{1}{1+\exp(-\eta_i)}$ e $\eta_i=\boldsymbol{X_i' \beta}$ e use optim()para estimar os parâmetros que o maximizam, como a seguir código de exemplo:

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, é claro, fornece as mesmas estimativas e maximiza a probabilidade de log para o mesmo valor. Para efeitos mistos, gostaríamos de algo como

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

Mas como podemos fazer o mesmo com nossa própria função? Como a probabilidade é

e a integral não tem expressão de forma fechada, precisamos usar a integração numérica como a quadratura gaussiana. Podemos usar o pacote statmodpara obter algumas quadraturas, digamos 10

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

ATUALIZAÇÃO: Usando esses locais em quadratura e pesos para ( aqui), podemos aproximar a integral sobre por uma soma dos termos com substituído por e por todo termo multiplicado pelos pesos respectivos . Assim, nossa função de probabilidade deve ser agora$g_r$$w_r$$r=1,...,R$$R=10$$b_j$$R$$g_r$$b_j$$w_r$

Além disso, precisamos levar em consideração a variação da parte aleatória. Li que isso pode ser alcançado substituindo o em nossa função por que , portanto, na função de probabilidade acima, substituímos 's por ' s e não 's.$b_j \sim N(0,\sigma_b^2)$$\eta$$\sigma_j \theta_j$$\theta_j \sim N(0,1)$$\theta$$g$$\beta$

Uma questão computacional que não entendo é como substituir os termos, pois os vetores não terão o mesmo comprimento. Mas provavelmente não entendo isso, porque estou perdendo algo crucial aqui ou entendi mal como esse método funciona.

Não vi como "os vetores não terão o mesmo comprimento", esclareça sua pergunta.

Primeiro, para a integral com dimensão menor que 4, os métodos numéricos diretos, como a quadratura, são mais eficientes que o MCMC. Estudei essas perguntas por um tempo e ficaria feliz em discutir esse problema com você.

Para regressão logística de efeitos mistos, o único Rcódigo explícito que encontrei é do livro do Prof. Demidenko, Modelos Mistos: Teoria e Aplicações , você pode fazer o download do código através da coluna "SOFTWARE E DADOS" na página da web. O logMLEgh()pode ser encontrado em \mixed_models_data.zip\MixedModels\Chapter07. Ele não usou o statmodpacote para obter as quadraturas, mas escreveu sua própria função gauher(). Existem alguns erros menores no código e eu os discuti com o autor, mas ainda é muito útil começar com o código e o livro dele. Eu posso fornecer a versão corrigida, se necessário.

Outra questão é que, se você deseja obter estimativas precisas, optim()não é suficiente, pode ser necessário usar métodos como a pontuação de Fisher, como em glm().