Dans quelle mesure les intervalles de confiance des objets lmer sont-ils fiables dans le paquet d'effets?

36

Effectspackage fournit un moyen très rapide et pratique pour tracer les résultats de modèle à effets mixtes linéaires obtenus par lme4package . leeffect fonction calcule très rapidement les intervalles de confiance (IC), mais dans quelle mesure ces intervalles de confiance sont-ils fiables?

Par exemple:

library(lme4)
library(effects)
library(ggplot)

data(Pastes)

fm1  <- lmer(strength ~ batch + (1 | cask), Pastes)
effs <- as.data.frame(effect(c("batch"), fm1))
ggplot(effs, aes(x = batch, y = fit, ymin = lower, ymax = upper)) + 
  geom_rect(xmax = Inf, xmin = -Inf, ymin = effs[effs$batch == "A", "lower"],
        ymax = effs[effs$batch == "A", "upper"], alpha = 0.5, fill = "grey") +
  geom_errorbar(width = 0.2) + geom_point() + theme_bw()

entrez la description de l'image ici

Selon les CI calculés avec le effectspackage, le lot "E" ne chevauche pas le lot "A".

Si j'essaie la même chose en utilisant confint.merModfunction et la méthode par défaut:

a <- fixef(fm1)
b <- confint(fm1)
# Computing profile confidence intervals ...
# There were 26 warnings (use warnings() to see them)

b <- data.frame(b)
b <- b[-1:-2,]

b1 <- b[[1]]
b2 <- b[[2]]

dt <- data.frame(fit   = c(a[1],  a[1] + a[2:length(a)]), 
                 lower = c(b1[1],  b1[1] + b1[2:length(b1)]), 
                 upper = c(b2[1],  b2[1] + b2[2:length(b2)]) )
dt$batch <- LETTERS[1:nrow(dt)]

ggplot(dt, aes(x = batch, y = fit, ymin = lower, ymax = upper)) +
  geom_rect(xmax = Inf, xmin = -Inf, ymin = dt[dt$batch == "A", "lower"], 
        ymax = dt[dt$batch == "A", "upper"], alpha = 0.5, fill = "grey") + 
  geom_errorbar(width = 0.2) + geom_point() + theme_bw()

entrez la description de l'image ici

Je constate que tous les IC se chevauchent. Je reçois également des avertissements indiquant que la fonction n'a pas pu calculer les CI dignes de confiance. Cet exemple, ainsi que mon ensemble de données actuel, me fait penser que le effectspackage prend des raccourcis dans le calcul de CI qui pourraient ne pas être entièrement approuvés par les statisticiens. Quelle est la fiabilité des CI renvoyés par la effectfonction deeffects package pour lmerobjets ?

Qu'ai-je essayé: En regardant dans le code source, j'ai remarqué que la effectfonction s'appuie sur la Effect.merModfonction, qui à son tour dirige à la Effect.merfonction, qui ressemble à ceci:

effects:::Effect.mer
function (focal.predictors, mod, ...) 
{
    result <- Effect(focal.predictors, mer.to.glm(mod), ...)
    result$formula <- as.formula(formula(mod))
    result
}
<environment: namespace:effects>

mer.to.glmLa fonction semble calculer la matrice variance-covariable à partir de l' lmerobjet:

effects:::mer.to.glm

function (mod) 
{
...
mod2$vcov <- as.matrix(vcov(mod))
...
mod2
}

Ceci, à son tour, est probablement utilisé en Effect.defaultfonction pour calculer les CI (j'ai peut-être mal compris cette partie):

effects:::Effect.default
...
     z <- qnorm(1 - (1 - confidence.level)/2)
        V <- vcov.(mod)
        eff.vcov <- mod.matrix %*% V %*% t(mod.matrix)
        rownames(eff.vcov) <- colnames(eff.vcov) <- NULL
        var <- diag(eff.vcov)
        result$vcov <- eff.vcov
        result$se <- sqrt(var)
        result$lower <- effect - z * result$se
        result$upper <- effect + z * result$se
...

Je ne connais pas suffisamment les LMM pour pouvoir juger s’il s’agit d’une bonne approche, mais compte tenu de la discussion sur le calcul de l’intervalle de confiance des LMM, cette approche semble étrangement simple.

Mikko
la source
1
Lorsque vous avez de longues lignes de code, j'apprécierais beaucoup si vous les divisez en plusieurs lignes afin que nous n'ayons pas à faire défiler l'écran pour tout voir.
vendredi
1
@rvl Le code devrait être plus lisible maintenant.
Mikko

Réponses:

52

Tous les résultats sont essentiellement les mêmes ( pour cet exemple particulier ). Quelques différences théoriques sont:

  • comme @rvl le fait remarquer, votre reconstruction des CI sans tenir compte de la covariance entre paramètres est tout simplement fausse (désolé)
  • intervalles de confiance pour les paramètres peuvent être basés sur des intervalles de confiance Wald ( en supposant une surface log-vraisemblance quadratique): lsmeans, effects, confint(.,method="Wald"); À l'exception de lsmeans, ces méthodes ignorent les effets de taille finie ("degrés de liberté"), mais dans ce cas, cela ne fait presque aucune différence ( df=40pratiquement impossible à distinguer de l'infini df)
  • ... ou sur les intervalles de confiance du profil (méthode par défaut; ignore les effets de taille finie mais autorise les surfaces non quadratiques)
  • ... ou sur amorçage paramétrique (l'étalon-or - suppose que le modèle est correct [les réponses sont normales, les effets aléatoires sont normalement distribués, les données sont conditionnellement indépendantes, etc.], mais ne pose que peu d'hypothèses)

Je pense que toutes ces approches sont raisonnables (certaines sont plus approximatives que d’autres), mais dans ce cas, cela ne fait guère de différence que vous utilisez. Si cela vous intéresse, essayez plusieurs méthodes contrastées sur vos données, ou sur des données simulées qui ressemblent aux vôtres, et voyez ce qui se passe ...

(PS: Je ne mettrais pas trop de poids sur le fait que les intervalles de confiance Aet Ene se chevauchent pas , vous auriez à faire une procédure de comparaison par paire appropriée pour tirer des conclusions fiables sur les différences entre ce. Particulier deux estimations. ..)

IC à 95%:

entrez la description de l'image ici

Code de comparaison:

library(lme4)
fm2 <- lmer(strength ~ batch - 1 + (1 | cask), Pastes)
c0 <- confint(fm2,method="Wald")
c1 <- confint(fm2)
c2 <- confint(fm2,method="boot")
library(effects)
library(lsmeans)
c3 <- with(effect("batch",fm2),cbind(lower,upper))
c4 <- with(summary(lsmeans(fm2,spec="batch")),cbind(lower.CL,upper.CL))
tmpf <- function(method,val) {
    data.frame(method=method,
               v=LETTERS[1:10],
               setNames(as.data.frame(tail(val,10)),
                        c("lwr","upr")))
}
library(ggplot2); theme_set(theme_bw())
allCI <- rbind(tmpf("lme4_wald",c0),
      tmpf("lme4_prof",c1),
      tmpf("lme4_boot",c2),
      tmpf("effects",c3),
               tmpf("lsmeans",c4))
ggplot(allCI,aes(v,ymin=lwr,ymax=upr,colour=method))+
    geom_linerange(position=position_dodge(width=0.8))

ggsave("pastes_confint.png",width=10)
Ben Bolker
la source
2
J'accepte cette réponse car elle va droit au but et donne une belle comparaison entre différentes méthodes. Cependant, consultez l'excellente réponse de rlv pour plus d'informations.
Mikko
Merci pour la réponse excellente et très utile. Si j'ai bien compris, on ne peut pas utiliser les CI pour comparer des groupes / lots, mais il est possible de comparer des effets. Disons que j'ai eu deux traitements, plusieurs individus et plusieurs mesures chez des individus. J'utiliserais des individus comme effets aléatoires, car chacun d'entre eux contiendrait x mesures. Ensuite, je voulais savoir si ces deux traitements entraînaient une réponse différente. Puis-je utiliser un effectschevauchement de package et de CI dans ce cas?
Mikko
5
Cette question est plus générale et concerne toute approche standard basée sur un modèle. Vaut peut-être une question distincte. (1) En général, on répond aux questions sur les différences entre les traitements en configurant le modèle de manière à ce que la différence entre les traitements focaux soit un contraste (c’est-à-dire un paramètre estimé) dans le modèle, puis en calculant une valeur p. ou vérifiez si les intervalles de confiance à un niveau alpha particulier incluent zéro. (suite)
Ben Bolker
4
(2) Le chevauchement des IC est au mieux un critère conservateur et approximatif des différences entre paramètres (plusieurs articles ont été publiés sur ce sujet). (3) Il existe un problème séparé / orthogonal avec les comparaisons par paires, selon lequel il faut bien contrôler la multiplicité et la non-indépendance des comparaisons (ceci peut être fait, par exemple, à l'aide des méthodes du multcomppaquet, mais cela nécessite au moins une un peu d'attention)
Ben Bolker
1
Pour quoi? Vous voudrez peut-être poser une nouvelle question.
Ben Bolker
20

Il semble que ce que vous avez fait dans la seconde méthode consiste à calculer des intervalles de confiance pour les coefficients de régression, puis à les transformer pour obtenir des CI pour les prévisions. Cela ignore les covariances entre les coefficients de régression.

Essayez d’ajuster le modèle sans interception, de sorte que les batcheffets correspondent bien aux prévisions et confintrenvoient les intervalles dont vous avez besoin.

Addendum 1

J'ai fait exactement ce que j'ai suggéré ci-dessus:

> fm2 <- lmer(strength ~ batch - 1 + (1 | cask), Pastes)
> confint(fm2)
Computing profile confidence intervals ...
           2.5 %    97.5 %
.sig01  0.000000  1.637468
.sigma  2.086385  3.007380
batchA 60.234772 64.298581
batchB 57.268105 61.331915
batchC 60.018105 64.081915
batchD 57.668105 61.731915
batchE 53.868105 57.931915
batchF 59.001439 63.065248
batchG 57.868105 61.931915
batchH 61.084772 65.148581
batchI 56.651439 60.715248
batchJ 56.551439 60.615248

Ces intervalles semblent correspondre aux résultats de effects.

Addendum 2

Une autre alternative est le paquetage lsmeans . Il obtient des degrés de liberté et une matrice de covariance ajustée à partir du paquet pbkrtest .

> library("lsmeans")
> lsmeans(fm1, "batch")
Loading required namespace: pbkrtest
 batch   lsmean       SE    df lower.CL upper.CL
 A     62.26667 1.125709 40.45 59.99232 64.54101
 B     59.30000 1.125709 40.45 57.02565 61.57435
 C     62.05000 1.125709 40.45 59.77565 64.32435
 D     59.70000 1.125709 40.45 57.42565 61.97435
 E     55.90000 1.125709 40.45 53.62565 58.17435
 F     61.03333 1.125709 40.45 58.75899 63.30768
 G     59.90000 1.125709 40.45 57.62565 62.17435
 H     63.11667 1.125709 40.45 60.84232 65.39101
 I     58.68333 1.125709 40.45 56.40899 60.95768
 J     58.58333 1.125709 40.45 56.30899 60.85768

Confidence level used: 0.95 

Celles-ci sont encore plus en ligne avec les effectrésultats: les erreurs-types sont identiques, mais effectutilise df différent. Les confintrésultats de l'addendum 1 sont encore plus étroits que ceux asymptotiques basés sur l'utilisation de±1,96×se. Alors maintenant, je pense que ceux-ci ne sont pas très fiables.

Les résultats de effectet lsmeanssont similaires, mais avec une situation multifactorielle non équilibrée, les lsmeansmoyennes par défaut sur les facteurs inutilisés ayant une pondération égale, alors que la effectpondération est fonction des fréquences observées (disponible en option dans lsmeans).

RVL
la source
Merci pour cette solution. Les intervalles sont maintenant plus similaires, mais pas exactement les mêmes. Votre réponse ne répond toujours pas à la question de savoir si les éléments de configuration du effectspackage peuvent être approuvés pour les lmerobjets. J'envisage d'utiliser les résultats dans une publication et je veux m'assurer que les CI sont calculés à l'aide d'une méthode approuvée pour les LMM.
Mikko
Pourriez-vous s'il vous plaît dire: dans votre Addendum 1 les deux premiers paramètres .sig01et .sigmaproduit par confint, est-ce que ces intervalles de confiance sont pour la variance ? ou intervalle de confiance de l' écart type ?
ABC
Ce sont des CI pour tous les paramètres étiquetés de cette façon dans le modèle. Vous devriez regarder la documentation lmerpour une réponse définitive. Cependant, les gens utilisent généralement des notations comme sigmase référer à des écarts-types et / sigma.squareou sigma^2à des variances.
vendredi
Est-il préférable d’utiliser lmertest, lsmeans ou mertools?
Skan