Saya mencoba untuk menganalisis efek Tahun pada variabel logInd untuk kelompok individu tertentu (saya punya 3 kelompok). Model paling sederhana:
> fix1 = lm(logInd ~ 0 + Group + Year:Group, data = mydata)
> summary(fix1)
Call:
lm(formula = logInd ~ 0 + Group + Year:Group, data = mydata)
Residuals:
Min 1Q Median 3Q Max
-5.5835 -0.3543 -0.0024 0.3944 4.7294
Coefficients:
Estimate Std. Error t value Pr(>|t|)
Group1 4.6395740 0.0466217 99.515 < 2e-16 ***
Group2 4.8094268 0.0534118 90.044 < 2e-16 ***
Group3 4.5607287 0.0561066 81.287 < 2e-16 ***
Group1:Year -0.0084165 0.0027144 -3.101 0.00195 **
Group2:Year 0.0032369 0.0031098 1.041 0.29802
Group3:Year 0.0006081 0.0032666 0.186 0.85235
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.7926 on 2981 degrees of freedom
Multiple R-squared: 0.9717, Adjusted R-squared: 0.9716
F-statistic: 1.705e+04 on 6 and 2981 DF, p-value: < 2.2e-16
Kita dapat melihat bahwa Group1 menurun secara signifikan, Groups2 dan 3 meningkat tetapi tidak begitu signifikan.
Jelas individu harus efek acak, jadi saya memperkenalkan efek intersepsi acak untuk setiap individu:
> mix1a = lmer(logInd ~ 0 + Group + Year:Group + (1|Individual), data = mydata)
> summary(mix1a)
Linear mixed model fit by REML
Formula: logInd ~ 0 + Group + Year:Group + (1 | Individual)
Data: mydata
AIC BIC logLik deviance REMLdev
4727 4775 -2356 4671 4711
Random effects:
Groups Name Variance Std.Dev.
Individual (Intercept) 0.39357 0.62735
Residual 0.24532 0.49530
Number of obs: 2987, groups: Individual, 103
Fixed effects:
Estimate Std. Error t value
Group1 4.6395740 0.1010868 45.90
Group2 4.8094268 0.1158095 41.53
Group3 4.5607287 0.1216522 37.49
Group1:Year -0.0084165 0.0016963 -4.96
Group2:Year 0.0032369 0.0019433 1.67
Group3:Year 0.0006081 0.0020414 0.30
Correlation of Fixed Effects:
Group1 Group2 Group3 Grp1:Y Grp2:Y
Group2 0.000
Group3 0.000 0.000
Group1:Year -0.252 0.000 0.000
Group2:Year 0.000 -0.252 0.000 0.000
Group3:Year 0.000 0.000 -0.252 0.000 0.000
Ini memiliki efek yang diharapkan - SE lereng (koefisien Group1-3: Tahun) sekarang lebih rendah dan residu SE juga lebih rendah.
Individu juga berbeda dalam kemiringan jadi saya juga memperkenalkan efek kemiringan acak:
> mix1c = lmer(logInd ~ 0 + Group + Year:Group + (1 + Year|Individual), data = mydata)
> summary(mix1c)
Linear mixed model fit by REML
Formula: logInd ~ 0 + Group + Year:Group + (1 + Year | Individual)
Data: mydata
AIC BIC logLik deviance REMLdev
2941 3001 -1461 2885 2921
Random effects:
Groups Name Variance Std.Dev. Corr
Individual (Intercept) 0.1054790 0.324775
Year 0.0017447 0.041769 -0.246
Residual 0.1223920 0.349846
Number of obs: 2987, groups: Individual, 103
Fixed effects:
Estimate Std. Error t value
Group1 4.6395740 0.0541746 85.64
Group2 4.8094268 0.0620648 77.49
Group3 4.5607287 0.0651960 69.95
Group1:Year -0.0084165 0.0065557 -1.28
Group2:Year 0.0032369 0.0075105 0.43
Group3:Year 0.0006081 0.0078894 0.08
Correlation of Fixed Effects:
Group1 Group2 Group3 Grp1:Y Grp2:Y
Group2 0.000
Group3 0.000 0.000
Group1:Year -0.285 0.000 0.000
Group2:Year 0.000 -0.285 0.000 0.000
Group3:Year 0.000 0.000 -0.285 0.000 0.000
Tapi sekarang, bertentangan dengan harapan, SE lereng (koefisien Group1-3: Tahun) sekarang jauh lebih tinggi, bahkan lebih tinggi daripada tanpa efek acak sama sekali!
Bagaimana ini mungkin? Saya akan berharap bahwa efek acak akan "memakan" variabilitas yang tidak dapat dijelaskan dan meningkatkan "kepastian" dari perkiraan!
Namun, sisa SE berperilaku seperti yang diharapkan - ini lebih rendah dari pada model intersep acak.
Ini data jika diperlukan.
Edit
Sekarang saya menyadari fakta yang mencengangkan. Jika saya melakukan regresi linier untuk setiap individu secara terpisah dan kemudian menjalankan ANOVA pada lereng yang dihasilkan, saya mendapatkan hasil yang persis sama dengan model lereng acak! Apakah kamu tahu mengapa?
indivSlope = c()
for (indiv in 1:103) {
mod1 = lm(logInd ~ Year, data = mydata[mydata$Individual == indiv,])
indivSlope[indiv] = coef(mod1)['Year']
}
indivGroup = unique(mydata[,c("Individual", "Group")])[,"Group"]
anova1 = lm(indivSlope ~ 0 + indivGroup)
summary(anova1)
Call:
lm(formula = indivSlope ~ 0 + indivGroup)
Residuals:
Min 1Q Median 3Q Max
-0.176288 -0.016502 0.004692 0.020316 0.153086
Coefficients:
Estimate Std. Error t value Pr(>|t|)
indivGroup1 -0.0084165 0.0065555 -1.284 0.202
indivGroup2 0.0032369 0.0075103 0.431 0.667
indivGroup3 0.0006081 0.0078892 0.077 0.939
Residual standard error: 0.04248 on 100 degrees of freedom
Multiple R-squared: 0.01807, Adjusted R-squared: -0.01139
F-statistic: 0.6133 on 3 and 100 DF, p-value: 0.6079
Ini data jika diperlukan.
sumber
Group
Group
:Year
adalah kemiringan dalam GruplogInd ~ Year*Group
, hanya koefisien-koefisien yang ada dalam bentuk yang berbeda, tidak lebih. Tergantung selera Anda dan bentuk koefisien apa yang Anda sukai, tidak lebih. Tidak ada pengecualian "Efek utama tahun" pada model pertama saya saat Anda menulis ...logInd ~ Year*Group
melakukan hal yang persis sama,Year
koefisiennya bukan efek utama, tetapi Group1: Tahun.Jawaban:
Saya pikir masalahnya adalah dengan harapan Anda :) Perhatikan bahwa ketika Anda menambahkan intersep acak untuk setiap individu, kesalahan standar intersep meningkat. Karena setiap individu dapat memiliki intersepnya sendiri, rata-rata grup kurang pasti. Hal yang sama terjadi dengan kemiringan acak: Anda tidak lagi memperkirakan satu kemiringan (dalam-kelompok) yang umum, tetapi rata-rata kemiringan yang bervariasi.
EDIT: Mengapa model yang lebih baik tidak memberikan perkiraan yang lebih tepat?
Mari kita pikirkan sebaliknya: mengapa model awal meremehkan kesalahan standar? Ini mengasumsikan independensi pengamatan yang tidak independen. Model kedua melonggarkan asumsi itu (dengan cara yang memengaruhi penyadapan), dan yang ketiga melemaskannya lebih jauh.
EDIT 2: hubungan dengan banyak model khusus pasien
Pengamatan Anda adalah properti yang dikenal (dan jika Anda hanya memiliki dua tahun, maka model efek acak akan setara dengan uji-t berpasangan). Saya tidak berpikir saya bisa mengelola bukti nyata, tetapi mungkin menulis dua model akan membuat hubungan lebih jelas. Mari kita abaikan variabel pengelompokan, karena hanya akan memperumit notasi. Saya akan menggunakan huruf Yunani untuk efek acak, dan huruf latin untuk efek tetap.
Model efek acak adalah (i - subjek, j - mereplikasi dalam subjek):
Bila Anda cocok dengan model terpisah untuk setiap subjek, maka
[Catatan: yang berikut ini hanya sekedar handwaving:]
Anda dapat melihat banyak kesamaan antara kedua model ini denganai sesuai dengan a+αi dan bi untuk b+βi . Rata-ratabi berkorespondensi dengan b , karena efek acak rata-rata menjadi 0. Korelasi tidak menentu dari intersep dan kemiringan acak mengarah pada fakta bahwa model hanya dapat dipasang secara terpisah. Saya tidak yakin bagaimana singleσ asumsi jala dengan subjek-spesifik σi , tapi saya akan menganggap itu αi mengambil perbedaannya.
sumber