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

36

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

Bildbeschreibung hier eingeben

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

Bildbeschreibung hier eingeben

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.

Mikko
quelle
1
Wenn Sie lange Codezeilen haben, würde ich es sehr begrüßen, wenn Sie sie in mehrere Zeilen aufteilen, damit wir nicht scrollen müssen, um alles zu sehen.
rvl
1
@rvl Der Code sollte jetzt besser lesbar sein.
Mikko

Antworten:

52

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:

Bildbeschreibung hier eingeben

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)
Ben Bolker
quelle
2
Ich akzeptiere diese Antwort als zutreffend und gebe einen schönen Vergleich zwischen verschiedenen Methoden. Weitere Informationen finden Sie in der hervorragenden Antwort von rlv.
Mikko
Vielen Dank für die ausgezeichnete und sehr hilfreiche Antwort. Verstehe ich richtig, dass man CIs nicht zum Vergleichen von Gruppen / Chargen verwenden kann, aber es ist möglich, Effekte zu vergleichen. Angenommen, ich hatte zwei Behandlungen, mehrere Einzelpersonen und mehrere Messungen innerhalb von Einzelpersonen. Ich würde Individuen als zufälligen Effekt verwenden, da jede von ihnen x Messungen enthält. Dann wollte ich wissen, ob diese beiden Behandlungen zu einer unterschiedlichen Reaktion führten. Könnte ich effectsin diesem Fall Paket- und CI-Überlappung verwenden?
Mikko
5
Dies ist eine allgemeinere Frage, die für jeden modellbasierten Standardansatz relevant ist. Könnte eine separate Frage wert sein. (1) Im Allgemeinen werden Fragen zu Unterschieden zwischen Behandlungen beantwortet, indem das Modell so eingerichtet wird, dass der Unterschied zwischen den Fokusbehandlungen ein Kontrast (dh ein geschätzter Parameter) im Modell ist, und dann ein p-Wert berechnet wird oder überprüfen Sie, ob die Konfidenzintervalle auf einem bestimmten Alpha-Level Null enthalten. (Fortsetzung)
Ben Bolker
4
(2) Überlappende CIs sind bestenfalls ein konservatives und ungefähres Kriterium für Unterschiede zwischen Parametern (es gibt mehrere veröffentlichte Artikel zu diesem Thema). (3) Bei paarweisen Vergleichen gibt es ein separates / orthogonales Problem, das darin besteht, dass man die Vielzahl und Nichtunabhängigkeit der Vergleiche angemessen kontrollieren muss (dies kann z. B. durch die im multcompPaket enthaltenen Methoden erfolgen, erfordert jedoch mindestens a wenig Pflege)
Ben Bolker
1
Für was? Vielleicht möchten Sie eine neue Frage stellen.
Ben Bolker
20

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.±1.96×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).

rvl
quelle
Vielen Dank für diese Lösung. Die Intervalle sind jetzt ähnlicher, obwohl sie nicht genau gleich sind. Ihre Antwort beantwortet immer noch nicht die Frage, ob CIs aus dem effectsPaket für lmerObjekte vertrauenswürdig sind. Ich erwäge, die Ergebnisse in einer Publikation zu verwenden, und möchte sichergehen, dass CIs mit einer für LMMs zugelassenen Methode berechnet werden.
Mikko
Würden Sie bitte sagen: In Ihrem Addendum 1 sind die ersten beiden Parameter .sig01und deren .sigmaErgebnis confint, das Konfidenzintervall für die Varianz ? oder Konfidenzintervall der Standardabweichung ?
ABC
Sie sind CIs für alle Parameter, die auf diese Weise im Modell gekennzeichnet sind. Sie sollten sich die Dokumentation ansehen, lmerum eine endgültige Antwort zu erhalten. In der Regel wird jedoch gern sigmaauf Standardabweichungen sigma.squareoder sigma^2auf Abweichungen Bezug genommen.
rvl
Ist es besser, lmertest, lsmeans oder mertools zu verwenden?
skan