Menguji asumsi normalitas untuk tindakan berulang anova? (dalam R)

Jadi dengan asumsi bahwa ada titik dalam menguji asumsi normalitas untuk anova (lihat 1 dan 2 )

Bagaimana bisa diuji dalam R?

Saya berharap untuk melakukan sesuatu seperti:

## From Venables and Ripley (2002) p.165.
utils::data(npk, package="MASS")
npk.aovE <- aov(yield ~  N*P*K + Error(block), npk)
residuals(npk.aovE)
qqnorm(residuals(npk.aov))

Yang tidak bekerja, karena "residu" tidak memiliki metode (atau memprediksi, dalam hal ini) untuk kasus anova tindakan berulang.

Jadi apa yang harus dilakukan dalam kasus ini?

Bisakah residu diekstraksi dari model pas yang sama tanpa istilah Kesalahan? Saya tidak cukup akrab dengan literatur untuk mengetahui apakah ini valid atau tidak, terima kasih sebelumnya atas sarannya.

Anda mungkin tidak mendapatkan respons sederhana residuals(npk.aovE)tetapi itu tidak berarti tidak ada residu pada objek tersebut. Lakukan strdan lihat bahwa di dalam level masih ada residu. Saya akan membayangkan Anda paling tertarik pada level "Dalam"

> residuals(npk.aovE$Within)
          7           8           9          10          11          12 
 4.68058815  2.84725482  1.56432584 -5.46900749 -1.16900749 -3.90234083 
         13          14          15          16          17          18 
 5.08903669  1.28903669  0.35570336 -3.27762998 -4.19422371  1.80577629 
         19          20          21          22          23          24 
-3.12755705  0.03910962  2.60396981  1.13730314  2.77063648  4.63730314 

Pelatihan dan praktik saya sendiri belum menggunakan pengujian normal, melainkan menggunakan plot QQ dan pengujian paralel dengan metode yang kuat.

Opsi lain adalah menggunakan lmefungsi nlmepaket (dan kemudian meneruskan model yang diperoleh ke anova). Anda dapat menggunakan residualspada outputnya.

Saya pikir asumsi normalitas dapat dinilai untuk setiap langkah yang diulang, sebelum melakukan analisis. Saya akan membentuk kembali kerangka data sehingga setiap kolom sesuai dengan ukuran yang diulang, dan kemudian melakukan shapiro.test untuk masing-masing kolom tersebut.

apply(cast(melt(npk,measure.vars="yield"), ...~N+P+K)[-c(1:2)],2,function(x) shapiro.test(x)$p.value)

Venables dan Ripley menjelaskan bagaimana melakukan diagnosa residual untuk desain tindakan berulang di kemudian hari dalam buku mereka (hal. 284), di bagian efek acak dan campuran.

Fungsi residual (atau residu) diimplementasikan untuk hasil aov untuk setiap strata:

dari contoh mereka: oats.aov <- aov(Y ~ N + V + Error(B/V), data=oats, qr=T)

Untuk mendapatkan nilai atau residu yang dipasang:

"Jadi fitted(oats.aov[[4]])dan resid(oats.aov[[4]])merupakan vektor dengan panjang 54 yang mewakili nilai dan residu yang pas dari strata terakhir."

Yang penting, mereka menambahkan:

"Tidak mungkin untuk mengaitkan mereka secara unik dengan plot percobaan asli."

Untuk diagnosa, mereka menggunakan proyeksi:

plot(fitted(oats.aov[[4]]), studres(oats.aov[[4]]))
abline(h=0, lty=2)
oats.pr <- proj(oats.aov)
qqnorm(oats.pr[[4]][, "Residuals"], ylab = "Stratum 4 residuals")
qqline(oats.pr[[4]][, "Residuals"])

Mereka juga menunjukkan bahwa model dapat dilakukan dengan menggunakan lme, seperti yang diposting pengguna lain.

Berikut adalah dua opsi, dengan aov dan dengan lme (saya pikir yang ke-2 lebih disukai):

require(MASS) ## for oats data set
require(nlme) ## for lme()
require(multcomp) ## for multiple comparison stuff

Aov.mod <- aov(Y ~ N * V + Error(B/V), data = oats)
the_residuals <- aov.out.pr[[3]][, "Residuals"]

Lme.mod <- lme(Y ~ N * V, random = ~1 | B/V, data = oats)
the_residuals <- residuals(Lme.mod)

Contoh asli datang tanpa interaksi ( Lme.mod <- lme(Y ~ N * V, random = ~1 | B/V, data = oats)) tetapi tampaknya berfungsi dengan itu (dan menghasilkan hasil yang berbeda, sehingga ia melakukan sesuatu).

Dan itu saja ...

tapi untuk kelengkapan:

1 - Ringkasan model

summary(Aov.mod)
anova(Lme.mod)

2 - Tes Tukey dengan tindakan berulang anova (3 jam mencari ini !!).

summary(Lme.mod)
summary(glht(Lme.mod, linfct=mcp(V="Tukey")))

3 - Plot normalitas dan homoseksualitas

par(mfrow=c(1,2)) #add room for the rotated labels
aov.out.pr <- proj(aov.mod)                                            
#oats$resi <- aov.out.pr[[3]][, "Residuals"]
oats$resi <- residuals(Lme.mod)
qqnorm(oats$resi, main="Normal Q-Q") # A quantile normal plot - good for checking normality
qqline(oats$resi)
boxplot(resi ~ interaction(N,V), main="Homoscedasticity", 
        xlab = "Code Categories", ylab = "Residuals", border = "white", 
        data=oats)
points(resi ~ interaction(N,V), pch = 1, 
       main="Homoscedasticity",  data=oats)