Comment obtenir une valeur de p et une taille d'effet «globales» pour un facteur catégoriel dans un modèle mixte (lme4)?

28

Je voudrais obtenir une valeur de p et une taille d'effet d'une variable catégorielle indépendante (avec plusieurs niveaux) - c'est-à-dire "global" et pas pour chaque niveau séparément, tout comme la sortie normale de lme4dans R. C'est comme la chose que les gens rapportent lorsqu'ils exécutent une ANOVA.

Comment puis-je l'obtenir?

user3288202
la source
Quelles statistiques voulez-vous exactement? Vous pouvez utiliser la anova()fonction pour obtenir une table anova avec des modèles mixtes linéaires comme avec les modèles linéaires.
smillig
J'ai essayé anova () mais cela me donne des valeurs Df, Sum Sq, Mean Sq et F. Je ne vois pas la taille de l'effet et la valeur p. Avez-vous des idées à ce sujet?
user3288202
1
Par taille d'effet, voulez-vous dire quelque chose comme un équivalent à ? En ce qui concerne les valeurs de p, il existe un débat long et substantiel autour de leur estimation et de leur mise en œuvre en . Jetez un œil à la discussion dans cette question pour plus de détails. R2lme4
smillig
Merci pour le lien, Smilig. Est-ce à dire que parce qu'il y a un problème avec le calcul de la valeur de p, la taille d'effet du facteur dans l'ensemble est également un problème?
user3288202
Ce ne sont pas des problèmes directement liés. Cependant, vous devez garder à l'esprit qu'un modèle mixte linéaire ne se comporte pas exactement comme un modèle linéaire sans effets aléatoires, donc une mesure qui peut être appropriée pour le modèle linéaire ne se généralise pas nécessairement aux modèles mixtes.
smillig

Réponses:

48

Les deux concepts que vous mentionnez (valeurs de p et tailles d'effet des modèles mixtes linéaires) ont des problèmes inhérents. En ce qui concerne la taille de l'effet , citant Doug Bates, l'auteur original de lme4,

En supposant que l'on veuille définir une mesure , je pense qu'un argument pourrait être avancé pour traiter la somme résiduelle pénalisée des carrés d'un modèle mixte linéaire de la même manière que nous considérons la somme résiduelle des carrés d'un modèle linéaire. Ou on pourrait utiliser uniquement la somme résiduelle des carrés sans pénalité ou la somme résiduelle minimale des carrés pouvant être obtenue à partir d'un ensemble de termes donné, ce qui correspond à une matrice de précision infinie. Je ne sais pas vraiment. Cela dépend de ce que vous essayez de caractériser.R2

Pour plus d'informations, vous pouvez consulter ce fil , ce fil et ce message . Fondamentalement, le problème est qu'il n'y a pas de méthode convenue pour l'inclusion et la décomposition de la variance des effets aléatoires dans le modèle. Cependant, quelques normes sont utilisées. Si vous jetez un œil à la liste de diffusion Wiki configurée pour / par la liste de diffusion r-sig-mixed-models , deux approches sont répertoriées.

L'une des méthodes suggérées examine la corrélation entre les valeurs ajustées et observées. Cela peut être implémenté dans R comme suggéré par Jarrett Byrnes dans l'un de ces threads:

r2.corr.mer <- function(m) {
  lmfit <-  lm(model.response(model.frame(m)) ~ fitted(m))
  summary(lmfit)$r.squared
}

Par exemple, supposons que nous estimons le modèle mixte linéaire suivant:

set.seed(1)
d <- data.frame(y = rnorm(250), x = rnorm(250), z = rnorm(250),
                g = sample(letters[1:4], 250, replace=T)       )
library(lme4)
summary(fm1 <- lmer(y ~ x + (z | g), data=d))
# Linear mixed model fit by REML ['lmerMod']
# Formula: y ~ x + (z | g)
#    Data: d
# REML criterion at convergence: 744.4
# 
# Scaled residuals: 
#     Min      1Q  Median      3Q     Max 
# -2.7808 -0.6123 -0.0244  0.6330  3.5374 
# 
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr 
#  g        (Intercept) 0.006218 0.07885       
#           z           0.001318 0.03631  -1.00
#  Residual             1.121439 1.05898       
# Number of obs: 250, groups: g, 4
# 
# Fixed effects:
#             Estimate Std. Error t value
# (Intercept)  0.02180    0.07795   0.280
# x            0.04446    0.06980   0.637
# 
# Correlation of Fixed Effects:
#   (Intr)
# x -0.005

Nous pouvons calculer la taille de l'effet en utilisant la fonction définie ci-dessus:

r2.corr.mer(fm1)
# [1] 0.0160841

Ω02

1-var(residuals(fm1))/(var(model.response(model.frame(fm1))))
# [1] 0.01173721  # Usually, it would be even closer to the value above

En ce qui concerne les valeurs p , c'est une question beaucoup plus litigieuse (au moins dans la communauté R / lme4). Voir les discussions dans les questions ici , ici et ici parmi beaucoup d'autres. En référençant à nouveau la page Wiki, il existe quelques approches pour tester des hypothèses sur les effets dans des modèles mixtes linéaires. Répertorié du "pire au meilleur" (selon les auteurs de la page Wiki qui comprend, je crois, Doug Bates ainsi que Ben Bolker qui y contribue beaucoup):

  • Wald Z-tests
  • Pour les LMM équilibrés et imbriqués où df peut être calculé: tests t de Wald
  • Test de rapport de vraisemblance, soit en configurant le modèle de manière à pouvoir isoler / supprimer le paramètre (via anovaou drop1), soit en calculant des profils de vraisemblance
  • MCMC ou intervalles de confiance de bootstrap paramétrique

Ils recommandent l'approche d'échantillonnage Monte Carlo de la chaîne de Markov et énumèrent également un certain nombre de possibilités de mise en œuvre à partir d'approches pseudo et entièrement bayésiennes, énumérées ci-dessous.

Pseudo-bayésien:

  • Échantillonnage post hoc, généralement (1) en supposant des a priori plats et (2) en partant de la MLE, en utilisant éventuellement l'estimation approximative de la variance-covariance pour choisir une distribution candidate
  • Via mcmcsamp(si disponible pour votre problème: c'est-à-dire LMM avec des effets aléatoires simples - pas des GLMM ou des effets aléatoires complexes)
    Via pvals.fncdans le languageRpackage, un wrapper pour mcmcsamp)
  • Dans AD Model Builder, éventuellement via le glmmADMBpackage (utilisez l' mcmc=TRUEoption) ou le R2admbpackage (écrivez votre propre définition de modèle dans AD Model Builder), ou en dehors de R
  • Via la simfonction du armpackage (simule le postérieur uniquement pour les coefficients bêta (à effet fixe)

Approches entièrement bayésiennes:

  • Via le MCMCglmmpackage
  • Utilisation glmmBUGS(d'une interface wrapper / R WinBUGS )
  • Utilisation de JAGS / WinBUGS / OpenBUGS etc., via les packages rjags/ r2jags/ R2WinBUGS/BRugs

À des fins d'illustration pour montrer à quoi cela pourrait ressembler, ci-dessous est une MCMCglmmestimation en utilisant le MCMCglmmpackage qui, vous verrez, donne des résultats similaires au modèle ci-dessus et a une sorte de valeurs p bayésiennes:

library(MCMCglmm)
summary(fm2 <- MCMCglmm(y ~ x, random=~us(z):g, data=d))
# Iterations = 3001:12991
# Thinning interval  = 10
#  Sample size  = 1000 
# 
#  DIC: 697.7438 
# 
#  G-structure:  ~us(z):g
# 
#       post.mean  l-95% CI u-95% CI eff.samp
# z:z.g 0.0004363 1.586e-17 0.001268    397.6
# 
#  R-structure:  ~units
# 
#       post.mean l-95% CI u-95% CI eff.samp
# units    0.9466   0.7926    1.123     1000
# 
#  Location effects: y ~ x 
# 
#             post.mean l-95% CI u-95% CI eff.samp pMCMC
# (Intercept)  -0.04936 -0.17176  0.07502     1000 0.424
# x            -0.07955 -0.19648  0.05811     1000 0.214

J'espère que cela aide quelque peu. Je pense que le meilleur conseil pour quelqu'un commençant avec des modèles mixtes linéaires et essayant de les estimer dans R est de lire les FAQ Wiki d'où la plupart de ces informations ont été tirées. C'est une excellente ressource pour toutes sortes de thèmes d'effets mixtes du basique au avancé et de la modélisation au traçage.

smillig
la source
Merci beaucoup smilig. Je ne peux donc pas signaler la taille de l'effet pour les paramètres globaux.
user3288202
r2
russellpierce
3
+6, d'une clarté impressionnante, complet et entièrement annoté.
gung - Rétablir Monica
1
en outre, vous pouvez jeter un œil au package afex et en particulier à la fonction mixte. voir ici
beginneR
6

En ce qui concerne le calcul de la signification ( valeurs de p ), Luke (2016), Évaluation de la signification dans les modèles linéaires à effets mixtes dans R, indique que la méthode optimale est l'approximation de Kenward-Roger ou Satterthwaite pour les degrés de liberté (disponible en R avec des packages tels que lmerTestou afex).

Abstrait

Les modèles à effets mixtes sont de plus en plus utilisés dans l'analyse des données expérimentales. Cependant, dans le paquet lme4 de R, les normes d'évaluation de l'importance des effets fixes dans ces modèles (c'est-à-dire l'obtention des valeurs de p) sont quelque peu vagues. Il y a de bonnes raisons à cela, mais comme les chercheurs qui utilisent ces modèles sont tenus dans de nombreux cas de déclarer des valeurs de p, une méthode pour évaluer la signification de la sortie du modèle est nécessaire. Cet article présente les résultats de simulations montrant que les deux méthodes les plus courantes pour évaluer la signification, en utilisant des tests de rapport de vraisemblance et en appliquant la distribution z aux valeurs de Wald t à partir de la sortie du modèle (t-as-z), sont quelque peu anti-conservatrices, spécialement pour les échantillons de plus petite taille. Autres méthodes d'évaluation de l'importance,Les résultats de ces simulations suggèrent que les taux d'erreur de type 1 sont les plus proches de 0,05 lorsque les modèles sont ajustés en utilisant REML et les valeurs de p sont dérivées en utilisant les approximations de Kenward-Roger ou Satterthwaite, car ces approximations ont toutes deux produit des taux d'erreur de type 1 acceptables même pour les plus petits échantillons.

(pas d'italique dans l'original)

Pablo Bernabeu
la source
4
+1 Merci d'avoir partagé ce lien. Je vais simplement commenter brièvement que l'approximation de Kenward-Roger est disponible dans le lmerTestpackage.
amibe dit Réintégrer Monica
5

J'utilise le lmerTestpackage. Cela inclut commodément une estimation de la valeur de p dans la anova()sortie pour mes analyses MLM, mais ne donne pas de taille d'effet pour les raisons données dans d'autres articles ici.

Bruna
la source
1
Dans mon cas, je préfère la comparaison par paire en utilisant lsmeans car elle me donne toutes les paires de contrastes, y compris les valeurs de p. Si j'utilise lmerTest, je devrai exécuter le modèle six fois avec différentes lignes de base pour voir toutes les paires de contrastes.
user3288202