R: Berechnung des Mittelwerts und des Standardfehlers des Mittelwerts für Faktoren mit lm () vs. direkte Berechnung - bearbeitet

8

Beim Umgang mit Daten mit Faktoren kann R verwendet werden, um die Mittelwerte für jede Gruppe mit der Funktion lm () zu berechnen. Dies gibt auch die Standardfehler für die geschätzten Mittelwerte an. Dieser Standardfehler unterscheidet sich jedoch von dem, was ich aus einer manuellen Berechnung erhalte.

Hier ist ein Beispiel (von hier aus Vorhersage des Unterschieds zwischen zwei Gruppen in R )

Berechnen Sie zuerst den Mittelwert mit lm ():

    mtcars$cyl <- factor(mtcars$cyl)
    mylm <- lm(mpg ~ cyl, data = mtcars)
    summary(mylm)$coef

                Estimate Std. Error   t value     Pr(>|t|)
  (Intercept)  26.663636  0.9718008 27.437347 2.688358e-22
  cyl6         -6.920779  1.5583482 -4.441099 1.194696e-04
  cyl8        -11.563636  1.2986235 -8.904534 8.568209e-10

Der Achsenabschnitt ist der Mittelwert für die erste Gruppe, die 4-Zylinder-Autos. Um die Mittel durch direkte Berechnung zu erhalten, benutze ich Folgendes:

  with(mtcars, tapply(mpg, cyl, mean))

         4        6        8 
    26.66364 19.74286 15.10000 

Um die Standardfehler für die Mittelwerte zu erhalten, berechne ich die Standardvariation der Stichprobe und dividiere durch die Anzahl der Beobachtungen in jeder Gruppe:

 with(mtcars, tapply(mpg, cyl, sd)/sqrt(summary(mtcars$cyl)) )

         4         6         8 
   1.3597642 0.5493967 0.6842016 

Die direkte Berechnung ergibt den gleichen Mittelwert, aber der Standardfehler ist für die beiden Ansätze unterschiedlich. Ich hatte erwartet, den gleichen Standardfehler zu erhalten. Was geht hier vor sich? Es hängt damit zusammen, dass lm () den Mittelwert für jede Gruppe und einen Fehlerterm anpasst.

Bearbeitet: Nach Svens Antwort (unten) kann ich meine Frage präziser und klarer formulieren.

Für kategoriale Daten können wir die Mittelwerte einer Variablen für verschiedene Gruppen berechnen, indem wir lm () ohne Achsenabschnitt verwenden.

  mtcars$cyl <- factor(mtcars$cyl)
  mylm <- lm(mpg ~ cyl, data = mtcars)
  summary(mylm)$coef

      Estimate Std. Error
  cyl4 26.66364  0.9718008
  cyl6 19.74286  1.2182168
  cyl8 15.10000  0.8614094

Wir können dies mit einer direkten Berechnung der Mittelwerte und ihrer Standardfehler vergleichen:

  with(mtcars, tapply(mpg, cyl, mean))

         4        6        8 
    26.66364 19.74286 15.10000 

  with(mtcars, tapply(mpg, cyl, sd)/sqrt(summary(mtcars$cyl)) )

         4         6         8 
   1.3597642 0.5493967 0.6842016 

Die Mittelwerte sind genau gleich, aber die Standardfehler unterscheiden sich für diese beiden Methoden (wie Sven ebenfalls bemerkt). Meine Frage ist, warum sie unterschiedlich und nicht gleich sind.

(Sollte ich beim Bearbeiten meiner Frage den Originaltext löschen oder meine Ausgabe wie bisher hinzufügen)

SRJ
quelle

Antworten:

7

Der Unterschied bei Standardfehlern besteht darin, dass Sie in der Regression eine kombinierte Schätzung der Varianz berechnen , während Sie in der anderen Berechnung separate Schätzungen der Varianz berechnen.

Glen_b -Reinstate Monica
quelle
2
Vielen Dank für die Klarstellung. Ich habe gerade eine sehr gute Antwort auf eine ähnliche Frage gefunden, mit einem gut ausgearbeiteten Beispiel: stats.stackexchange.com/questions/29479/…
SRJ
Ja, das sieht relevant aus. Gut erkannt.
Glen_b -Reinstate Monica
5

Die lmFunktion schätzt nicht Mittelwerte und Standardfehler der Faktorstufen, sondern der mit den Faktorstufen verbundenen Kontrate.

Wenn kein Kontrast manuell angegeben wird, werden in R Behandlungskontraste verwendet. Dies ist die Standardeinstellung für kategoriale Daten.

Der Faktor mtcars$cylhat drei Ebenen (4,6 und 8). Standardmäßig wird die erste Ebene 4 als Referenzkategorie verwendet. Der Achsenabschnitt des linearen Modells entspricht dem Mittelwert der abhängigen Variablen in der Referenzkategorie. Die anderen Effekte ergeben sich jedoch aus einem Vergleich einer Faktorstufe mit der Referenzkategorie. Daher cyl6beziehen sich die Schätzung und der Standardfehler für auf die Differenz zwischen cyl == 6und cyl == 4. Der Effekt cyl8hängt mit dem Unterschied zwischen cyl == 8und zusammen cyl == 4.

Wenn die lmFunktion die Mittelwerte der Faktorstufen berechnen soll, müssen Sie den Intercept-Term ( 0 + ...) ausschließen:

summary(lm(mpg ~ 0 + as.factor(cyl), mtcars))

Call:
lm(formula = mpg ~ 0 + as.factor(cyl), data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.2636 -1.8357  0.0286  1.3893  7.2364 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
as.factor(cyl)4  26.6636     0.9718   27.44  < 2e-16 ***
as.factor(cyl)6  19.7429     1.2182   16.21 4.49e-16 ***
as.factor(cyl)8  15.1000     0.8614   17.53  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1 

Residual standard error: 3.223 on 29 degrees of freedom
Multiple R-squared: 0.9785, Adjusted R-squared: 0.9763 
F-statistic: 440.9 on 3 and 29 DF,  p-value: < 2.2e-16 

Wie Sie sehen können, sind diese Schätzungen mit den Mitteln der Faktorstufen identisch. Beachten Sie jedoch, dass die Standardfehler der Schätzungen nicht mit den Standardfehlern der Daten identisch sind.

Übrigens: Daten können einfach mit der aggregateFunktion aggregiert werden :

aggregate(mpg ~ cyl, mtcars, function(x) c(M = mean(x), SE = sd(x)/sqrt(length(x))))

  cyl      mpg.M     mpg.SE
1   4 26.6636364  1.3597642
2   6 19.7428571  0.5493967
3   8 15.1000000  0.6842016
Sven Hohenstein
quelle
Danke für die Antwort. Ich weiß bereits, dass Koeffizienten nicht das Mittel sind, da ich geschrieben habe, dass der Achsenabschnitt der Mittelwert der ersten Ebene ist, die anderen Koeffizienten sind die Differenz des Mittelwerts der anderen Ebenen zu dieser Ebene. Sie bemerken auch, dass mit Ihrer Bemerkung "Standardfehler der Schätzungen nicht mit den Standardfehlern der Daten identisch sind". Bedeutet das, dass lm () die Mittelwerte schätzt und Standardfehler für diese Schätzungen berechnet
SRJ
Hoppla, ich wollte diesen Kommentar aus Gründen der Übersichtlichkeit bearbeiten, wusste aber nicht, dass ich ihn nur 5 Minuten lang bearbeiten kann. Kann ich einen Kommentar löschen? Ich wusste nicht, dass ich durch Weglassen des Abschnitts direkt mittlere Schätzungen erhalten konnte, danke für diesen Tipp. Wenn ich Sie richtig verstehe, stimmen die Standardfehler der geschätzten Mittelwerte nicht mit den Standardfehlern überein, die direkt aus den Daten berechnet wurden. Handelt es sich jeweils um einen anderen Satz von Gleichungen? Und was sind diese Gleichungen? Ich hätte gerne mehr Details, um den Unterschied besser zu verstehen
SRJ
1

Zusätzlich zu dem, was Sven Hohenstein sagte, sind die mtcarsDaten nicht ausgewogen . Normalerweise verwendet man aovlm mit kategorialen Daten (für die es sich nur um einen Wrapper handelt lm), in denen speziell angegeben ist ?aov:

aov wurde für ausgewogene Designs entwickelt, und die Ergebnisse können ohne Ausgewogenheit schwer zu interpretieren sein: Beachten Sie, dass fehlende Werte in den Antworten wahrscheinlich das Ausgewogenheit verlieren.

Ich denke, Sie können dies auch an den seltsamen Korrelationen der Modellmatrix erkennen:

mf <- model.matrix(mpg ~ cyl, data = mtcars)
cor(mf)
            (Intercept)       cyl6       cyl8
(Intercept)           1         NA         NA
cyl6                 NA  1.0000000 -0.4666667
cyl8                 NA -0.4666667  1.0000000
Warning message:
In cor(mf) : the standard deviation is zero

Daher sind die von aov(oder lm) erhaltenen Standardfehler wahrscheinlich falsch (Sie können dies überprüfen, wenn Sie sie mit Standardfehlern vergleichen lmeoder vergleichen lmer.

Henrik
quelle
Wie würden Sie mich hier bewerben?
SRJ
Die Korrelationen der Modellmatrixwerte sind nicht seltsam. Da die Konstante (Achsenabschnitt) von Natur aus gleich eins ist, gibt es keine Variation zwischen ihren Werten. Aus diesem Grund können Sie keinen Korrelationskoeffizienten zwischen einer Variablen und der Konstanten berechnen.
Sven Hohenstein
-1
Y = matrix(0,5,6)
Y[1,] = c(1250, 980, 1800, 2040, 1000, 1180)
Y[2,] = c(1700, 3080,1700,2820,5760,3480)
Y[3,] = c(2050,3560,2800,1600,4200,2650)
Y[4,] = c(4690,4370,4800,9070,3770,5250)
Y[5,] = c(7150,3480,5010,4810,8740,7260)

n = ncol(Y)
R = rowMeans(Y)
M = mean(R)

s = mean(apply(Y,1,var))

v = var(R)  -s/n


#z = n/(n+(E(s2)/var(m)))
Q = 6/(6+(s/v))
t = Q*R[1] + (1-Z)*M
user257426
quelle
Dies ist unlesbar und es fehlt jeglicher Kommentar darüber, was es bedeutet oder tut. Können Sie es bearbeiten, um mehr Klarheit zu schaffen?
Mdewey
Es ist nicht möglich, die Antwort zu verstehen. Es bedarf einiger Kommentare, um zu erklären, was Sie tun.
Michael R. Chernick