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

8

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

logL(β)=i=1nyilogΛ(ηi)+(1yi)log(1Λ(ηi))

onde Λ(ηi)=11+exp(ηi) e ηi=Xiβ 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 é

L=j=1JPr(y1j,...,ynjj|x,bj)f(bj)dbj

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 agoragrwrr=1,...,RR=10bjRgrbjwr

L=j=1Jr=1RPr(y1j,...,ynjj|x,bj=gr)wr

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.bjN(0,σb2)ησjθjθjN(0,1)θgβ

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.

Steve
fonte
Não tenho tempo para dar uma boa olhada agora, mas isso parece ser um bom uso para o MCMC.
shadowtalker
@ssdecontrol obrigado, o MCMC é uma boa alternativa. Mas eu gostaria de aplicar a abordagem clássica.
Steve Steve
O que não é clássico sobre o MCMC para avaliar uma integral?
shadowtalker 22/08
@ssdecontrol bem, não tenho certeza ... Mas todos os livros que verifiquei e os pacotes ordinais R do lme4, não usam o MCMC. Então, eu gostaria de continuar com isso. Pelo menos no começo.
Steve
1
Você já perguntou isso na lista R-sig-ME ([email protected])? Algumas pessoas podem ajudá-lo ainda mais. Além disso: eu recomendaria que você estude o artigo Algoritmos de Quadratura Gaussiana Eficiente de Laplaciano e Adaptativo para Modelos Mistos Lineares Generalizados Multinível de Pinheiro e Chao. Ele contém resultados sobre o desempenho do AGQ, bem como a álgebra linear por trás deles. Se você quiser codificá-los ... bem, prepare-se para algumas decomposições de matriz esparsas sérias. : D
usεr11852

Respostas:

2

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().

Randel
fonte
o livro parece ter informações valiosas sobre o que estou trabalhando. O código em si não diz muito - apenas pedi o livro, então tenho que esperar por isso. A questão dos vetores é que, se substituirmos os termos, os 2 bpelos 10 nós, como poderemos multiplicar as matrizes Ze g? Ou eu entendi completamente errado?
Steve
Sei que devo ir além para obter estimativas precisas, mas esperava entender primeiro o GQ como um primeiro passo.
Steve
Você pode visualizar as duas edições primeiro no Google Livros . No seu código, é uma matriz ? Mas apenas um escalar? Você pode começar pelo modelo de interceptação aleatória primeiro. Se a dimensão dos efeitos aleatórios for dois, você precisará totalmente de 100 nós bidimensionais, 10 nós em cada dimensão. Zn×2bZ = rep(1,n)
Randel
desculpe, mas quanto mais eu estou pensando, mais isso me confunde. Se eu fizer Z=rep(1,n)isso, obteria uma interceptação aleatória para cada linha, o que significa que cada indivíduo é um grupo? No meu exemplo, eu tenho dois grupos, portanto, temos e para fornecer o que precisamos. Não? Z%*%b n×22×1n×1
Steve
Ah, acabei de notar que isso Zé usado para separar a interceptação aleatória de cada cluster, não a matriz de design para o efeito aleatório. Então você está certo, mas deve avaliar a integral e utilizar a quadratura separadamente para cada cluster. Você não precisa Zmais, basta avaliar a integral para cada cluster e, em seguida, somar. A matriz de design para a interceptação aleatória é apenas o vetor de 1s.
Randel