kiểm tra ý nghĩa của tương tác trong các mô hình hỗn hợp tuyến tính trong nlme trong R

Question

Tôi sử dụng hàm lme trong gói nlme R để kiểm tra xem các mức hệ số items có tương tác đáng kể với các mức hệ số condition hay không. Hệ số condition có hai cấp độ: Control và Treatment và hệ số items có 3 cấp: E1,...,E3. Tôi sử dụng mã sau:

f.lme = lme(response ~ 0 + factor(condition) * factor(items), random = ~1|subject) 

trong đó subject là hiệu ứng ngẫu nhiên. Bằng cách này, khi tôi chạy:

summary(f.lme)$tTable 

tôi sẽ nhận được kết quả như sau:

factor(condition)Control 
factor(condition)Treatment 
factor(items)E2 
factor(items)E3 
factor(condition)Treatment:factor(items)E2 
factor(condition)Treatment:factor(items)E3 

cùng với Value, Std.Error, DF, t-value, p-value cột. Tôi có hai câu hỏi:

1. Nếu tôi muốn so sánh Control vs Treatment thì tôi chỉ sử dụng estimable() chức năng trong gmodels và làm cho một sự tương phản của (-1,1,0,0,0,0)?

2. Tôi quan tâm đến mức độ biểu lộ items, tức là E1, E2, E3 khác nhau trên khắp condition, vì vậy tôi đang quan tâm đến việc liệu các điều khoản tương tác có ý nghĩa (bằng cách chỉ kiểm tra cột p-value ??):

   factor(condition)Treatment:factor(items)E2 factor(condition)Treatment:factor(items)E3

Tuy nhiên, làm cách nào để biết liệu factor(condition)Treatment:factor(items)E1 có quan trọng hay không? Nó không được hiển thị trong bản tóm tắt đầu ra và tôi nghĩ rằng nó có một cái gì đó để làm với độ tương phản được sử dụng trong R ... Cảm ơn rất nhiều!

Answer 1

Tôi trân trọng không đồng ý với @ sven-hohenstein

Trong R, mã hóa mặc định cho các biến categorial là xử lý độ tương phản mã hóa. Trong điều trị tương phản, mức đầu tiên là mức tham chiếu. Tất cả các mức yếu tố còn lại được so sánh với mức tham chiếu.

Đầu tiên, các hiệu ứng cố định được chỉ định ở đây với điểm chặn bằng không, ... ~ 0 + .... Điều này có nghĩa là mã hóa condition không còn là contr.treatment. Nếu tôi không nhầm, những tác động chính của Control và Treatment hiện nay interpretable như độ lệch tương ứng của họ khỏi nhóm có nghĩa là ...

Trong mô hình của bạn, các mục yếu tố có ba cấp độ: E1, E2, và E3. Hai sự tương phản kiểm tra sự khác nhau giữa (a) E2 và E1, và (b) E3 và E1. Các hiệu ứng chính của những tương phản này được ước tính cho việc kiểm soát mức độ của điều kiện yếu tố, vì đây là danh mục tham chiếu của yếu tố này.

... khi giá trị của items ở mức tham chiếu là E1!Do đó:

• hiệu lực thi hành chính Control = bao nhiêu Control:E1 quan sát đi chệch khỏi giá trị trung bình của mặt hàng E1.
• Hiệu ứng chính Treatment = số tiền quan sát Treatment:E1 chênh lệch giá trị trung bình của mặt hàng E1.
• Hiệu ứng chính E2 = bao nhiêu Control:E2 quan sát lệch khỏi giá trị của mục E2.
• Hiệu ứng chính E3 = số lượng Control quan sát lệch khỏi giá trị trung bình của mặt hàng E3.
• Interaction Treatment:E2 = bao nhiêu Treatment:E2 quan sát đi chệch khỏi giá trị trung bình của mặt hàng E2
• Interaction Treatment:E3 = bao nhiêu Treatment:E3 quan sát đi chệch khỏi giá trị trung bình của mặt hàng E3.

Cảm ơn con trỏ đến estimable, tôi chưa thử trước đây. Đối với các tương phản tùy chỉnh, tôi đã là (ab)using glht từ gói multcomp.

Answer 2

Cách thông thường để kiểm tra xem tương tác có ý nghĩa hay không là thực hiện kiểm tra tỷ lệ khả năng (ví dụ: xem discussion on R-Sig-ME).

Để làm điều đó bạn có cũng để ước lượng một mô hình mà không cần sự tương tác và bạn cũng sẽ phải sử dụng method="ML":

f0 = lme(response ~ 0 + factor(condition) * factor(items), 
    random = ~1|subject, method="ML") 
f1 = lme(response ~ 0 + factor(condition) + factor(items), 
    random = ~1|subject, method="ML") 

Sau đó bạn có thể so sánh bằng anova:

anova(f0, f1) 

Xem thêm this blog post

Answer 3

Bạn cần giải quyết câu hỏi thứ hai của mình về tương tác trước. Bạn chắc chắn có thể thiết lập kiểm tra tỷ lệ khả năng như trong câu trả lời của Jan van der Laan. Bạn cũng có thể sử dụng trực tiếp anova trên đối tượng lme được trang bị. Xem trang trợ giúp cho anova.lme để biết thêm thông tin.

Về mặt diễn giải các hệ số của mình, tôi thường thấy rằng đôi khi mất một bảng tóm tắt của nhóm có nghĩa là để tìm ra kết hợp tuyến tính của các hệ số trong mô hình đại diện cho từng nhóm. Tôi sẽ trình bày một ví dụ với việc loại bỏ đánh chặn như trong câu hỏi của bạn, mặc dù tôi thấy điều này hiếm khi giúp tôi tìm ra các hệ số của tôi khi tôi có hai yếu tố trong một mô hình. Dưới đây là ví dụ về ý nghĩa của tôi với tập dữ liệu Orthodont (mà tôi quyết định cân bằng):

require(nlme) 

# Make dataset balanced 
Orthodont2 = Orthodont[-c(45:64),] 
# Factor age 
Orthodont2$fage = factor(Orthodont2$age) 

# Create a model with an interaction using lme; remove the intercept 
fit1 = lme(distance ~ Sex*fage - 1, random = ~1|Subject, data = Orthodont2) 
summary(fit1) 

Dưới đây là các hiệu ứng cố định ước tính. Nhưng mỗi hệ số này đại diện cho cái gì?

Fixed effects: distance ~ Sex * fage - 1 
        Value Std.Error DF t-value p-value 
SexMale   23.636364 0.7108225 20 33.25213 0.0000 
SexFemale  21.181818 0.7108225 20 29.79903 0.0000 
fage10   0.136364 0.5283622 61 0.25809 0.7972 
fage12   2.409091 0.5283622 61 4.55954 0.0000 
fage14   3.727273 0.5283622 61 7.05439 0.0000 
SexFemale:fage10 0.909091 0.7472171 61 1.21664 0.2284 
SexFemale:fage12 -0.500000 0.7472171 61 -0.66915 0.5059 
SexFemale:fage14 -0.818182 0.7472171 61 -1.09497 0.2778 

Tóm tắt các phương tiện nhóm giúp tìm ra điều này.

require(plyr) 
ddply(Orthodont2, .(Sex, age), summarise, dist = mean(distance)) 

    Sex fage  dist 
1 Male 8 23.63636 
2 Male 10 23.77273 
3 Male 12 26.04545 
4 Male 14 27.36364 
5 Female 8 21.18182 
6 Female 10 22.22727 
7 Female 12 23.09091 
8 Female 14 24.09091 

Lưu ý rằng hệ số hiệu ứng cố định đầu tiên, được gọi là SexMale, là khoảng cách trung bình cho nam 8 tuổi. Hiệu ứng cố định SexFemale là khoảng cách trung bình của nữ 8 tuổi.Đó là những thứ dễ thấy nhất (tôi luôn bắt đầu với những cái dễ), nhưng phần còn lại không quá tệ để tìm ra. Khoảng cách trung bình của nam 10 tuổi là hệ số đầu tiên cộng với hệ số thứ ba (fage10). Khoảng cách trung bình cho nữ 10 tuổi là tổng của các hệ số SexFemale, fage10 và SexFemale:fage10. Phần còn lại theo cùng các dòng.

Khi bạn biết cách tạo kết hợp tuyến tính của các hệ số cho nhóm, bạn có thể sử dụng estimable để tính toán bất kỳ sự so sánh nào. Tất nhiên có một loạt các caveats ở đây về hiệu ứng chính, bằng chứng thống kê của một tương tác, lý do lý do để lại tương tác trong, vv Đó là để bạn có thể quyết định. Nhưng nếu tôi rời khỏi tương tác trong mô hình (lưu ý không có bằng chứng thống kê về tương tác, xem anova(fit1)) và muốn so sánh giá trị trung bình tổng thể của Male đến Female, tôi sẽ ghi lại các kết hợp tuyến tính sau đây:

# male/age group means 
male8 = c(1, 0, 0, 0, 0, 0, 0, 0) 
male10 = c(1, 0, 1, 0, 0, 0, 0, 0) 
male12 = c(1, 0, 0, 1, 0, 0, 0, 0) 
male14 = c(1, 0, 0, 0, 1, 0, 0, 0) 

# female/age group means 
female8 = c(0, 1, 0, 0, 0, 0, 0, 0) 
female10 = c(0, 1, 1, 0, 0, 1, 0, 0) 
female12 = c(0, 1, 0, 1, 0, 0, 1, 0) 
female14 = c(0, 1, 0, 0, 1, 0, 0, 1) 

# overall male group mean 
male = (male8 + male10 + male12 +male14)/4 
# overall female group mean 
female = (female8 + female10 + female12 + female14)/4 

require(gmodels) 
estimable(fit1, rbind(male - female)) 

Bạn có thể kiểm tra nhóm tổng thể của mình có nghĩa là đảm bảo bạn đã thực hiện kết hợp tuyến tính của các hệ số một cách chính xác.

ddply(Orthodont2, .(Sex), summarise, dist = mean(distance)) 

    Sex  dist 
1 Male 25.20455 
2 Female 22.64773