phân vùng của anova và so sánh (orthogonal single df) trong r

Question

Tôi muốn thực hiện tương phản trực giao đơn df trong anova (mô hình cố định hoặc hỗn hợp). Đây chỉ là ví dụ:

require(nlme) 
data (Alfalfa) 
    Variety: a factor with levels Cossack, Ladak, and Ranger 
    Date : a factor with levels None S1 S20 O7 
    Block: a factor with levels 1 2 3 4 5 6 
    Yield : a numeric vector 

Những dữ liệu này được mô tả trong Snedecor và Cochran (1980) là một ví dụ của một thiết kế tách âm mưu. Cấu trúc xử lý được sử dụng trong thí nghiệm là giai thừa 3 lần, với ba giống cỏ linh lăng và bốn ngày ngày cắt thứ ba vào năm 1943. Các đơn vị thí nghiệm được bố trí thành sáu khối, mỗi khối được chia thành 4 lô. Các giống cỏ linh lăng (Cossac, Ladak, và Ranger) được gán ngẫu nhiên vào các khối và ngày cắt thứ ba (Không, S1 - 1 tháng 9, S20 - 20 tháng 9, và O7 — ngày 7 tháng 10) được gán ngẫu nhiên vào cốt truyện. Tất cả bốn ngày được sử dụng trên mỗi khối.

model<-with (Alfalfa, aov(Yield~Variety*Date +Error(Block/Date/Variety))) 

    > summary(model) 

Error: Block 
      Df Sum Sq Mean Sq F value Pr(>F) 
Residuals 5 4.15 0.83 

Error: Block:Date 
      Df Sum Sq Mean Sq F value Pr(>F) 
Date  3 1.9625 0.6542 17.84 3.29e-05 *** 
Residuals 15 0.5501 0.0367 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Error: Block:Date:Variety 
      Df Sum Sq Mean Sq F value Pr(>F) 
Variety  2 0.1780 0.08901 1.719 0.192 
Variety:Date 6 0.2106 0.03509 0.678 0.668 
Residuals 40 2.0708 0.05177 

Tôi muốn thực hiện một số so sánh (tương phản trực giao trong một nhóm), ví dụ cho ngày, hai sự tương phản:

(a) S1 vs others (S20 O7) 
    (b) S20 vs 07, 

Đối với yếu tố đa dạng hai tương phản:

(c) Cossack vs others (Ladak and Ranger) 
    (d) Ladak vs Ranger 

Vì vậy đầu ra anova sẽ trông giống như:

Error: Block 
       Df Sum Sq Mean Sq F value Pr(>F) 
Residuals 5 4.15 0.83 

Error: Block:Date 
      Df Sum Sq Mean Sq F value Pr(>F) 
Date  3 1.9625 0.6542 17.84 3.29e-05 *** 
     (a) S1 vs others ?  ? 
     (b) S20 vs 07  ?  ? 
Residuals 15 0.5501 0.0367 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Error: Block:Date:Variety 
      Df Sum Sq Mean Sq F value Pr(>F) 
Variety  2 0.1780 0.08901 1.719 0.192 
    (c) Cossack vs others ?  ? ? 
    (d) Ladak vs Ranger ?  ?  ? 
Variety:Date 6 0.2106 0.03509 0.678 0.668 
Residuals 40 2.0708 0.05177 

Tôi có thể thực hiện điều này bằng cách nào? ....................

Answer 1

Trước hết, tại sao lại sử dụng ANOVA? Bạn có thể sử dụng lme từ gói nlme và ngoài các thử nghiệm giả thuyết aov cung cấp cho bạn, bạn cũng nhận được ước tính có thể diễn giải về kích thước hiệu ứng và hướng của các hiệu ứng. Ở bất kỳ mức nào, hai phương pháp tiếp cận đến với tâm trí:

• Chỉ định tương phản trên các biến theo cách thủ công, như được giải thích here.
• Cài đặt gói multcomp và sử dụng glht.

glht là một chút ý kiến ​​về các mô hình đa biến trong dự đoán của chúng. Mặc dù vậy, nếu bạn tạo một ma trận đường chéo cm0 với cùng kích thước và tên mờ như vcov của mô hình của bạn (giả sử nó là lme phù hợp gọi là model0), thì summary(glht(model0,linfct=cm0)) phải cung cấp cùng một ước tính, SE và thử nghiệm thống kê dưới dạng summary(model0)$tTable (nhưng không chính xác giá trị p). Bây giờ, nếu bạn lộn xộn với các kết hợp tuyến tính của các hàng từ cm0 và tạo các ma trận mới với cùng số cột như cm0 nhưng các kết hợp tuyến tính này là hàng, cuối cùng bạn sẽ tìm ra mẫu để tạo ma trận. ước tính cho mỗi ô (kiểm tra nó với predict(model0,level=0)). Bây giờ, một ma trận khác với các hàng khác nhau của ma trận này sẽ cho bạn sự khác biệt giữa các nhóm tương ứng. Cách tiếp cận tương tự nhưng với các giá trị số được đặt thành 1 thay vì 0 có thể được sử dụng để lấy ước tính độ dốc cho mỗi ô. Sau đó, sự khác biệt giữa các ước tính độ dốc này có thể được sử dụng để có sự khác biệt về độ dốc giữa các nhóm.

Ba điều cần lưu ý:

• Như tôi đã nói các p-giá trị sẽ là sai lầm cho các mô hình khác hơn lm, (có thể, đã không cố gắng) aov, và các mô hình tồn tại nhất định. Điều này là do glht giả sử phân phối z thay vì phân phối t theo mặc định (ngoại trừ lm). Để có được đúng p-giá trị, mất kiểm định thống kê glht tính toán và tự làm 2*pt(abs(STAT),df=DF,lower=F) để có được hai đuôi p-giá trị nơi STAT là kiểm định thống kê được trả về bởi glht và DF là số df từ loại tương phản mặc định tương ứng trong summary(model0)$tTable.
• Tương phản của bạn có thể không còn thử nghiệm các giả thuyết độc lập nữa và cần phải hiệu chỉnh nhiều lần kiểm tra nếu chưa thực hiện. Chạy các giá trị p thông qua p.adjust.
• Đây là chưng cất của riêng tôi về rất nhiều sự trao tay từ các giáo sư và đồng nghiệp, và rất nhiều đọc Crossvalidated và Stackoverflow về các chủ đề liên quan. Tôi có thể sai theo nhiều cách, và nếu tôi, hy vọng một người nào đó hiểu biết hơn sẽ sửa cả hai chúng tôi.