Comment obtenir la valeur p (vérification de la signification) d'un effet dans un modèle mixte lme4?

56

J'utilise lme4 in R pour s'adapter au modèle mixte

lmer(value~status+(1|experiment)))

où la valeur est continue, le statut et l'expérience sont des facteurs, et je reçois

Linear mixed model fit by REML 
Formula: value ~ status + (1 | experiment) 
  AIC   BIC logLik deviance REMLdev
 29.1 46.98 -9.548    5.911    19.1
Random effects:
 Groups     Name        Variance Std.Dev.
 experiment (Intercept) 0.065526 0.25598 
 Residual               0.053029 0.23028 
Number of obs: 264, groups: experiment, 10

Fixed effects:
            Estimate Std. Error t value
(Intercept)  2.78004    0.08448   32.91
statusD      0.20493    0.03389    6.05
statusR      0.88690    0.03583   24.76

Correlation of Fixed Effects:
        (Intr) statsD
statusD -0.204       
statusR -0.193  0.476

Comment puis-je savoir que l'effet du statut est significatif? R ne rapporte que des valeurs et non des valeurs p .tp

ECII
la source
1
tH0:β=βnullF

Réponses:

61

Il y a beaucoup d'informations sur ce sujet dans la FAQ GLMM . Cependant, dans votre cas particulier, je suggérerais d'utiliser

library(nlme)
m1 <- lme(value~status,random=~1|experiment,data=mydata)
anova(m1)

parce que vous n’avez pas besoin de tout ce que lmervous proposez (vitesse supérieure, traitement des effets aléatoires croisés, GLMM, etc.). lmedevrait vous donner exactement les mêmes estimations de coefficient et de variance, mais calculera également les valeurs de df et de p pour vous (ce qui est logique dans une conception "classique" telle que vous semblez avoir). Vous pouvez également envisager le terme aléatoire ~status|experiment(en permettant une variation des effets de statut d’un bloc à l’autre ou en incluant de manière équivalente une interaction statut par expérience). Les affiches ci-dessus indiquent également que vos tstatistiques sont si grandes que votre p-value sera certainement <0,05, mais je peux imaginer que vous aimeriez de "vraies" p-valeurs.

Ben Bolker
la source
3
Je ne sais pas pour cette réponse. lmerpourrait tout aussi bien signaler les mêmes types de valeurs p, mais pas pour des raisons valables. J'imagine que c'est le commentaire qu'il y a de "vraies" valeurs-p ici qui me dérange. Vous pourriez faire valoir que vous pouvez trouver un seuil possible et que tout seuil raisonnable est dépassé. Mais vous ne pouvez pas dire qu'il existe une vraie valeur p.
John
11
Pour un plan classique (équilibré, imbriqué, etc.), je pense que je peux effectivement affirmer qu'il existe une valeur réelle, c'est-à-dire une probabilité d'obtenir une estimation de bêta d'une magnitude observée ou supérieure si l'hypothèse nulle (bêta = 0) were false ... Je pense que lme4 ne fournit pas ce dénominateur df, car il est plus difficile de détecter en général à partir d'une structure de modèle lme4 lorsque le modèle spécifié est celui où une méthode heuristique permettant de calculer un dénominateur classique df fonctionnerait ...
Ben Bolker
essayez summary(m1)plutôt (je l'utilise avec nlme package)
jena
36

Vous pouvez utiliser le paquet lmerTest . Vous venez de l'installer / charger et les modèles Lmer sont étendus. Donc par exemple

library(lmerTest)
lmm <- lmer(value~status+(1|experiment)))
summary(lmm)
anova(lmm)

vous donnerait des résultats avec des valeurs de p. Si les valeurs p sont la bonne indication, c'est un peu contesté, mais si vous voulez les avoir, c'est le moyen de les obtenir.

pbx101
la source
28

Si vous pouvez gérer l'abandon des valeurs p ( et vous devriez le faire ), vous pouvez calculer un rapport de vraisemblance qui représenterait le poids de la preuve de l'effet du statut via:

#compute a model where the effect of status is estimated
unrestricted_fit = lmer(
    formula = value ~ (1|experiment) + status
    , REML = F #because we want to compare models on likelihood
)
#next, compute a model where the effect of status is not estimated
restricted_fit = lmer(
    formula = value ~ (1|experiment)
    , REML = F #because we want to compare models on likelihood
)
#compute the AIC-corrected log-base-2 likelihood ratio (a.k.a. "bits" of evidence)
(AIC(restricted_fit)-AIC(unrestricted_fit))*log2(exp(1))
Mike Lawrence
la source
16
Notez que les ratios de vraisemblance sont asymptotiques, c'est-à-dire qu'ils ne tiennent pas compte de l'incertitude dans l'estimation de la variance résiduelle ...
Ben Bolker
5
Je suis intéressé par votre dernière ligne. Quelle est l'interprétation du résultat? Y a-t-il des sources sur lesquelles je peux jeter un coup d'œil?
mguzmann
13

Le problème est que le calcul des valeurs p pour ces modèles n’est pas trivial. Voir la discussion ici afin que les auteurs du lme4package aient délibérément choisi de ne pas inclure les valeurs p dans la sortie. Vous pouvez trouver une méthode pour les calculer, mais elles ne seront pas nécessairement correctes.

Michelle
la source
9

Considérez ce que vous demandez. Si vous voulez simplement savoir si la valeur p globale pour l'effet de status transmet une valeur limite quelconque, telle que 0,05, alors c'est facile. Tout d'abord, vous voulez connaître l'effet global. Vous pouvez l'obtenir de anova.

m <- lmer(...) #just run your lmer command but save the model
anova(m)

Maintenant, vous avez une valeur F. Vous pouvez prendre cela et le chercher dans certaines tables F. Il suffit de choisir le plus bas denom possible. degrés de liberté. Le seuil sera d'environ 20. Votre F peut être plus grand que cela, mais je peux me tromper. Même si ce n'est pas le cas, examinez ici le nombre de degrés de liberté d'un calcul ANOVA classique en utilisant le nombre d'expériences que vous avez. Si vous tenez compte de cette valeur, vous êtes réduit à environ 5 pour une coupure. Maintenant, vous le passez facilement dans votre étude. Le «vrai» df de votre modèle sera quelque chose de plus élevé que cela, car vous modélisez chaque point de données par opposition aux valeurs agrégées qu'une ANOVA modéliserait.

Si vous voulez réellement une valeur p exacte, il n'y a rien de tel sauf si vous êtes prêt à faire une déclaration théorique à ce sujet. Si vous lisez Pinheiro & Bates (2001, et peut-être quelques autres livres sur le sujet ... voir d'autres liens dans ces réponses) et que vous vous en sortez avec un argument pour une df spécifique, vous pouvez l'utiliser. Mais vous ne recherchez pas réellement une p-valeur exacte de toute façon. Je mentionne cela parce que vous ne devez donc pas déclarer une valeur p exacte, mais simplement que votre seuil est dépassé.

Vous devriez vraiment considérer la réponse de Mike Lawrence, car l'idée de ne conserver qu'un point de passage pour les valeurs p en tant qu'information finale et la plus importante à extraire de vos données est généralement erronée (mais peut-être pas dans votre cas puisque nous ne le faisons pas). t vraiment avoir assez d’informations pour savoir). Mike utilise une version familière du calcul de la LR qui est intéressante, mais il peut être difficile de trouver beaucoup de documentation à ce sujet. Si vous examinez la sélection et l'interprétation des modèles à l'aide de l'AIC, cela peut vous plaire.

John
la source
9

Edit: Cette méthode n'est plus prise en charge dans les nouvelles versions de lme4. Utilisez le paquet lmerTest comme suggéré dans cette réponse par pbx101 .

L'auteur de lme4 a publié un article sur la liste R expliquant pourquoi les valeurs p ne sont pas affichées. Il suggère d'utiliser des exemples MCMC à la place, à l'aide du fichier pvals.fnc du paquet languageR:

library("lme4")
library("languageR")
model=lmer(...)
pvals.fnc(model)

Voir http://www2.hawaii.edu/~kdrager/MixedEffectsModels.pdf pour un exemple et des détails.

Jeff
la source
3
lme4 ne supporte plus cela. Cet article pourrait être mis à jour pour éviter aux personnes de le découvrir comme je viens de le faire.
timothy.s.lau
5

Voulez-vous savoir si l'effet combiné de statusa un effet significatif sur value? Si c'est le cas, vous pouvez utiliser la Anovafonction dans le carpackage (à ne pas confondre avec la anovafonction dans la base R).

dat <- data.frame(
  experiment = sample(c("A","B","C","D"), 264, replace=TRUE), 
  status = sample(c("D","R","A"), 264, replace=TRUE), 
  value = runif(264)   
)
require(lme4)
(fm <- lmer(value~status+(1|experiment), data=dat))

require(car)
Anova(fm)

Regardez ?Anovaaprès le chargement du carpaquet.

Smillig
la source
Avez-vous une idée de la façon dont vous car::Anova()évitez les problèmes épineux entourant le calcul des valeurs p que Michelle lie?
Mike Lawrence
Je ne le fais pas, mais je suppose que cela évite les problèmes épineux en les ignorant! Après avoir relu le message original, j’ai le sentiment que j’ai peut-être mal compris la question. Si le PO veut des valeurs p exactes pour les paramètres des effets fixes, il est en difficulté. Mais si le PO veut juste savoir s’ils sont significatifs, je pense que les valeurs t sont plus grandes que toute incertitude sur la manière dont la valeur p exacte serait calculée. (En d'autres termes, ils sont significatifs.)
smillig
1
Je pense que c’était vraiment une bonne idée de rediriger vers un calcul ANOVA pour connaître l’effet global des statistiques, mais je ne suis pas sûr que l’affinement des valeurs p soit bon. La anovacommande régulière vous donnera des F.
John
Je pense que c'est un peu plus collant qu'apparent. Exécuter ANOVA est valable lorsque vous voulez minimiser la variance, mais d'après le libellé de la question, je pense que OP veut établir l'effet marginal des variables, à savoir les coefficients de test par rapport à zéro.
Firebug
0

La fonction pvals.fncn'est plus supportée par lme4. En utilisant le paquetage lmerTest, il est possible d’utiliser une autre méthode pour calculer la valeur p, telle que les approximations de Kenward-Roger

model=lmer(value~status+1|experiment)
anova(model, ddf="Kenward-Roger")
Stefano Cacciatore
la source
0

Le simple chargement du paquet afex imprimera les valeurs p dans la sortie de la fonction lmer à partir du paquet lme4 (vous n'avez pas besoin d'utiliser l'afex; il suffit de le charger):

library(lme4)  #for mixed model
library(afex)  #for p-values

Cela ajoutera automatiquement une colonne de valeur p à la sortie du lmer (votre modèle) pour les effets fixes.

Maryam Nasseri
la source