混合效应模型的多重比较

机器算法验证 r 方差分析 混合模式 多重比较 重复测量
2022-02-07 19:52:56

我正在尝试使用混合效应模型分析一些数据。我收集的数据代表了一些不同基因型的年轻动物随时间的体重。

我正在使用这里提出的方法： https ://gribblelab.wordpress.com/2009/03/09/repeated-measures-anova-using-r/

特别是我正在使用解决方案＃2

所以我有类似的东西

require(nlme)
model <- lme(weight ~ time * Genotype, random = ~1|Animal/time, 
         data=weights)    
av <- anova(model)

现在，我想进行一些多重比较。使用multcomp我可以做到：

require(multcomp)
comp.geno <- glht(model, linfct=mcp(Genotype="Tukey"))
print(summary(comp.geno))

当然，我也可以用时间做同样的事情。

我有两个问题：

  1. 我如何使用mcp来查看时间和基因型之间的相互作用？

  2. 当我运行时，glht我收到此警告：

    covariate interactions found -- default contrast might be inappropriate

    这是什么意思？我可以放心地忽略它吗？或者我应该怎么做才能避免它？

编辑： 我发现这个 PDF说：

因为在这种情况下不可能自动确定感兴趣的参数，所以 multcomp 中的 mcp() 默认情况下将只生成主效应的比较，而忽略协变量和交互作用。从 1.1-2 版本开始，可以分别使用参数 interaction_average = TRUE 和 covariate_average = TRUE 指定对交互项和协变量进行平均，而早于 1.0-0 的版本会自动对交互项进行平均。然而，我们建议用户手动写出他们想要的对比集。每当对默认对比测量值存在疑问时，都应该这样做，这通常发生在具有更高阶交互项的模型中。我们参考 Hsu (1996), Chapter~7, 和 Searle (1971), Chapter~7.3，关于这个问题的进一步讨论和例子。

我无权访问这些书，但也许这里有人有？

1个回答

如果time和Genotype看起来都是分类预测变量，并且您有兴趣将所有时间/基因型对相互比较，那么您可以只创建一个交互变量，并在其上使用 Tukey 对比：

weights$TimeGeno <- interaction(weigths$Time, weights$Geno)
model <- lme(weight ~ TimeGeno, random = ~1|Animal/time, data=weights) 
comp.timegeno <- glht(model, linfct=mcp(TimeGeno="Tukey"))

如果您对其他对比感兴趣，那么您可以使用linfct参数可以为对比采用系数矩阵这一事实 - 这样您就可以准确地设置您想要的比较。

编辑

评论中出现了一些担忧，即装有预测器的模型与装有TimeGeno预测器的原始模型不同Time * Genotype。情况并非如此，模型是等效的。唯一的区别在于固定效果的参数化，它的设置是为了更容易使用该glht功能。

我使用了一个内置数据集（它有饮食而不是基因型）来证明这两种方法具有相同的可能性、预测值等：

> # extract a subset of a built-in dataset for the example
> data(BodyWeight)
> ex <- as.data.frame(subset(BodyWeight, Time %in% c(1, 22, 44)))
> ex$Time <- factor(ex$Time)
> 
> #create interaction variable
> ex$TimeDiet <- interaction(ex$Time, ex$Diet)
    > 
    > model1 <- lme(weight ~ Time * Diet, random = ~1|Rat/Time,  data=ex)    
    > model2 <- lme(weight ~ TimeDiet, random = ~1|Rat/Time, data=ex)    
    > 
    > # the degrees of freedom, AIC, BIC, log-likelihood are all the same 
    > anova(model1, model2)
           Model df      AIC      BIC    logLik
    model1     1 12 367.4266 387.3893 -171.7133
    model2     2 12 367.4266 387.3893 -171.7133
    Warning message:
    In anova.lme(model1, model2) :
      fitted objects with different fixed effects. REML comparisons are not meaningful.
    > 
    > # the second model collapses the main and interaction effects of the first model
    > anova(model1)
                numDF denDF   F-value p-value
    (Intercept)     1    26 1719.5059  <.0001
    Time            2    26   28.9986  <.0001
    Diet            2    13   85.3659  <.0001
    Time:Diet       4    26    1.7610  0.1671
    > anova(model2)
                numDF denDF   F-value p-value
    (Intercept)     1    24 1719.5059  <.0001
    TimeDiet        8    24   29.4716  <.0001
    > 
    > # they give the same predicted values
    > newdata <- expand.grid(Time=levels(ex$Time), Diet=levels(ex$Diet))
    > newdata$TimeDiet <- interaction(newdata$Time, newdata$Diet)
> newdata$pred1 <- predict(model1, newdata=newdata, level=0)
    > newdata$pred2 <- predict(model2, newdata=newdata, level=0)
> newdata
  Time Diet TimeDiet   pred1   pred2
1    1    1      1.1 250.625 250.625
2   22    1     22.1 261.875 261.875
3   44    1     44.1 267.250 267.250
4    1    2      1.2 453.750 453.750
5   22    2     22.2 475.000 475.000
6   44    2     44.2 488.750 488.750
7    1    3      1.3 508.750 508.750
8   22    3     22.3 518.250 518.250
9   44    3     44.3 530.000 530.000

唯一的区别是假设很容易检验。例如，在第一个模型中，很容易测试两个预测变量是否相互作用，在第二个模型中，没有明确的测试。另一方面，两个预测变量的联合效应很容易在第二个模型中测试，但不是第一个模型。其他假设是可以检验的，只是建立这些假设需要更多的工作。

其它你可能感兴趣的问题
上一篇交叉验证滥用（报告最佳超参数值的性能） 下一篇GBM、XGBoost、LightGBM、CatBoost 之间的数学差异？