Probabilité et estimations des effets mixtes Régression logistique

8

Simulons d'abord quelques données pour une régression logistique avec des parties fixes et aléatoires:

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 nous voulions simplement ajuster une régression logistique sans parties aléatoires, nous pourrions utiliser la glmfonction:

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

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

Ou construire notre propre fonction de log-vraisemblance

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

Λ(ηi)=11+exp(ηi) et ηi=Xiβ et utilisez optim()pour estimer les paramètres qui le maximisent, comme ci-dessous exemple de code:

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

ce qui bien sûr donne les mêmes estimations et maximise la log-vraisemblance pour la même valeur. Pour les effets mixtes, nous voudrions quelque chose comme

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

Mais comment pouvons-nous faire de même avec notre propre fonction? Étant donné que la probabilité est

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

et l'intégrale n'a pas d'expression de forme fermée, nous devons utiliser l'intégration numérique comme la quadrature gaussienne. Nous pouvons utiliser le package statmodpour obtenir des quadratures, disons 10

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

MISE À JOUR: En utilisant ces emplacements de quadrature et les poids pour les ( ici), nous pouvons approximer l'intégrale sur par une somme des termes avec substitué à et l'ensemble terme multiplié par les poids respectifs . Ainsi, notre fonction de vraisemblance devrait être maintenantgrwrr=1,...,RR=10bjRgrbjwr

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

De plus, nous devons tenir compte de la variance de la partie aléatoire, j'ai lu que cela peut être réalisé en remplaçant le dans notre fonction par où , donc dans la fonction de vraisemblance ci-dessus, nous remplaçons en fait les par des et non des .bjN(0,σb2)ησjθjθjN(0,1)θgβ

Un problème de calcul que je ne comprends pas est de savoir comment remplacer les termes car les vecteurs ne seront pas de la même longueur. Mais je ne comprends probablement pas cela, car je manque quelque chose de crucial ici, ou j'ai mal compris comment cette méthode fonctionne.

Steve
la source
Je n'ai pas le temps de bien regarder maintenant, mais cela semble être une bonne utilisation pour MCMC.
shadowtalker
@ssdecontrol merci, MCMC est une bonne alternative. Mais je voudrais appliquer l'approche classique.
Steve
Qu'est-ce qui n'est pas classique dans MCMC pour évaluer une intégrale?
shadowtalker du
@ssdecontrol bien, je ne suis pas sûr ... Mais tous les livres que j'ai vérifiés et les paquets lme4, ordinaux R, n'utilisent pas MCMC. Je voudrais donc m'en tenir à cela. Au moins au début.
Steve
1
Avez-vous posé cette question dans la liste R-sig-ME ([email protected])? Certaines personnes pourraient vous aider davantage. De plus: je vous exhorte fortement à étudier le papier Algorithmes efficaces en laplacien et en quadrature gaussienne adaptative pour les modèles mixtes linéaires généralisés à plusieurs niveaux de Pinheiro et Chao. Il contient des résultats concernant les performances AGQ ainsi que l'algèbre linéaire derrière eux. Si vous voulez les coder ... eh bien, préparez-vous à de sérieuses décompositions matricielles clairsemées. : D
usεr11852

Réponses:

2

Je n'ai pas vu comment "les vecteurs ne seront pas de la même longueur", veuillez clarifier votre question.

Tout d'abord, pour l'intégrale de dimension inférieure à 4, les méthodes numériques directes comme la quadrature sont plus efficaces que MCMC. J'ai étudié ces questions pendant un certain temps et je serais heureux de discuter de ce problème avec vous.

Pour la régression logistique à effets mixtes, le seul Rcode explicite que j'ai trouvé provient du livre du professeur Demidenko, Mixed Models: Theory and Applications , vous pouvez télécharger le code via la colonne "LOGICIELS ET DONNÉES" sur la page Web. Le logMLEgh()peut être trouvé dans \mixed_models_data.zip\MixedModels\Chapter07. Il n'a pas utilisé le statmodpackage pour obtenir les quadratures, mais a écrit sa propre fonction gauher(). Il y a quelques erreurs mineures dans le code et j'en ai discuté avec l'auteur, mais il est toujours très utile de partir de son code et de son livre. Je peux fournir la version corrigée si nécessaire.

Un autre problème est que, si vous voulez obtenir des estimations précises, ce optim()n'est pas suffisant, vous devrez peut-être utiliser des méthodes comme la notation de Fisher, comme dans glm().

Randel
la source
le livre semble avoir de riches informations sur ce sur quoi je travaille. Le code lui-même ne dit pas grand-chose - je viens de commander le livre, donc je dois attendre cela. La chose à propos des vecteurs est que si nous substituons les termes, les 2 baux 10 nœuds, comment pourrons-nous multiplier les matrices Zet g? Ou je me trompe complètement?
Steve
Je sais que je devrais aller plus loin pour obtenir des estimations précises, mais j'espérais d'abord comprendre le GQ dans un premier temps.
Steve
Vous pouvez d'abord prévisualiser les deux éditions sur Google Livres . Dans votre code, est une matrice ? Mais juste un scalaire? Vous pouvez commencer par le modèle d'interception aléatoire . Si la dimension des effets aléatoires est de deux, vous avez alors besoin de 100 nœuds à deux dimensions, 10 nœuds sur chaque dimension. Zn×2bZ = rep(1,n)
Randel
désolé mais plus j'y pense, plus ça m'embrouille. Si je fais Z=rep(1,n)alors, j'obtiendrais une interception aléatoire pour chaque ligne, ce qui signifie que chaque individu est un groupe? Dans mon exemple, j'ai deux groupes, nous avons donc et pour donner les nous avons besoin. Non? Z%*%b n×22×1n×1
Steve
Oh, je viens de remarquer que Zc'est utilisé pour séparer l'interception aléatoire pour chaque cluster, pas la matrice de conception pour l'effet aléatoire. Alors vous avez raison, mais vous devez évaluer l'intégrale et utiliser la quadrature séparément pour chaque cluster. Vous n'en avez plus besoin Z, évaluez simplement l'intégrale pour chaque cluster, puis additionnez-les. La matrice de conception pour l'ordonnée à l'origine aléatoire n'est que le vecteur de 1s.
Randel