我使用R包lme中的函数nlme来测试因子的水平是否与因子水平items有显着的相互作用condition.该因子condition有两个级别:Control和Treatment,因子items有3个级别:E1,...,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列一起.我有两个问题:
如果我想比较Control与Treatment应我只是使用estimable()功能gmodels,使的反差(-1,1,0,0,0,0)？
我感兴趣的水平是否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中使用的对比度有关...非常感谢!
我恭敬地不同意@ sven-hohenstein
在R中,分类变量的默认编码是治疗对比度编码.在治疗对比中,第一级是参考水平.将所有剩余因子水平与参考水平进行比较.
首先,这里用零截距指定固定效果... ~ 0 + ....这意味着condition不再编码contr.treatment.如果我没有记错,主要作用Control和Treatment现在解释从组各自偏差的意思是...
在您的模型中,因子项有三个级别:E1,E2和E3.两个对比测试了(a)E2和E1,以及(b)E3和E1之间的差异.这些对比的主要影响是对因子条件的控制水平进行估计,因为这是该因子的参考类别.
......当价值items处于其参考水平时E1!因此:
Control= Control:E1观察偏离项目平均值的程度E1.Treatment= Treatment:E1观察偏离项目平均值的程度E1.E2= Control:E2观察偏离项目平均值的程度E2.E3= Control观察偏离项目平均值的程度E3.Treatment:E2= Treatment:E2观察偏离项目平均值的程度E2Treatment:E3= Treatment:E3观察偏离项目平均值的程度E3.感谢指针estimable,我以前没试过.对于自定义的对比,我一直使用(AB)glht从multcomp包.
您需要先解决有关互动的第二个问题.您当然可以像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
请注意,第一个固定效应系数称为SexMale8岁男性的平均距离.固定效应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