Wie vertrauenswürdig sind die Konfidenzintervalle für ältere Objekte durch das Effektpaket?

Effectspackage bietet eine sehr schnelle und bequeme Möglichkeit , lineare Mischeffekt-Modellergebnisse zu zeichnen, die mit lme4package erhalten wurden . Die effectFunktion berechnet Konfidenzintervalle (CIs) sehr schnell, aber wie vertrauenswürdig sind diese Konfidenzintervalle?

Beispielsweise:

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

Gemäß den mit effectspackage berechneten CIs überschneidet sich Charge "E" nicht mit Charge "A".

Wenn ich dasselbe mit der confint.merModFunktion und der Standardmethode versuche :

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

Ich sehe, dass sich alle CIs überschneiden. Ich erhalte auch Warnungen, die darauf hinweisen, dass die Funktion keine vertrauenswürdigen CIs berechnen konnte. Dieses Beispiel und mein tatsächlicher Datensatz lassen vermuten, dass das effectsPaket Verknüpfungen in der CI-Berechnung enthält, die möglicherweise nicht vollständig von den Statistikern genehmigt wurden. Wie vertrauenswürdig sind die von der effectFunktion zurückgegebenen CIs aus dem effectsPaket für lmerObjekte?

Was habe ich versucht: Beim Blick in den Quellcode ist mir aufgefallen, dass die effectFunktion von der Funktion abhängt Effect.merMod, die wiederum zur Effect.merFunktion führt, die so aussieht:

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.glmDie Funktion scheint die Varianz-Covariate-Matrix aus dem lmerObjekt zu berechnen :

effects:::mer.to.glm

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

Dies wiederum wird wahrscheinlich in der Effect.defaultFunktion zur Berechnung von CIs verwendet (ich habe diesen Teil möglicherweise falsch verstanden):

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
...

Ich weiß nicht genug über LMMs, um beurteilen zu können, ob dies ein richtiger Ansatz ist, aber angesichts der Diskussion um die Konfidenzintervallberechnung für LMMs erscheint dieser Ansatz verdächtig einfach.

Alle Ergebnisse sind im Wesentlichen gleich ( für dieses spezielle Beispiel ). Einige theoretische Unterschiede sind:

• Wie @rvl hervorhebt, ist Ihre Rekonstruktion von CIs ohne Berücksichtigung der Kovarianz zwischen den Parametern einfach falsch. (Entschuldigung)
• Konfidenzintervalle für Parameter basierend auf Wald Konfidenzintervalle werden kann (eine quadratische log-likelihood Oberfläche vorausgesetzt): lsmeans, effects, confint(.,method="Wald"); abgesehen lsmeansdavon ignorieren diese Methoden Effekte endlicher Größe ("Freiheitsgrade"), aber in diesem Fall macht es kaum einen Unterschied ( df=40ist praktisch nicht von unendlich zu unterscheiden df)
• ... oder in Profilvertrauensbereichen (die Standardmethode; ignoriert Effekte mit endlicher Größe, berücksichtigt jedoch nicht quadratische Flächen)
• ... oder beim parametrischen Bootstrapping (der Goldstandard - setzt voraus, dass das Modell korrekt ist [Antworten sind normal, zufällige Effekte sind normal verteilt, Daten sind bedingt unabhängig usw.], aber ansonsten werden nur wenige Annahmen getroffen)

Ich denke, all diese Ansätze sind vernünftig (einige sind näherungsweise als andere), aber in diesem Fall macht es kaum einen Unterschied, welchen Sie verwenden. Wenn Sie Bedenken haben, probieren Sie verschiedene Kontrastmethoden für Ihre Daten oder für simulierte Daten aus, die Ihren eigenen ähneln, und sehen Sie, was passiert ...

(PS: Ich würde nicht zu viel Gewicht auf die Tatsache legen, dass sich die Konfidenzintervalle von Aund Enicht überlappen. Sie müssten ein korrektes paarweises Vergleichsverfahren durchführen, um verlässliche Rückschlüsse auf die Unterschiede zwischen diesem bestimmten Schätzungspaar zu ziehen. ..)

95% CIs:

Vergleichscode:

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)

Es sieht so aus, als hätten Sie bei der zweiten Methode Konfidenzintervalle für die Regressionskoeffizienten berechnet und diese dann transformiert, um CIs für die Vorhersagen zu erhalten. Dies ignoriert die Kovarianzen zwischen den Regressionskoeffizienten.

Versuchen Sie, das Modell ohne Schnittpunkt anzupassen, damit die batchAuswirkungen tatsächlich die Vorhersagen sind und confintdie Intervalle zurückgeben, die Sie benötigen.

Anhang 1

Ich habe genau das gemacht, was ich oben vorgeschlagen habe:

> 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

Diese Intervalle scheinen mit den Ergebnissen von zu übereinstimmen effects.

Anhang 2

Eine weitere Alternative ist das lsmeans- Paket. Es erhält Freiheitsgrade und eine angepasste Kovarianzmatrix aus dem Paket 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 

Diese stimmen noch mehr mit den effectErgebnissen überein: Die Standardfehler sind identisch, effectverwenden jedoch unterschiedliche df. Die confintErgebnisse in Anhang 1 sind noch enger als die asymptotischen, basierend auf der Verwendung von . Also jetzt denke ich, dass diese nicht sehr vertrauenswürdig sind.$\pm1.96\times\mbox{se}$

Ergebnisse aus effectund lsmeanssind ähnlich, jedoch mit einer nicht ausgeglichenen Multi-Faktor-Situation. lsmeansStandardmäßig werden nicht verwendete Faktoren mit gleicher Gewichtung gemittelt, wohingegen effectdie beobachteten Häufigkeiten gewichtet werden (optional in verfügbar lsmeans).