R: Menghitung rata-rata dan kesalahan standar rata-rata untuk faktor-faktor dengan lm () vs. perhitungan langsung -ditandai

8

Ketika berhadapan dengan data dengan faktor R dapat digunakan untuk menghitung rata-rata untuk setiap kelompok dengan fungsi lm (). Ini juga memberikan kesalahan standar untuk estimasi cara. Tetapi kesalahan standar ini berbeda dari apa yang saya dapatkan dari perhitungan dengan tangan.

Ini adalah contoh (diambil dari sini. Memprediksi perbedaan antara dua grup dalam R )

Pertama menghitung mean dengan 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

Mencegat adalah rata-rata untuk kelompok pertama, 4 mobil silinder. Untuk mendapatkan sarana dengan perhitungan langsung, saya menggunakan ini:

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

         4        6        8 
    26.66364 19.74286 15.10000 

Untuk mendapatkan kesalahan standar untuk cara saya menghitung variasi standar sampel dan membaginya dengan jumlah pengamatan di setiap kelompok:

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

         4         6         8 
   1.3597642 0.5493967 0.6842016 

Perhitungan langsung memberikan rata-rata yang sama tetapi kesalahan standar berbeda untuk 2 pendekatan, saya berharap mendapatkan kesalahan standar yang sama. Apa yang terjadi disini? Ini terkait dengan lm () yang cocok dengan rata-rata untuk setiap kelompok dan istilah kesalahan?

Diedit: Setelah jawaban Svens (di bawah) saya dapat merumuskan pertanyaan saya lebih ringkas dan jelas.

Untuk data kategorikal kita dapat menghitung rata-rata variabel untuk kelompok yang berbeda adalah dengan menggunakan lm () tanpa intersep.

  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

Kita dapat membandingkan ini dengan perhitungan langsung dari cara dan kesalahan standar mereka:

  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 

Berarti persis sama tetapi kesalahan standar berbeda untuk 2 metode ini (seperti Sven juga perhatikan). Pertanyaan saya adalah mengapa mereka berbeda dan tidak sama?

(saat mengedit pertanyaan saya, haruskah saya menghapus teks asli atau menambahkan edisi saya seperti yang saya lakukan)

SRJ
sumber

Jawaban:

7

Perbedaan dalam kesalahan standar adalah karena dalam regresi Anda menghitung estimasi gabungan varians, sedangkan dalam perhitungan lain Anda menghitung estimasi varians terpisah.

Glen_b -Reinstate Monica
sumber
2
Terima kasih telah menjelaskan ini. Saya baru saja menemukan jawaban yang sangat bagus untuk pertanyaan serupa di sini, dengan contoh yang bagus: stats.stackexchange.com/questions/29479/…
SRJ
Ya, itu terlihat relevan. Terlihat dengan baik.
Glen_b -Reinstate Monica
5

The lmFungsi tidak memperkirakan sarana dan kesalahan standar dari tingkat faktor tapi contrats terkait dengan tingkat faktor.

Jika tidak ada kontras yang ditentukan secara manual, kontras pengobatan digunakan dalam R. Ini adalah default untuk data kategorikal.

Faktor ini mtcars$cylmemiliki tiga tingkatan (4,6, dan 8). Secara default, level pertama, 4, digunakan sebagai kategori referensi. Intersep model linier berhubungan dengan rata-rata variabel dependen dalam kategori referensi. Tetapi efek lainnya dihasilkan dari perbandingan satu tingkat faktor dengan kategori referensi. Oleh karena itu, estimasi dan kesalahan standar cyl6terkait dengan perbedaan antara cyl == 6dan cyl == 4. Efeknya cyl8terkait dengan perbedaan antara cyl == 8dan cyl == 4.

Jika Anda ingin lmfungsi untuk menghitung rata-rata tingkat faktor, Anda harus mengecualikan istilah intersep ( 0 + ...):

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 

Seperti yang Anda lihat, perkiraan ini identik dengan rata-rata tingkat faktor. Tetapi perhatikan bahwa kesalahan standar estimasi tidak identik dengan kesalahan standar data.

Omong-omong: Data dapat dikumpulkan dengan mudah dengan aggregatefungsi:

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
sumber
Terima kasih atas tanggapannya. Saya sudah tahu bahwa koefisien bukan berarti, karena saya menulis intersep adalah rata-rata dari tingkat pertama, koefisien lainnya adalah perbedaan rata-rata dari tingkat lain ke tingkat itu. Anda juga memperhatikan bahwa dengan komentar Anda, "kesalahan standar estimasi tidak identik dengan kesalahan standar data." Apakah itu berarti bahwa lm () memperkirakan rata-rata dan menghitung kesalahan standar untuk perkiraan tersebut
SRJ
Ooops, saya ingin mengedit komentar itu untuk kejelasan tetapi tidak tahu saya hanya dapat mengedit selama 5 menit, dapatkah saya menghapus komentar? Saya tidak menyadari bahwa saya bisa mendapatkan perkiraan rata-rata secara langsung dengan menghilangkan intersepsi, terima kasih atas tipnya. Jika saya mengerti Anda dengan benar, kesalahan standar dari perkiraan berarti tidak sama dengan kesalahan standar yang dihitung langsung dari data. Apakah persamaan yang digunakan dalam setiap kasus berbeda? Dan apa persamaan ini? Saya ingin memiliki lebih banyak detail untuk memahami perbedaannya lebih baik
SRJ
1

Selain apa yang dikatakan Sven Hohenstein, mtcarsdata tidak seimbang . Biasanya satu menggunakan aovuntuk lm dengan data kategorikal (yang hanya pembungkus untuk lm) yang secara khusus mengatakan pada ?aov:

aov dirancang untuk desain yang seimbang, dan hasilnya bisa sulit untuk diartikan tanpa keseimbangan: berhati-hatilah bahwa nilai yang hilang dalam respons akan kehilangan keseimbangan.

Saya pikir Anda juga bisa melihat ini pada korelasi aneh dari matriks model:

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

Oleh karena itu, kesalahan standar yang diperoleh dari aov(atau lm) kemungkinan akan palsu (Anda dapat memeriksa ini jika Anda membandingkan dengan lmeatau lmerkesalahan standar.

Henrik
sumber
Bagaimana Anda akan melamar saya di sini?
SRJ
Korelasi dari nilai-nilai model matriks tidak aneh. Karena konstanta (intersep) secara inheren sama dengan satu, tidak ada variasi antara nilainya. Karena ini, Anda tidak dapat menghitung koefisien korelasi antara variabel dan konstanta.
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
pengguna257426
sumber
Ini tidak dapat dibaca dan tidak memiliki komentar apa pun tentang apa artinya atau apa artinya. Bisakah Anda mengeditnya untuk memberikan lebih banyak kejelasan?
mdewey
Tidak mungkin untuk memahami jawabannya. Memang diperlukan beberapa komentar untuk menjelaskan apa yang Anda lakukan.
Michael R. Chernick