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 glm
fonction:
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
où et
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
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 statmod
pour 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 maintenant
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 .
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.
Réponses:
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
R
code 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. LelogMLEgh()
peut être trouvé dans\mixed_models_data.zip\MixedModels\Chapter07
. Il n'a pas utilisé lestatmod
package pour obtenir les quadratures, mais a écrit sa propre fonctiongauher()
. 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 dansglm()
.la source
b
aux 10 nœuds, comment pourrons-nous multiplier les matricesZ
etg
? Ou je me trompe complètement?Z = rep(1,n)
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
Z
c'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 besoinZ
, é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.