使用 nlme 的阻塞因子重复测量方差分析

Question

我想确保这是正确的，尽管我认为它与 Whosebug 上的其他版本相似但不完全相同。

经验设计：

街区 - 北田和南田
治疗 - 参考，treat_1，treat_2
时间为月 - 3、4、5、6
响应变量是硝酸盐 - no3

北场有2个重复，南场有1个重复。重复是 2 英亩的田地，我们在土壤中随时间测量硝酸盐，因为它对处理有反应。

包是：

library(tidyverse) 
library(car)
library(multcompView)
library(nlme)
library(emmeans)

下面是一个简化的数据框。

no3.df <- structure(list(month = c(3, 3, 3, 4, 5, 5, 5, 5, 6, 3, 3, 3, 
                4, 5, 5, 5, 5, 6, 3, 4, 5, 5, 5, 5, 6, 3, 5, 5, 5, 5, 6, 3, 3, 
                3, 4, 6, 3, 3, 3, 4, 5, 5, 5, 3, 3, 4, 5, 5, 5, 5, 6, 3, 3, 3, 
                4, 5, 5, 5, 5, 6, 3, 3, 3, 4, 5, 5, 5, 5, 6), 
                block = c("north", "north", "north", "north", "north", "north", "north", "north", 
                        "north", "north", "north", "north", "north", "north", "north", 
                        "north", "north", "north", "south", "south", "south", "south", 
                        "south", "south", "south", "north", "north", "north", "north", 
                        "north", "north", "north", "north", "north", "north", "north", 
                        "south", "south", "south", "south", "south", "south", "south", 
                        "north", "north", "north", "north", "north", "north", "north", 
                        "north", "north", "north", "north", "north", "north", "north", 
                        "north", "north", "north", "south", "south", "south", "south", 
                        "south", "south", "south", "south", "south"), 
                plot = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 8, 8, 8, 8, 8, 
                        8, 8, 3, 3, 3, 3, 3, 3, 5, 5, 5, 5, 5, 9, 9, 9, 9, 9, 9, 9, 2, 
                        2, 2, 2, 2, 2, 2, 2, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 
                        7, 7, 7, 7), 
                treatment = c("treat_1", "treat_1", "treat_1", "treat_1", 
                       "treat_1", "treat_1", "treat_1", "treat_1", "treat_1", "treat_1", 
                       "treat_1", "treat_1", "treat_1", "treat_1", "treat_1", "treat_1", 
                       "treat_1", "treat_1", "treat_1", "treat_1", "treat_1", "treat_1", 
                       "treat_1", "treat_1", "treat_1", "treat_2", "treat_2", "treat_2", 
                       "treat_2", "treat_2", "treat_2", "treat_2", "treat_2", "treat_2", 
                       "treat_2", "treat_2", "treat_2", "treat_2", "treat_2", "treat_2", 
                       "treat_2", "treat_2", "treat_2", "reference", "reference", "reference", 
                       "reference", "reference", "reference", "reference", "reference", 
                       "reference", "reference", "reference", "reference", "reference", 
                       "reference", "reference", "reference", "reference", "reference", 
                       "reference", "reference", "reference", "reference", "reference", 
                       "reference", "reference", "reference"), 
                no3 = c(36.8, 20.4925, 21.03333333, 16.33, 7.723, 1.566333333, 0.533333333, 0.189, 0.31, 
                     25.8, 16.13333333, 24.86666667, 3.979, 1.814, 0.34635, 0.244666667, 
                     0.247333333, 0.97675, 14.305, 11.91, 12.4, 6.79, 7.26825, 8.4615, 
                     3.43575, 22.225, 0.3243, 0.1376, 0.6244, 0.962233333, 1.36675, 
                     8.27, 14.96, 19.62, 44.7, 9.197, 15.6, 13.85, 17.76, 14.84, 17.8, 
                     23.06, 12.19333333, 19.06, 22.675, 27.47, 18.295, 16.5425, 18.7375, 
                     22.25333333, 24.63125, 21.75, 23.73333333, 13.09, 20.54, 17.1, 
                     10.58666667, 17.5565, 20.5, 25.575, 19.8, 15.76666667, 18.25333333, 
                     15.93, 11.89, 10.791, 22.65, 22.025, 23.93333333)), 
           row.names = c(NA, -69L), class = c("tbl_df", "tbl", "data.frame"))

读入数据并生成因子

no3.df <- no3.df %>% 
  mutate( 
         treatment = as.factor(treatment),
         plot=as.factor(plot),
         month=as.factor(month))

我正在使用 nlme 来指定 covariance/variance 结构。最终我会用其他协方差和方差结构来尝试这个，并查看 AIC 看看什么是最好的，但现在我认为这种方法可能最适合作为 AR1 矩阵。

lme_fitno3.block <- lme(fixed =no3 ~ treatment * month ,  
                    random = ~1|plot/block, 
                    method='REML',
                    corr = corAR1( form= ~1|plot/block),
                    data = no3.df)
summary(lme_fitno3.block)
Anova(lme_fitno3.block, type="III")

模型结果为

Analysis of Deviance Table (Type III tests)

Response: no3
                  Chisq Df Pr(>Chisq)    
(Intercept)     50.8817  1  9.810e-13 ***
treatment        1.9561  2      0.376    
month            3.4219  3      0.331    
treatment:month 29.7859  6  4.317e-05 ***

我从中得出治疗和月份之间存在显着的相互作用，然后进行后续测试。

marginal = emmeans(lme_fitno3.block, 
                   ~ treatment:month)

plot(marginal, comparisons = TRUE)

emminteraction = emmeans(lme_fitno3.block, 
                         pairwise ~ treatment:month,
                         adjust="bonferroni",
                         alpha=0.5)
emminteraction$contrasts
multcomp::cld(marginal,
              Letters = letters,
              adjust="bonferroni")

我不会 post 结果，因为它们很广泛。

Answer 1

对于那些可能感兴趣的人，我确实在 website for repeated measures and random factors

找到了 nlme 和 lme4 的一个很好的介绍和示例 R 代码

使用 nlme 的阻塞因子重复测量方差分析

Repeated Measures ANOVA with blocking factor using nlme

r

anova

mixed-models

nlme

emmeans