10

我使用R 包lme中的函数来测试 factor 的水平是否与 factor的水平有显着的交互作用。该因子有两个水平：和，该因子有 3 个水平：。我使用以下代码：nlmeitemsconditionconditionControlTreatmentitemsE1,...,E3

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

哪里subject是随机效应。这样，当我运行时：

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列。我有两个问题：

  1. 如果我想比较Controlvs. Treatment，我应该只使用estimable()函数gmodels并进行对比(-1,1,0,0,0,0)吗？

  2. 我对 的水平是否不同感兴趣items，即在 之间是否E1, E2, E3不同condition，所以我对交互项是否显着感兴趣（只需检查p-value列？？）：

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

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

4

3 回答 3

5

您需要首先解决有关交互的第二个问题。您当然可以像 Jan van der Laan 的回答那样设置似然比检验。您也可以anova直接在适合的 lme 对象上使用。anova.lme有关更多信息，请参阅帮助页面。

在解释您的系数方面，我经常发现有时需要制作组均值的汇总表，以便适当地确定模型中哪个系数的线性组合代表每个组。我将在你的问题中展示一个截距去除的例子，尽管我发现一旦我在模型中有两个因素，这很少能帮助我计算出我的系数。这是我对 Orthodont 数据集（我决定平衡）的意思的一个例子：

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)

这是估计的固定效应。但是这些系数中的每一个代表什么？

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

总结组均值有助于弄清楚这一点。

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。其余的遵循相同的路线。

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

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

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

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

     Sex     dist
1   Male 25.20455
2 Female 22.64773
于 2013-07-31T16:38:45.993 回答
5

我恭敬地不同意@sven-hohenstein

在 R 中，分类变量的默认编码是治疗对比编码。在治疗对比中，第一水平是参考水平。所有剩余的因子水平都与参考水平进行比较。

首先，这里用零截距指定固定效应，... ~ 0 + ...。这意味着condition编码不再contr.treatment。如果我没记错的话， 和 的主要影响Control现在Treatment可以解释为它们各自与组的偏差意味着......

在您的模型中，因子项具有三个级别：E1、E2 和 E3。这两个对比测试了 (a) E2 和 E1 以及 (b) E3 和 E1 之间的差异。这些对比的主要影响是针对因子条件的控制水平估计的，因为这是该因子的参考类别。

...当 的值items处于其参考水平时E1! 所以：

  • 主效应Control=Control:E1观察值偏离 item 平均值的程度E1。
  • 主效应Treatment=Treatment:E1观察值偏离 item 平均值的程度E1。
  • 主效应E2=Control:E2观察值偏离 item 平均值的程度E2。
  • 主效应E3=Control观察值偏离 item 平均值的程度E3。
  • 交互作用Treatment:E2=Treatment:E2观察值偏离项目均值的程度E2
  • 交互作用Treatment:E3=Treatment:E3观察值偏离 item 平均值的程度E3。

感谢您的指针estimable，我以前没有尝试过。对于自定义对比，我一直（ab）glht从multcomp包装中使用。

于 2013-07-25T20:59:45.153 回答
3

测试交互是否显着的常用方法是进行似然比测试（例如，参见关于 R-Sig-ME 的讨论）。

为此，您还必须估计一个没有交互的模型，您还必须使用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")

然后，您可以使用以下方法进行比较anova：

anova(f0, f1)

另请参阅此博客文章

于 2013-07-31T09:50:06.413 回答