lmerTest::lmer：手动创建交互项

问题描述 投票：0回答：1

精简版

factorC <- paste(factorA, factorB)

factorA*factorB
相当于线性混合效应模型公式中的 
factorA+factorB+factorC
吗？

背景

我使用包

ancombc2()
中的函数
ANCOMBC
。 
ANCOMBC
是一个用于使用物种频率数据进行差异丰度测试的软件包。有关详细信息，请参阅此处：https://www.nature.com/articles/s41592-023-02092-7但是，我希望我们不必针对这个问题处理此包的特殊性，因为
ANCOMBC
正在使用 
lmer()
/ 
lmerTest
中的 
lme4
。

我的问题是关于“手动”包含交互术语。 ANCOMBC不支持交互术语。然而，作者写道：



问：ancombc 或 ancombc2 函数可以处理分析中的交互项吗？


A：不幸的是，在 ancombc 或 ancombc2 的 fix_formula 参数中包含交互项可能会导致多组比较的复杂性和潜在的混乱。为了解决这个问题，建议在公式之外手动创建感兴趣的交互项并进行相应的分析。


通过手动创建交互项，可以确保分析准确捕捉变量之间的交互作用。创建交互项后，您可以将其包含在 fix_formula 参数或分析的任何其他相关部分中，具体取决于您的具体研究问题和设计。


来自：
https://github.com/FrederickHuangLin/ANCOMBC



我的解释和可重现的例子 
lmer()
我使用 
lmer()
 中的 
lmeTest
 进行测试。我将使用以下示例数据展示我如何解释作者的解释：
在一项实验中，对来自两个治疗组和两个年龄组的样本进行了观察。此外，该实验在两年内重复进行。因此，实验包括以下因素：




“年份”（随机），级别为“1”和“2”
    

  • “治疗”（固定），级别为“对照”和“治疗”
    • 

  • “年龄”（固定），级别为“旧”和“新”
    • 

  • 

    • 这是一个虚拟数据集：


# create example data
set.seed(123)
dat_year1 <- data.frame(year = factor(1),
                        treatment = c(rep("control", 20),
                                      rep("treated", 20)),
                        age = c(rep("old", 10),
                                rep("new", 10),
                                rep("old", 10),
                                rep("new", 10)),
                        value = c(rnorm(n = 10, mean = 10, sd = 2),
                                  rnorm(n = 10, mean = 10, sd = 2),
                                  rnorm(n = 10, mean = 15, sd = 2),
                                  rnorm(n = 10, mean = 16, sd = 2)))

set.seed(321)
dat_year2 <- data.frame(year = factor(2),
                        treatment = c(rep("control", 20),
                                      rep("treated", 20)),
                        age = c(rep("old", 10),
                                rep("new", 10),
                                rep("old", 10),
                                rep("new", 10)),
                        value = c(rnorm(n = 10, mean = 10, sd = 2),
                                  rnorm(n = 10, mean = 10, sd = 2),
                                  rnorm(n = 10, mean = 15, sd = 2),
                                  rnorm(n = 10, mean = 18, sd = 2)))

dat <- rbind(dat_year1, dat_year2)



及其可视化：



我有兴趣测试“治疗”和“年龄”之间的相互作用。对于 
lmer()
 我会这样做：
# run lmer()
test <- lmerTest::lmer(value~treatment*age+(1|year), data = dat)
summary(test)



返回：


Linear mixed model fit by REML. t-tests use
  Satterthwaite's method [lmerModLmerTest]
Formula: value ~ treatment * age + (1 | year)
   Data: dat

REML criterion at convergence: 322.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.41222 -0.69075  0.07493  0.44019  2.22631 

Random effects:
 Groups   Name        Variance Std.Dev.
 year     (Intercept) 0.05503  0.2346  
 Residual             3.44757  1.8568  
Number of obs: 80, groups:  year, 2

Fixed effects:
                        Estimate Std. Error      df t value
(Intercept)              10.6492     0.4471  7.6713  23.819
treatmenttreated          6.6717     0.5872 75.0000  11.363
ageold                   -0.4259     0.5872 75.0000  -0.725
treatmenttreated:ageold  -2.6641     0.8304 75.0000  -3.208
                        Pr(>|t|)    
(Intercept)             1.81e-08 ***
treatmenttreated         < 2e-16 ***
ageold                   0.47044    
treatmenttreated:ageold  0.00196 ** 
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) trtmnt ageold
tretmnttrtd -0.657              
ageold      -0.657  0.500       
trtmnttrtd:  0.464 -0.707 -0.707



我将创建交互项的建议解释为应手动创建一个新因素，该因素结合了要测试交互的因素的级别：


dat$treatment_age <- paste(dat$treatment, "_", dat$age, sep = "")



现在我将其作为另一个固定因子包含在模型公式中：


test_manual <- lmerTest::lmer(value~treatment+age+treatment_age+(1|year),
                              data = dat)
summary(test_manual)



返回：


Linear mixed model fit by REML. t-tests use
  Satterthwaite's method [lmerModLmerTest]
Formula: 
value ~ treatment + age + treatment_age + (1 | year)
   Data: dat

REML criterion at convergence: 322.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.41222 -0.69075  0.07493  0.44019  2.22631 

Random effects:
 Groups   Name        Variance Std.Dev.
 year     (Intercept) 0.05503  0.2346  
 Residual             3.44757  1.8568  
Number of obs: 80, groups:  year, 2

Fixed effects:
                         Estimate Std. Error      df
(Intercept)               10.6492     0.4471  7.6713
treatmenttreated           6.6717     0.5872 75.0000
ageold                    -3.0900     0.5872 75.0000
treatment_agecontrol_old   2.6641     0.8304 75.0000
                         t value Pr(>|t|)    
(Intercept)               23.819 1.81e-08 ***
treatmenttreated          11.363  < 2e-16 ***
ageold                    -5.263 1.30e-06 ***
treatment_agecontrol_old   3.208  0.00196 ** 
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) trtmnt ageold
tretmnttrtd -0.657              
ageold       0.000 -0.500       
trtmnt_gcn_ -0.464  0.707 -0.707
fit warnings:
fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients



对于“治疗”，所有统计数据都是相同的，对于“年龄”，它们有很大不同，对于“治疗：年龄”，估计值和 t 值现在为负值。此外，还有一个拟合警告“固定效应模型矩阵存在秩不足，因此删除 2 列/系数”。


问题


如何改进以手动正确包含交互项？


与
ancombc2
一起使用

当我尝试将 
ancombc2()
 与手动包含的交互项一起使用时，它不会执行但返回错误：
library(ANCOMBC)
library(phyloseq)

ancombc2(data = phyloseq_object,
         tax_level = NULL,
         fix_formula = "treatment+age+treatment_age",
         verbose = TRUE)



返回：


Obtaining initial estimates ...

Error: Estimation failed for the following covariates:
treated_old, treated_new
Please ensure that these covariates do not have missing values and check for multicollinearity before re-estimating the model


    



    
    r
    
    lme4
    
    interaction
    
    lmertest
    
    
1个回答

    

      
      
      

        

          

            

              
0
投票

              
            

          

        

        
factorC
本身相当于
factorA*factorB
 - 您不需要包含
factorA+factorB
。这就是为什么您在较长版本中收到排名不足警告的原因，
treatment_age
中的两个级别已被
treatment
和
age
主要效果覆盖。
    



      

      
    

  

  

    
最新问题


  





  

    © www.soinside.com 2019 - 2024. All rights reserved.