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

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()

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()

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.

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%:

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)

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$\pm1.96\times\mbox{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).