Apa yang dilakukan perintah anova () dengan objek model lmer?

Mudah-mudahan ini adalah pertanyaan yang seseorang di sini dapat menjawab untuk saya tentang sifat dekomposisi jumlah kuadrat dari model efek campuran cocok dengan lmer(dari paket lme4 R).

Pertama saya harus mengatakan bahwa saya menyadari kontroversi dengan menggunakan pendekatan ini, dan dalam praktiknya saya akan lebih cenderung menggunakan LRT bootstrap untuk membandingkan model (seperti yang disarankan oleh Faraway, 2006). Namun, saya bingung bagaimana mereplikasi hasil, dan jadi untuk kewarasan saya sendiri saya pikir saya akan bertanya di sini.

Pada dasarnya, saya mulai terbiasa menggunakan model efek campuran yang sesuai dengan lme4paket. Saya tahu bahwa Anda dapat menggunakan anova()perintah untuk memberikan ringkasan pengujian efek tetap dalam model secara berurutan. Sejauh yang saya tahu ini adalah apa yang Faraway (2006) sebut sebagai pendekatan 'kotak Rata-rata yang diharapkan'. Yang ingin saya ketahui adalah bagaimana jumlah kuadrat dihitung?

Saya tahu bahwa saya dapat mengambil nilai estimasi dari model tertentu (menggunakan coef()), berasumsi bahwa mereka sudah diperbaiki, dan kemudian melakukan tes menggunakan jumlah kuadrat residual model dengan dan tanpa faktor-faktor yang menarik. Ini bagus untuk model yang mengandung faktor dalam-subjek tunggal. Namun, ketika menerapkan desain petak-petak, jumlah nilai kuadrat yang saya dapatkan setara dengan nilai yang dihasilkan oleh R aov()dengan Error()peruntukan yang sesuai . Namun, ini tidak sama dengan jumlah kuadrat yang dihasilkan oleh anova()perintah pada objek model, terlepas dari kenyataan bahwa rasio-F adalah sama.

Tentu saja ini masuk akal karena tidak perlu untuk Error()strata dalam model campuran. Namun, ini harus berarti bahwa jumlah kuadrat dihukum entah bagaimana dalam model campuran untuk memberikan rasio-F yang tepat. Bagaimana ini dicapai? Dan bagaimana model tersebut entah bagaimana mengoreksi jumlah kuadrat antar-plot tetapi tidak mengoreksi jumlah kuadrat dalam-plot. Jelas ini adalah sesuatu yang diperlukan untuk ANOVA split-plot klasik yang dicapai dengan merancang nilai kesalahan yang berbeda untuk efek yang berbeda, jadi bagaimana model efek campuran memungkinkan untuk ini?

Pada dasarnya, saya ingin dapat mereplikasi hasil dari anova()perintah yang diterapkan pada objek model lmer sendiri untuk memverifikasi hasil dan pemahaman saya, namun, saat ini saya dapat mencapai ini untuk desain dalam-subjek normal tetapi tidak untuk pemisahan. desain plot dan sepertinya saya tidak tahu mengapa ini terjadi.

Sebagai contoh:

library(faraway)
library(lme4)
data(irrigation)

anova(lmer(yield ~ irrigation + variety + (1|field), data = irrigation))

Analysis of Variance Table
           Df Sum Sq Mean Sq F value
irrigation  3 1.6605  0.5535  0.3882
variety     1 2.2500  2.2500  1.5782

summary(aov(yield ~ irrigation + variety + Error(field/irrigation), data = irrigation))

Error: field
           Df Sum Sq Mean Sq F value Pr(>F)
irrigation  3  40.19   13.40   0.388  0.769
Residuals   4 138.03   34.51               

Error: Within
          Df Sum Sq Mean Sq F value Pr(>F)
variety    1   2.25   2.250   1.578  0.249
Residuals  7   9.98   1.426               

Seperti dapat dilihat di atas, semua rasio-F setuju. Jumlah kuadrat untuk varietas juga setuju. Namun, jumlah kuadrat untuk irigasi tidak setuju, namun tampaknya output yang lebih kecil diskalakan. Jadi, apa sebenarnya yang dilakukan perintah anova ()?

Gunakan Sumbernya, Luke. Kita dapat mengintip ke dalam fungsi ANOVA dengan melakukan getAnywhere(anova.Mermod). Bagian pertama dari fungsi itu adalah untuk membandingkan dua model yang berbeda. Anova pada efek tetap datang di elseblok besar di babak kedua:

 dc <- getME(object, "devcomp")
        X <- getME(object, "X")
        asgn <- attr(X, "assign")
        stopifnot(length(asgn) == (p <- dc$dims[["p"]]))
            ss <- as.vector(object@pp$RX() %*% object@beta)^2
        names(ss) <- colnames(X)
        terms <- terms(object)
        nmeffects <- attr(terms, "term.labels")[unique(asgn)]
        if ("(Intercept)" %in% names(ss)) 
            nmeffects <- c("(Intercept)", nmeffects)
        ss <- unlist(lapply(split(ss, asgn), sum))
        stopifnot(length(ss) == length(nmeffects))
        df <- vapply(split(asgn, asgn), length, 1L)
        ms <- ss/df
        f <- ms/(sigma(object)^2)
        table <- data.frame(df, ss, ms, f)
        dimnames(table) <- list(nmeffects, c("Df", "Sum Sq", 
            "Mean Sq", "F value"))
        if ("(Intercept)" %in% nmeffects) 
            table <- table[-match("(Intercept)", nmeffects), 
                ]
        attr(table, "heading") <- "Analysis of Variance Table"
        class(table) <- c("anova", "data.frame")
        table

objectadalah output lmer. Kita mulai menghitung jumlah kuadrat di baris 5: ss <- as.vector ...Kode ini mengalikan parameter tetap (in beta) dengan matriks segitiga atas; lalu kuadratkan setiap istilah. Inilah matriks segitiga atas untuk contoh irigasi. Setiap baris sesuai dengan salah satu dari lima parameter efek tetap (intersep, 3 derajat kebebasan untuk irigasi, 1 df untuk variasi).

zapsmall(irrigation.lmer@pp$RX(), digits = 3)
      [,1]  [,2]   [,3]   [,4]  [,5]
[1,] 0.813 0.203  0.203  0.203 0.407
[2,] 0.000 0.352 -0.117 -0.117 0.000
[3,] 0.000 0.000  0.332 -0.166 0.000
[4,] 0.000 0.000  0.000  0.287 0.000
[5,] 0.000 0.000  0.000  0.000 2.000

Baris pertama memberi Anda jumlah kuadrat untuk intersep dan yang terakhir memberi Anda SS untuk efek variasi dalam-bidang. Baris 2-4 hanya melibatkan 3 parameter untuk tingkat irigasi, jadi pra-penggandaan memberi Anda tiga potong SS untuk irigasi.

Potongan-potongan ini tidak menarik dalam diri mereka sendiri karena mereka berasal dari kontras pengobatan standar dalam R, tetapi sejalan ss <- unlist(lapply(split ....Bates mengambil bit dari jumlah kuadrat sesuai dengan jumlah level dan faktor mana yang mereka rujuk. Ada banyak pembukuan yang terjadi di sini. Kami mendapatkan derajat kebebasan juga (yaitu 3 untuk irigasi). Kemudian, ia mendapatkan kotak rata-rata yang muncul di cetakan anova. Akhirnya, ia membagi semua kuadrat berarti dengan varians residual dalam kelompok sigma(object)^2,.

lmeraovlmerRX$R_{00}$$\sqrt{\sigma^2/\sigma_f^2}$$\sigma_f^2$

Asimptotik, perkiraan efek tetap memiliki distribusi:

$R_{00}\hat{\beta}$$\beta=0$$\sigma^2$$\sigma^2$$\sigma^2$$R_{00}$$\sigma^2$

Perhatikan bahwa Anda tidak akan mendapatkan statistik F yang sama seandainya data tidak seimbang. Anda juga tidak akan memperoleh statistik F yang sama jika Anda telah menggunakan ML dan bukan REML.

aov$\sigma^2$$\sigma_f^2$$\sigma^2$$\sigma_f^2$

Menariknya, Bates dan Pinheiro merekomendasikan penggunaan ANOVA pada pemasangan dua model dan melakukan uji rasio kemungkinan. Yang terakhir cenderung anti-konservatif.

$R_{00}$

zapsmall(fit2@pp$RX(), digits = 3)
      [,1]  [,2]   [,3]   [,4]   [,5]
[1,] 0.816 0.205  0.205  0.205  0.457
[2,] 0.000 0.354 -0.119 -0.119 -0.029
[3,] 0.000 0.000  0.334 -0.168 -0.040
[4,] 0.000 0.000  0.000  0.288 -0.071
[5,] 0.000 0.000  0.000  0.000  1.874

Seperti yang Anda lihat, jumlah kuadrat untuk parameter irigasi sekarang mengandung beberapa varietyefek juga.