Dans un modèle à plusieurs niveaux, quelles sont les implications pratiques de l'estimation et de la non-estimation des paramètres de corrélation à effet aléatoire?

Dans un modèle à plusieurs niveaux, quelles sont les implications pratiques et liées à l'interprétation de l'estimation et de la non-estimation des paramètres de corrélation à effet aléatoire? La raison pratique de demander ceci est que dans le cadre de lmer dans R, il n'y a pas de méthode implémentée pour estimer les valeurs de p via des techniques MCMC lorsque des estimations sont faites dans le modèle des corrélations entre paramètres.

Par exemple, en regardant cet exemple (parties citées ci-dessous), quelles sont les implications pratiques de M2 ​​par rapport à M3. De toute évidence, dans un cas, P5 ne sera pas estimé et dans l'autre, il le sera.

Des questions

1. Pour des raisons pratiques (le désir d'obtenir une valeur de p grâce aux techniques MCMC), on peut vouloir ajuster un modèle sans corrélation entre les effets aléatoires même si P5 est sensiblement non nul. Si l'on fait cela, puis si l'on estime les valeurs de p via la technique MCMC, les résultats sont-ils interprétables? (Je sais que @Ben Bolker a déjà mentionné que "combiner les tests de signification avec MCMC est un peu incohérent, statistiquement, bien que je comprenne l'envie de le faire (obtenir des intervalles de confiance est plus supportable)" , donc si cela vous fera mieux dormir la nuit, je fais semblant d'avoir dit des intervalles de confiance.)
2. Si l'on ne parvient pas à estimer P5, est-ce la même chose que d'affirmer qu'il est égal à 0?
3. Si P5 est vraiment non nul, de quelle manière les valeurs estimées de P1-P4 sont-elles affectées?
4. Si P5 est vraiment non nul, de quelle manière les estimations d'erreur pour P1-P4 sont-elles affectées?
5. Si P5 est vraiment différent de zéro, de quelles manières les interprétations d'un modèle ne parviennent-elles pas à inclure P5?

Empruntant à la réponse de @Mike Lawrence (ceux qui sont mieux informés que moi sont libres de remplacer cela par une notation complète du modèle, je ne suis pas entièrement convaincu de pouvoir le faire avec une fidélité raisonnable):

M2: V1 ~ (1|V2) + V3 + (0+V3|V2)(estimations P1 - P4)

M3: V1 ~ (1+V3|V2) + V3(estimations P1-P5)

Paramètres pouvant être estimés:

P1 : Une interception globale

P2 : interceptions à effet aléatoire pour V2 (c'est-à-dire pour chaque niveau de V2, écart de l'interception de ce niveau par rapport à l'interception globale)

P3 : Une seule estimation globale de l'effet (pente) de V3

P4 : L'effet de V3 à l'intérieur de chaque niveau de V2 (plus précisément, le degré auquel l'effet V3 à l'intérieur d'un niveau donné s'écarte de l'effet global de V3), tout en imposant une corrélation nulle entre les écarts d'interception et les écarts d'effet V3 à travers les niveaux de V2.

P5 : La corrélation entre les écarts d'interception et les écarts V3 à travers les niveaux de V2

Les réponses dérivées d'une simulation suffisamment large et large ainsi que le code d'accompagnement dans R utilisant lmer seraient acceptables.

Considérez les données de l'étude du sommeil, incluses dans lme4. Bates en parle dans son livre en ligne sur lme4. Dans le chapitre 3, il considère deux modèles pour les données.

et

L'étude a porté sur 18 sujets, étudiés sur une période de 10 jours sans sommeil. Les temps de réaction ont été calculés au départ et les jours suivants. Il y a un effet clair entre le temps de réaction et la durée de la privation de sommeil. Il existe également des différences significatives entre les sujets. Le modèle A permet la possibilité d'une interaction entre l'interception aléatoire et les effets de pente: imaginons, par exemple, que les personnes ayant de faibles temps de réaction souffrent plus fortement des effets de la privation de sommeil. Cela impliquerait une corrélation positive dans les effets aléatoires.

Dans l'exemple de Bates, il n'y avait aucune corrélation apparente à partir du tracé du réseau et aucune différence significative entre les modèles. Cependant, pour enquêter sur la question posée ci-dessus, j'ai décidé de prendre les valeurs ajustées de l'étude du sommeil, d'augmenter la corrélation et d'examiner les performances des deux modèles.

Comme vous pouvez le voir sur l'image, les longs temps de réaction sont associés à une plus grande perte de performances. La corrélation utilisée pour la simulation était de 0,58

J'ai simulé 1000 échantillons, en utilisant la méthode de simulation dans lme4, sur la base des valeurs ajustées de mes données artificielles. J'ai adapté M0 et Ma à chacun et j'ai regardé les résultats. L'ensemble de données d'origine comportait 180 observations (10 pour chacun des 18 sujets) et les données simulées ont la même structure.

L'essentiel est qu'il y a très peu de différence.

1. Les paramètres fixes ont exactement les mêmes valeurs sous les deux modèles.
2. Les effets aléatoires sont légèrement différents. Il existe 18 effets aléatoires d'interception et 18 de pente pour chaque échantillon simulé. Pour chaque échantillon, ces effets sont forcés de s'ajouter à 0, ce qui signifie que la différence moyenne entre les deux modèles est (artificiellement) 0. Mais les variances et les covariances diffèrent. La covariance médiane sous MA était de 104, contre 84 sous M0 (valeur réelle, 112). Les variances des pentes et des intersections étaient plus grandes sous M0 que MA, probablement pour obtenir la marge de manœuvre supplémentaire nécessaire en l'absence d'un paramètre de covariance libre.
3. La méthode ANOVA pour lmer donne une statistique F pour comparer le modèle Slope à un modèle avec seulement une interception aléatoire (aucun effet dû à la privation de sommeil). De toute évidence, cette valeur était très élevée dans les deux modèles, mais elle était généralement (mais pas toujours) plus élevée sous MA (moyenne 62 vs moyenne de 55).
4. La covariance et la variance des effets fixes sont différentes.
5. Environ la moitié du temps, il sait que MA est correct. La valeur p médiane pour comparer M0 à MA est de 0,0442. Malgré la présence d'une corrélation significative et de 180 observations équilibrées, le modèle correct ne serait choisi qu'environ la moitié du temps.
6. Les valeurs prévues diffèrent sous les deux modèles, mais très légèrement. La différence moyenne entre les prédictions est de 0, avec sd de 2,7. Le sd des valeurs prédites elles-mêmes est de 60,9

Alors pourquoi cela se produit-il? @gung a deviné, raisonnablement, que le fait de ne pas inclure la possibilité d'une corrélation force les effets aléatoires à ne pas être corrélés. Peut-être que cela devrait; mais dans cette implémentation, les effets aléatoires peuvent être corrélés, ce qui signifie que les données sont capables de tirer les paramètres dans la bonne direction, quel que soit le modèle. La fausseté du mauvais modèle apparaît dans la probabilité, c'est pourquoi vous pouvez (parfois) distinguer les deux modèles à ce niveau. Le modèle à effets mixtes ajuste essentiellement les régressions linéaires à chaque sujet, influencées par ce que le modèle pense qu'elles devraient être. Le mauvais modèle force l'ajustement de valeurs moins plausibles que celles obtenues avec le bon modèle. Mais les paramètres, en fin de compte, sont régis par l'ajustement aux données réelles.

Voici mon code quelque peu maladroit. L'idée était d'ajuster les données de l'étude sur le sommeil, puis de créer un ensemble de données simulées avec les mêmes paramètres, mais une plus grande corrélation pour les effets aléatoires. Cet ensemble de données a été alimenté dans simulate.lmer () pour simuler 1000 échantillons, chacun étant ajusté dans les deux sens. Une fois que j'avais apparié des objets ajustés, je pouvais retirer différentes caractéristiques de l'ajustement et les comparer, en utilisant des tests t, ou autre chose.

    # Fit a model to the sleep study data, allowing non-zero correlation
fm01 <- lmer(Reaction ~ 1 + Days +(1+Days|Subject), data=sleepstudy, REML=FALSE)
# Now use this to build a similar data set with a correlation = 0.9
# Here is the covariance function for the random effects
# The variances come from the sleep study. The covariance is chosen to give a larger correlation
sigma.Subjects <- matrix(c(565.5,122,122,32.68),2,2) 
# Simulate 18 pairs of random effects
ranef.sim <- mvrnorm(18,mu=c(0,0),Sigma=sigma.Subjects)
# Pull out the pattern of days and subjects.
XXM <- model.frame(fm01) 
n <- nrow(XXM) # Sample size
# Add an intercept to the model matrix.
XX.f <- cbind(rep(1,n),XXM[,2])
# Calculate the fixed effects, using the parameters from the sleep study. 
yhat <- XX.f %*%  fixef(fm01 )
# Simulate a random intercept for each subject
intercept.r <- rep(ranef.sim[,1], each=10) 
# Now build the random slopes
slope.r <- XXM[,2]*rep(ranef.sim[,2],each=10)
# Add the slopes to the random intercepts and fixed effects
yhat2 <- yhat+intercept.r+slope.r
# And finally, add some noise, using the variance from the sleep study
y <- yhat2 + rnorm(n,mean=0,sd=sigma(fm01))
# Here is new "sleep study" data, with a stronger correlation.
new.data <- data.frame(Reaction=y,Days=XXM$Days,Subject=XXM$Subject)
# Fit the new data with its correct model
fm.sim <- lmer(Reaction ~ 1 + Days +(1+Days|Subject), data=new.data, REML=FALSE)
# Have a look at it
xyplot(Reaction ~ Days | Subject, data=new.data, layout=c(6,3), type=c("p","r"))
# Now simulate 1000 new data sets like new.data and fit each one
# using the right model and zero correlation model.
# For each simulation, output a list containing the fit from each and
# the ANOVA comparing them.
n.sim <- 1000
    sim.data <- vector(mode="list",)
    tempReaction <- simulate(fm.sim, nsim=n.sim)
    tempdata <- model.frame(fm.sim)
    for (i in 1:n.sim){
        tempdata$Reaction <- tempReaction[,i]
			output0 <- lmer(Reaction ~ 1 + Days +(1|Subject)+(0+Days|Subject), data = tempdata, REML=FALSE)
			output1 <- lmer(Reaction ~ 1 + Days +(Days|Subject), data=tempdata, REML=FALSE)
			temp <- anova(output0,output1)
			pval <- temp$`Pr(>Chisq)`[2]
        sim.data[[i]] <- list(model0=output0,modelA=output1, pvalue=pval)
    }

Placidia a déjà fourni une réponse complète en utilisant des données simulées basées sur l' sleepstudyensemble de données. Voici une autre réponse (moins rigoureuse) qui utilise également les sleepstudydonnées.

Nous voyons que l'on peut affecter la corrélation estimée entre l'ordonnée à l'origine aléatoire et la pente aléatoire en «décalant» la variable de prédiction aléatoire. Regardez les résultats des modèles fm1et fm2ci - dessous:

library(lmer)

#Fit Models
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
k <- 3 # Shift "Days" by an arbitrary amount
fm2 <- lmer(Reaction ~ I(Days + k) + (I(Days + k)| Subject), sleepstudy)

fm1 # Model Output
# Linear mixed model fit by REML ['lmerMod']
# Formula: Reaction ~ Days + (Days | Subject)
# Data: sleepstudy
# REML criterion at convergence: 1743.628
# Random effects:
#   Groups   Name        Std.Dev. Corr
# Subject  (Intercept) 24.740       
# Days         5.922   0.07
# Residual             25.592       
# Number of obs: 180, groups:  Subject, 18
# Fixed Effects:
#   (Intercept)         Days  
# 251.41        10.47

fm2 # Model Output
# Linear mixed model fit by REML ['lmerMod']
# Formula: Reaction ~ I(Days + k) + (I(Days + k) | Subject)
# Data: sleepstudy
# REML criterion at convergence: 1743.628
# Random effects:
#   Groups   Name        Std.Dev. Corr 
# Subject  (Intercept) 29.498        
# I(Days + k)  5.922   -0.55
# Residual             25.592        
# Number of obs: 180, groups:  Subject, 18
# Fixed Effects:
#   (Intercept)  I(Days + k)  
# 220.00        10.47

# Random effects from both models
cbind(ranef(fm1)$Subject,ranef(fm2)$Subject)
# (Intercept)        Days (Intercept) I(Days + k)
# 308   2.2585654   9.1989719 -25.3383538   9.1989727
# 309 -40.3985769  -8.6197032 -14.5394628  -8.6197043
# 310 -38.9602458  -5.4488799 -22.6136027  -5.4488807
# 330  23.6904985  -4.8143313  38.1334933  -4.8143315
# 331  22.2602027  -3.0698946  31.4698868  -3.0698946
# 332   9.0395259  -0.2721707   9.8560377  -0.2721706
# 333  16.8404311  -0.2236244  17.5113040  -0.2236243
# 334  -7.2325792   1.0745761 -10.4563076   1.0745761
# 335  -0.3336958 -10.7521591  31.9227854 -10.7521600
# 337  34.8903508   8.6282840   9.0054946   8.6282850
# 349 -25.2101104   1.1734142 -28.7303527   1.1734141
# 350 -13.0699567   6.6142050 -32.9125736   6.6142054
# 351   4.5778352  -3.0152572  13.6236077  -3.0152574
# 352  20.8635924   3.5360133  10.2555505   3.5360138
# 369   3.2754530   0.8722166   0.6588028   0.8722167
# 370 -25.6128694   4.8224646 -40.0802641   4.8224648
# 371   0.8070397  -0.9881551   3.7715053  -0.9881552
# 372  12.3145393   1.2840297   8.4624492   1.2840300

De la sortie du modèle, nous voyons que la corrélation de variance aléatoire a changé. Cependant, les pentes (fixes et aléatoires) sont restées les mêmes, tout comme l'estimation de la variance résiduelle. Les estimations des intersections (fixes et aléatoires) ont changé en réponse à la variable décalée.

Décorrélation covariance interception pente aléatoire pour LMM est discuté dans les notes de cours du Dr Jack Weiss ici . Weiss note que la réduction de la corrélation de variance de cette manière peut parfois aider, entre autres, à la convergence des modèles.

L'exemple ci-dessus fait varier la corrélation aléatoire (paramètre "P5"). Concernant partiellement le Q3 de l'OP, nous voyons dans la sortie ci-dessus que:

#   Parameter           Status
=================================
P1  Fixed Intercept     Affected
P2  Random Intercepts   Affected
P3  Fixed Slope         Not Affected
P4  Random Slopes       Not Affected
P5  Random Correlation  Affected