Укороченная версия

factorC <- paste(factorA, factorB)
Эквивалентен ли factorA*factorBfactorA+factorB+factorC в формуле линейной модели смешанных эффектов?

Фон

Я использую функцию ancombc2() из пакета ANCOMBC. ANCOMBC — это пакет, используемый для дифференциального тестирования численности с данными о частоте видов. Подробности смотрите здесь: https://www.nature.com/articles/s41592-023-02092-7 Однако я надеюсь, что для этого вопроса нам не придется сталкиваться со спецификой этого пакета, поскольку ANCOMBC есть используя lmer() из lmerTest/lme4.

Мой вопрос касается ручного включения термина взаимодействия. ANCOMBC не поддерживает условия взаимодействия. Однако авторы пишут следующее:

Вопрос: Может ли функция ancombc или ancombc2 обрабатывать условия взаимодействия в анализе?

О: К сожалению, включение терминов взаимодействия в аргумент fix_formula команд ancombc или ancombc2 может привести к сложностям и потенциальной путанице при сравнении нескольких групп. Чтобы решить эту проблему, рекомендуется вручную создать интересующий член взаимодействия вне формулы и выполнить соответствующий анализ.

Создав вручную термин взаимодействия, вы можете гарантировать, что анализ точно уловит эффект взаимодействия между переменными. После создания термина взаимодействия вы можете включить его в аргумент 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