ANCOVA in R schlägt unterschiedliche Abschnitte vor, aber die 95% CIs überlappen sich… wie ist das möglich?

Wir haben einen Datensatz mit zwei Kovariaten und einer kategorialen Gruppierungsvariablen und möchten wissen, ob es signifikante Unterschiede zwischen der Steigung oder dem Achsenabschnitt zwischen den Kovariaten gibt, die den verschiedenen Gruppierungsvariablen zugeordnet sind. Wir haben anova () und lm () verwendet, um die Anpassungen von drei verschiedenen Modellen zu vergleichen: 1) mit einer einzelnen Steigung und einem Achsenabschnitt, 2) mit unterschiedlichen Abschnitten für jede Gruppe und 3) mit einer Steigung und einem Achsenabschnitt für jede Gruppe . Nach dem allgemeinen linearen Test von anova () ist das zweite Modell das am besten geeignete der drei. Es gibt eine signifikante Verbesserung des Modells, indem für jede Gruppe ein separater Abschnitt eingefügt wird. Wenn wir uns jedoch die 95% -Konfidenzintervalle für diese Abschnitte ansehen, überlappen sich alle, was darauf hindeutet, dass es keine signifikanten Unterschiede zwischen den Abschnitten gibt. Wie können diese beiden Ergebnisse in Einklang gebracht werden? Wir dachten, eine andere Art der Interpretation der Ergebnisse der Modellauswahlmethode wäre, dass es mindestens einen signifikanten Unterschied zwischen den Abschnitten geben muss ... aber vielleicht ist dies nicht korrekt?

Unten finden Sie den R-Code zum Replizieren dieser Analyse. Wir haben die Funktion dput () verwendet, damit Sie mit genau den Daten arbeiten können, mit denen wir uns auseinandersetzen.

# Begin R Script
# > dput(data)
structure(list(Head = c(1.92, 1.93, 1.79, 1.94, 1.91, 1.88, 1.91, 
1.9, 1.97, 1.97, 1.95, 1.93, 1.95, 2, 1.87, 1.88, 1.97, 1.88, 
1.89, 1.86, 1.86, 1.97, 2.02, 2.04, 1.9, 1.83, 1.95, 1.87, 1.93, 
1.94, 1.91, 1.96, 1.89, 1.87, 1.95, 1.86, 2.03, 1.88, 1.98, 1.97, 
1.86, 2.04, 1.86, 1.92, 1.98, 1.86, 1.83, 1.93, 1.9, 1.97, 1.92, 
2.04, 1.92, 1.9, 1.93, 1.96, 1.91, 2.01, 1.97, 1.96, 1.76, 1.84, 
1.92, 1.96, 1.87, 2.1, 2.17, 2.1, 2.11, 2.17, 2.12, 2.06, 2.06, 
2.1, 2.05, 2.07, 2.2, 2.14, 2.02, 2.08, 2.16, 2.11, 2.29, 2.08, 
2.04, 2.12, 2.02, 2.22, 2.22, 2.2, 2.26, 2.15, 2, 2.24, 2.18, 
2.07, 2.06, 2.18, 2.14, 2.13, 2.2, 2.1, 2.13, 2.15, 2.25, 2.14, 
2.07, 1.98, 2.16, 2.11, 2.21, 2.18, 2.13, 2.06, 2.21, 2.08, 1.88, 
1.81, 1.87, 1.88, 1.87, 1.79, 1.99, 1.87, 1.95, 1.91, 1.99, 1.85, 
2.03, 1.88, 1.88, 1.87, 1.85, 1.94, 1.98, 2.01, 1.82, 1.85, 1.75, 
1.95, 1.92, 1.91, 1.98, 1.92, 1.96, 1.9, 1.86, 1.97, 2.06, 1.86, 
1.91, 2.01, 1.73, 1.97, 1.94, 1.81, 1.86, 1.99, 1.96, 1.94, 1.85, 
1.91, 1.96, 1.9, 1.98, 1.89, 1.88, 1.95, 1.9, 1.94, NA, 1.84, 
1.83, 1.84, 1.96, 1.74, 1.91, 1.84, 1.88, 1.83, 1.93, 1.78, 1.88, 
1.93, 2.15, 2.16, 2.23, 2.09, 2.36, 2.31, 2.25, 2.29, 2.3, 2.04, 
2.22, 2.19, 2.25, 2.31, 2.3, 2.28, 2.25, 2.15, 2.29, 2.24, 2.34, 
2.2, 2.24, 2.17, 2.26, 2.18, 2.17, 2.34, 2.23, 2.36, 2.31, 2.13, 
2.2, 2.27, 2.27, 2.2, 2.34, 2.12, 2.26, 2.18, 2.31, 2.24, 2.26, 
2.15, 2.29, 2.14, 2.25, 2.31, 2.13, 2.09, 2.24, 2.26, 2.26, 2.21, 
2.25, 2.29, 2.15, 2.2, 2.18, 2.16, 2.14, 2.26, 2.22, 2.12, 2.12, 
2.16, 2.27, 2.17, 2.27, 2.17, 2.3, 2.25, 2.17, 2.27, 2.06, 2.13, 
2.11, 2.11, 1.97, 2.09, 2.06, 2.11, 2.09, 2.08, 2.17, 2.12, 2.13, 
1.99, 2.08, 2.01, 1.97, 1.97, 2.09, 1.94, 2.06, 2.09, 2.04, 2, 
2.14, 2.07, 1.98, 2, 2.19, 2.12, 2.06, 2, 2.02, 2.16, 2.1, 1.97, 
1.97, 2.1, 2.02, 1.99, 2.13, 2.05, 2.05, 2.16, 2.02, 2.02, 2.08, 
1.98, 2.04, 2.02, 2.07, 2.02, 2.02, 2.02), Site = structure(c(2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L), .Label = c("ANZ", "BC", "DV", "MC", 
"RB", "WW"), class = "factor"), Leg = c(2.38, 2.45, 2.22, 2.23, 
2.26, 2.32, 2.28, 2.17, 2.39, 2.27, 2.42, 2.33, 2.31, 2.32, 2.25, 
2.27, 2.38, 2.28, 2.33, 2.24, 2.21, 2.22, 2.42, 2.23, 2.36, 2.2, 
2.28, 2.23, 2.33, 2.35, 2.36, 2.26, 2.26, 2.3, 2.23, 2.31, 2.27, 
2.23, 2.37, 2.27, 2.26, 2.3, 2.33, 2.34, 2.27, 2.4, 2.22, 2.25, 
2.28, 2.33, 2.26, 2.32, 2.29, 2.31, 2.37, 2.24, 2.26, 2.36, 2.32, 
2.32, 2.15, 2.2, 2.29, 2.37, 2.26, 2.24, 2.23, 2.24, 2.26, 2.18, 
2.11, 2.23, 2.31, 2.25, 2.15, 2.3, 2.33, 2.35, 2.21, 2.36, 2.27, 
2.24, 2.35, 2.24, 2.33, 2.32, 2.24, 2.35, 2.36, 2.39, 2.28, 2.36, 
2.19, 2.27, 2.39, 2.23, 2.29, 2.32, 2.3, 2.32, NA, 2.25, 2.24, 
2.21, 2.37, 2.21, 2.21, 2.27, 2.27, 2.26, 2.19, 2.2, 2.25, 2.25, 
2.25, NA, 2.24, 2.17, 2.2, 2.2, 2.18, 2.14, 2.17, 2.27, 2.28, 
2.27, 2.29, 2.23, 2.25, 2.33, 2.22, 2.29, 2.19, 2.15, 2.24, 2.24, 
2.26, 2.25, 2.09, 2.27, 2.18, 2.2, 2.25, 2.24, 2.18, 2.3, 2.26, 
2.18, 2.27, 2.12, 2.18, 2.33, 2.13, 2.28, 2.23, 2.16, 2.2, 2.3, 
2.31, 2.18, 2.33, 2.29, 2.26, 2.21, 2.22, 2.27, 2.32, 2.24, 2.25, 
2.17, 2.2, 2.26, 2.27, 2.24, 2.25, 2.09, 2.25, 2.21, 2.24, 2.21, 
2.22, 2.13, 2.24, 2.21, 2.3, 2.34, 2.35, 2.32, 2.46, 2.43, 2.42, 
2.41, 2.32, 2.25, 2.33, 2.19, 2.45, 2.32, 2.4, 2.38, 2.35, 2.39, 
2.29, 2.35, 2.43, 2.29, 2.33, 2.31, 2.28, 2.38, 2.32, 2.43, 2.27, 
2.4, 2.37, 2.27, 2.41, 2.32, 2.38, 2.23, 2.33, 2.21, 2.34, 2.19, 
2.34, 2.35, 2.35, 2.31, 2.33, 2.41, 2.53, 2.39, 2.17, 2.16, 2.38, 
2.34, 2.33, 2.33, 2.29, 2.43, 2.28, 2.34, 2.38, 2.3, 2.29, 2.43, 
2.36, 2.24, 2.35, 2.38, 2.4, 2.36, 2.42, 2.28, 2.45, 2.33, 2.32, 
2.33, 2.31, 2.44, 2.37, 2.4, 2.35, 2.33, 2.31, 2.36, 2.43, 2.38, 
2.4, 2.38, 2.46, 2.33, 2.38, 2.23, 2.24, 2.39, 2.36, 2.19, 2.32, 
2.37, 2.39, 2.34, 2.39, 2.23, 2.25, 2.29, 2.39, 2.35, NA, 2.28, 
2.35, 2.38, 2.34, 2.17, 2.29, NA, 2.26, NA, NA, NA, 2.24, 2.33, 
2.23, 2.28, 2.29, 2.23, 2.2, 2.27, 2.31, 2.31, 2.26, 2.28)), .Names = c("Head", 
"Site", "Leg"), class = "data.frame", row.names = c(NA, -312L
)) 

# plot graph
library(ggplot2)

qplot(Head, Leg, 
    color=Site, 
    data=data) + 
        stat_smooth(method="lm", alpha=0.2) + 
        theme_bw()

# create linear models
lm.1 <- lm(Leg ~ Head, data)
lm.2 <- lm(Leg ~ Head + Site, data)
lm.3 <- lm(Leg ~ Head*Site, data)

# evaluate linear models
anova(lm.1, lm.2, lm.3)
anova(lm.1, lm.2)

# > anova(lm.1, lm.2)
# Analysis of Variance Table
# Model 1: Leg.3.1 ~ Head.W1
# Model 2: Leg.3.1 ~ Head.W1 + Site
  # Res.Df     RSS Df Sum of Sq     F    Pr(>F)    
# 1    302 1.25589                                 
# 2    297 0.91332  5   0.34257 22.28 < 2.2e-16 ***


# examining the multiple-intercepts model (lm.2)
summary(lm.2)
coef(lm.2)
confint(lm.2)

# extracting the intercepts
intercepts <- coef(lm.2)[c(1, 3:7)]
intercepts.1 <- intercepts[1]
intercepts <- intercepts.1 + intercepts
intercepts[1] <- intercepts.1
intercepts

# extracting the confidence intervals
ci <- confint(lm.2)[c(1, 3:7),]
ci[2:6,] <- ci[2:6,] + confint(lm.2)[1,]
ci[,1]

# putting everything together in a dataframe
labels <- c("ANZ", "BC", "DV", "MC", "RB", "WW")
ci.dataframe <- data.frame(Site=labels, Intercept=intercepts, CI.low = ci[,1], CI.high = ci[,2])
ci.dataframe

# plotting intercepts and 95% CI
qplot(Site, Intercept, geom=c("point", "errorbar"), ymin=CI.low, ymax=CI.high, data=ci.dataframe, ylab="Intercept & 95% CI")

Zusammenfassend lässt sich sagen, dass sich die 95% -KI für die Abschnitte alle überlappen. Die Modellauswahlmethode legt jedoch nahe, dass das beste Modell für verschiedene Abschnitte geeignet ist. Daher neige ich dazu zu glauben, dass entweder unsere Modellauswahlmethode fehlerhaft ist oder die 95% CIs für die Intercept-Schätzungen falsch berechnet wurden. Irgendwelche Gedanken wären sehr dankbar!

Denken Sie daran, dass der Unterschied zwischen signifikant und nicht signifikant nicht (immer) statistisch signifikant ist

Genauer gesagt heißt Modell 1 gepoolte Regression und Modell 2 nicht gepoolte Regression. Wie Sie bereits bemerkt haben, gehen Sie bei der gepoolten Regression davon aus, dass die Gruppen nicht relevant sind. Dies bedeutet, dass die Varianz zwischen den Gruppen auf Null gesetzt wird.

In der nicht gepoolten Regression mit einem Achsenabschnitt pro Gruppe setzen Sie die Varianz auf unendlich.

Im Allgemeinen würde ich eine Zwischenlösung bevorzugen, bei der es sich um ein hierarchisches Modell oder eine teilweise gepoolte Regression (oder einen Schrumpfungsschätzer) handelt. Sie können dieses Modell mit dem lmer4-Paket in R einbauen.

Schauen Sie sich abschließend dieses Papier von Gelman an , in dem er argumentiert, warum hierarchische Modelle bei den Problemen mit mehreren Vergleichen helfen (in Ihrem Fall sind die Koeffizienten pro Gruppe unterschiedlich? Wie korrigieren wir einen p-Wert für mehrere Vergleiche).

Zum Beispiel in Ihrem Fall

library(lme4)
summary(lmer( leg ~ head + (1 | site)) # varying intercept model

Wenn Sie eine unterschiedliche Steigung mit variierendem Achsenabschnitt (das dritte Modell) anpassen möchten, führen Sie einfach aus

summary(lmer( leg ~ head + (1 | site) + (0+head|site) )) # varying intercept, varying-slope model

Dann können Sie einen Blick auf die Gruppenvarianz werfen und feststellen, ob sie sich von Null (die gepoolte Regression ist nicht das bessere Modell) und weit von der Unendlichkeit (nicht gepoolte Regression) unterscheidet.

Update: Nach den Kommentaren (siehe unten) habe ich beschlossen, meine Antwort zu erweitern.

Der Zweck eines hierarchischen Modells, insbesondere in solchen Fällen, besteht darin, die Variation nach Gruppen (in diesem Fall Sites) zu modellieren. Anstatt eine ANOVA auszuführen, um zu testen, ob sich ein Modell von einem anderen unterscheidet, würde ich mir die Vorhersagen meines Modells ansehen und prüfen, ob die Vorhersagen nach Gruppen in den hierarchischen Modellen besser sind als in der gepoolten Regression (klassische Regression). .

Jetzt habe ich meine Vorschläge oben ausgeführt und das bestätigt

ranef(lmer( leg ~ head + (1 | site) + (0+head|site) )

Würde Null als Schätzung der variierenden Steigung zurückgeben (variierender Effekt des Kopfes nach Standort). dann rannte ich

ranef(lmer( leg ~ head + (head| site))

Und ich habe eine Schätzung ungleich Null für die unterschiedliche Wirkung des Kopfes erhalten. Ich weiß noch nicht, warum das passiert ist, da ich es zum ersten Mal gefunden habe. Dieses Problem tut mir wirklich leid, aber zu meiner Verteidigung habe ich mich nur an die Spezifikation gehalten, die in der Hilfe der lmer-Funktion beschrieben ist. (Siehe das Beispiel mit der Datenschlafstudie). Ich werde versuchen zu verstehen, was passiert, und ich werde hier berichten, wenn (wenn) ich verstehe, was passiert.

Vor jeder Moderatorintervention können Sie sich ansehen

library(car)

crPlots(lm.2,terms=~Site)

Dies sind Komponenten + Restdiagramme (Teilreste)

Ich denke unter anderem, dass Sie die Konfidenzintervalle falsch berechnen. Hier sind zwei Möglichkeiten, es zu betrachten:

(1) Unterschiede zwischen den einzelnen Standorten und dem ANZ-Standort (Baseline) [Sie können auch Unterschiede zum Gesamtmittelwert berechnen, indem Sie zu Kontrasten von Summe zu Null wechseln

library(coefplot2)  ## on r-forge
coefplot2(lm.2)

oder (2) alle paarweisen Vergleiche (ich mag diesen Ansatz nicht, aber er ist üblich):

library(multcomp)
ci <- confint(glht(lm.2, linfct = mcp(Site = "Tukey")))
ggplot(fortify(ci),aes(lhs,estimate,ymin=lwr,ymax=upr))+
    geom_pointrange()+theme_bw()+geom_hline(yintercept=0,col="red")

Beachten Sie, dass alle Ihre HeadWerte im Bereich von 1,7 bis 2,4 liegen, während die Abschnitte versuchen, den LegWert auf zu schätzen Head=0. Dies ist eine große Extrapolation, daher gibt es viele Unsicherheiten. Wenn Sie die HeadWerte zentrieren und diese Analyse wiederholen, werden die Konfidenzintervalle viel enger.

Darüber hinaus bedeuten überlappende 95% -Konfidenzintervalle keinen Mangel an statistisch signifikanten Unterschieden. Tatsächlich nähern sich für zwei Gruppen nicht überlappende 84% -Konfidenzintervalle ungefähr signifikanten Unterschieden bei 5% an. Aufgrund mehrerer Tests funktioniert dies natürlich nicht mit mehreren Gruppen.

Zusätzlich zu den anderen Antworten finden Sie hier einige Links der Cornell Statistical Consulting Unit, die für überlappende Konfidenzintervalle relevant sind und als gute, kurze Erinnerung daran dienen, was sie tun und was nicht

http://www.cscu.cornell.edu/news/statnews/stnews73.pdf http://www.cscu.cornell.edu/news/statnews/Stnews73insert.pdf

Hier ist der Hauptpunkt:

Wenn zwei Statistiken nicht überlappende Konfidenzintervall haben, sind sie notwendigerweise signifikant unterschiedlich , aber wenn sie Konfidenzintervall überlappende haben, ist es nicht unbedingt wahr , dass sie nicht signifikant verschieden sind.

Hier ist der relevante Text aus dem ersten Link:

Wir können dies anhand eines einfachen Beispiels veranschaulichen. Angenommen, wir möchten Mittel aus zwei unabhängigen Stichproben vergleichen. Der Mittelwert der ersten Stichprobe beträgt 9 und der Mittelwert der zweiten Stichprobe 17. Nehmen wir an, dass die beiden Gruppenmittelwerte die gleichen Standardfehler von 2,5 aufweisen. Das 95-Prozent-Konfidenzintervall für den ersten Gruppenmittelwert kann wie folgt berechnet werden: ± × 5,296,19 wobei 1,96 der kritische t-Wert ist. Das Konfidenzintervall für den Mittelwert der ersten Gruppe beträgt somit (4.1, 13.9). In ähnlicher Weise beträgt das Konfidenzintervall für den Mittelwert für die zweite Gruppe (12,1, 21,9). Beachten Sie, dass sich die beiden Intervalle überschneiden. Die t-Statistik zum Vergleichen zweier Mittelwerte lautet jedoch:

t = (17-9)/√(2.5² + 2.5²) = 2.26

was widerspiegelt, dass die Nullhypothese, dass die Mittelwerte der beiden Gruppen gleich sind, auf der Ebene von α = 0,05 verworfen werden sollte. Um die obige Schlussfolgerung zu überprüfen, betrachten Sie das 95-Prozent-Konfidenzintervall für die Differenz zwischen den beiden Gruppenmitteln: (17-9) ± 1,96 x √ (2,5² + 2,5²), was (1,09, 14,91) ergibt. Das Intervall enthält keine Null, daher lehnen wir die Nullhypothese ab, dass die Gruppenmittelwerte gleich sind.

Im Allgemeinen ist es beim Vergleich zweier Parameterschätzungen immer richtig, dass sich die Statistiken statistisch signifikant unterscheiden, wenn sich die Konfidenzintervalle nicht überschneiden. Das Gegenteil ist jedoch nicht der Fall. Das heißt, es ist falsch, die statistische Signifikanz der Differenz zwischen zwei Statistiken basierend auf überlappenden Konfidenzintervallen zu bestimmen. Eine Erklärung, warum dies für den Vergleich von Mitteln mit zwei Stichproben gilt, finden Sie unter folgendem Link: http://www.cscu.cornell.edu/news/statnews/Stnews73insert.pdf

Hier sind die Informationen vom anderen Link: