在 R 中的 nlme 中测试线性混合模型中交互作用的显着性 [英] test for significance of interaction in linear mixed models in nlme in R

问题描述

我使用 nlme R 包中的 lme 函数来测试因子 items 的水平是否与因子 的水平有显着的交互作用条件.因素condition有两个级别:Control和Treatment，因素items有3个级别:E1,...,E3.我使用以下代码:

I use lme function in the nlme R package to test if levels of factor items has significant interaction with levels of factor condition. The factor condition has two levels: Control and Treatment, and the factor items has 3 levels: E1,...,E3. I use the following code:

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

其中 subject 是随机效应.这样，当我运行时:

where subject is the random effect. In this way, when I run:

summary(f.lme)$tTable

我会得到以下输出:

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

连同 Value、Std.Error、DF、t-value、p-value 列.我有两个问题:

together with Value, Std.Error, DF, t-value, p-value columns. I have two questions:

1. 如果我想比较 Control 和 Treatment，我应该在 gmodels 中使用 estimable() 函数吗? 和 (-1,1,0,0,0,0)?

1. If I want to compare Control vs. Treatment, shall I just use estimable() function in gmodels and make a contrast of (-1,1,0,0,0,0)?

我感兴趣的是items的等级，即E1, E2, E3在condition中是否不同，所以我是对交互项是否重要感兴趣(只需检查 p-value 列??):

I am interested in whether levels of items, i.e. E1, E2, E3 are different across condition, so I am interested in whether the interaction terms are significant (by just checking the p-value column??):

factor(condition)处理:factor(items)E2factor(condition)Treatment:factor(items)E3

但是，我如何判断 factor(condition)Treatment:factor(items)E1 是否重要?它没有显示在摘要输出中，我认为这与 R 中使用的对比度有关......非常感谢！

However, how can I tell if factor(condition)Treatment:factor(items)E1 is significant or not? It is not shown in the summary output and I think it has something to do with the contrast used in R... Thanks a lot!

推荐答案

您需要先解决关于交互的第二个问题.您当然可以按照 Jan van der Laan 的回答设置似然比检验.您也可以直接在拟合的 lme 对象上使用 anova.有关详细信息，请参阅 anova.lme 的帮助页面.

You'll need to address your second question about the interaction first. You can certainly set up the likelihood ratio test as in Jan van der Laan's answer. You can also use anova directly on a fitted lme object. See the help page for anova.lme for more information.

在解释系数方面，我经常发现有时需要制作组均值的汇总表，以便适当地找出模型中系数的哪些线性组合代表每个组.我将展示一个像您的问题中那样去除截距的例子，尽管我发现一旦我在模型中有两个因素，这很少能帮助我计算出我的系数.以下是我对 Orthodont 数据集(我决定对其进行平衡)的含义的示例:

In terms of interpreting your coefficients, I often find that it sometimes takes making a summary table of the group means in order to appropriately figure out which linear combination of coefficients in the model represents each group. I will show an example with the intercept removal like in your question, although I find this rarely helps me figure out my coefficients once I have two factors in a model. Here is an example of what I mean with the Orthodont dataset (which I decided to make balanced):

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)

这是估计的固定效应.但这些系数分别代表什么?

Here are the estimated fixed effects. But what do each of these coefficients represent?

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

总结组均值有助于解决这个问题.

Making a summary of group means helps figure this out.

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

请注意，第一个固定效应系数，称为 SexMale，是 8 岁男性的平均距离.固定效应 SexFemale 是 8 岁女性的平均距离.这些是最容易看到的(我总是从简单的开始)，但其余的都不太好弄明白.10 岁男性的平均距离是第一个系数加上第三个系数 (fage10).10 岁女性的平均距离是系数 SexFemale、fage10 和 SexFemale:fage10 的总和.其余的都遵循相同的路线.

Notice that the first fixed effect coefficient, called SexMale, is the mean distance for the age 8 males. The fixed effect SexFemale is the age 8 female mean distance. Those are the easiest to see (I always start with the easy ones), but the rest aren't too bad to figure out. The mean distance for age 10 males is the first coefficient plus the third coefficient (fage10). The mean distance for age 10 females is the sum of the coefficients SexFemale, fage10, and SexFemale:fage10. The rest follow along the same lines.

一旦您知道如何为组均值创建系数的线性组合，您就可以使用 estimable 来计算任何感兴趣的比较.当然，这里有很多关于主效应、相互作用的统计证据、保留相互作用的理论原因等的警告.这由你来决定.但是，如果我将交互留在模型中(请注意，没有交互的统计证据，请参阅 anova(fit1))，但仍想比较 Male 的总体平均值> 到 Female，我会写出以下系数的线性组合:

Once you know how to create linear combinations of coefficients for the group means, you can use estimable to calculate any comparisons of interest. Of course there are a bunch of caveats here about main effects, statistical evidence of an interactions, theoretical reasons to leave interactions in, etc. That is for you to decide. But if I was leaving the interaction in the model (note there is no statistical evidence of an interaction, see anova(fit1)) and yet wanted to compare the overall mean of Male to Female, I would write out the following linear combinations of coefficients:

# 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))

您可以检查您的整体组均值，以确保您正确地进行了系数的线性组合.

You can check your overall group means to make sure you made your linear combinations of coefficients correctly.

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

     Sex     dist
1   Male 25.20455
2 Female 22.64773

这篇关于在 R 中的 nlme 中测试线性混合模型中交互作用的显着性的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！