我进行了一项实验,在该实验中,我抚养了来自两个不同来源人群的不同家庭。每个家庭都被分配了两种治疗方法之一。实验结束后,我测量了每个人的几个特征。为了测试治疗或来源的影响以及它们的相互作用,我使用了一个以家庭为随机因素的线性混合效应模型,即
lme(fixed=Trait~Treatment*Source,random=~1|Family,method="ML")
到目前为止一切顺利,现在我必须计算相对方差分量,即由处理或来源以及交互作用解释的变异百分比。
如果没有随机效应,我可以轻松地使用平方和 (SS) 来计算每个因素解释的方差。但是对于混合模型(使用 ML 估计),没有 SS,因此我认为我也可以使用 Treatment 和 Source 作为随机效应来估计方差,即
lme(fixed=Trait~1,random=~(Treatment*Source)|Family, method="REML")
但是,在某些情况下, lme 不会收敛,因此我使用了 lme4 包中的 lmer:
lmer(Trait~1+(Treatment*Source|Family),data=DATA)
我使用汇总函数从模型中提取方差:
model<-lmer(Trait~1+(Treatment*Source|Family),data=regrexpdat)
results<-VarCorr(model)
variances<-results[,3]
我得到与 VarCorr 函数相同的值。然后我使用这些值来计算变化的实际百分比,将总和作为总变化。
我正在努力解释初始 lme 模型的结果(治疗和来源作为固定效应)和随机模型来估计方差分量(治疗和来源作为随机效应)。我发现在大多数情况下,每个因素解释的方差百分比并不对应于固定效应的显着性。
例如对于特征 HD,最初的 lme 表明了相互作用的趋势以及治疗的意义。使用反向程序,我发现治疗具有接近显着的趋势。但是,在估计方差分量时,我发现 Source 的方差最高,占总方差的 26.7%。
我:
anova(lme(fixed=HD~as.factor(Treatment)*as.factor(Source),random=~1|as.factor(Family),method="ML",data=test),type="m")
numDF denDF F-value p-value
(Intercept) 1 426 0.044523 0.8330
as.factor(Treatment) 1 426 5.935189 0.0153
as.factor(Source) 1 11 0.042662 0.8401
as.factor(Treatment):as.factor(Source) 1 426 3.754112 0.0533
和 lmer:
summary(lmer(HD~1+(as.factor(Treatment)*as.factor(Source)|Family),data=regrexpdat))
Linear mixed model fit by REML
Formula: HD ~ 1 + (as.factor(Treatment) * as.factor(Source) | Family)
Data: regrexpdat
AIC BIC logLik deviance REMLdev
-103.5 -54.43 63.75 -132.5 -127.5
Random effects:
Groups Name Variance Std.Dev. Corr
Family (Intercept) 0.0113276 0.106431
as.factor(Treatment) 0.0063710 0.079819 0.405
as.factor(Source) 0.0235294 0.153393 -0.134 -0.157
as.factor(Treatment)L:as.factor(Source) 0.0076353 0.087380 -0.578 -0.589 -0.585
Residual 0.0394610 0.198648
Number of obs: 441, groups: Family, 13
Fixed effects:
Estimate Std. Error t value
(Intercept) -0.02740 0.03237 -0.846
因此我的问题是,我在做什么正确吗?或者我应该使用另一种方法来估计每个因素(即治疗、来源及其相互作用)解释的方差量。例如,效果大小会是更合适的方法吗?