基于nlme的带区组因素重复测量ANOVA代码正确性验证

整体分析逻辑符合实验设计要求：明确嵌套层级、选用nlme处理重复测量自相关、交互效应显著后开展事后多重比较的思路是正确的，仅存在两处关键笔误需要修正。

需修正的问题 #

• 随机效应嵌套顺序写反：你的实验设计为田块(plot)嵌套在区组(block)下，原代码~1|plot/block的嵌套逻辑完全颠倒，区组是更高层级的分组因素
• 事后检验显著性水平笔误：原代码alpha=0.5不符合常规统计标准，应为0.05

实验设计说明 #

• 区组：北田、南田，北田2个重复、南田1个重复，每个重复为2英亩田块，定期测定土壤硝酸盐含量
• 处理：reference、treat_1、treat_2
• 时间：3、4、5、6月
• 响应变量：硝酸盐含量no3

依赖包加载 #

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)) 

修正后模型拟合代码 #

lme_fitno3.block <- lme(fixed =no3 ~ treatment * month ,  
                    random = ~1|block/plot, # 修正嵌套顺序，区组下嵌套田块
                    method='REML',
                    corr = corAR1( form= ~1|block/plot), # 同步修正相关结构的嵌套层级
                    data = no3.df)
summary(lme_fitno3.block)
Anova(lme_fitno3.block, type="III")

修正后事后检验代码 #

marginal = emmeans(lme_fitno3.block, ~ treatment:month)
plot(marginal, comparisons = TRUE)
emminteraction = emmeans(lme_fitno3.block, 
                         pairwise ~ treatment:month,
                         adjust="bonferroni",
                         alpha=0.05) # 修正显著性水平笔误
emminteraction$contrasts
multcomp::cld(marginal,
              Letters = letters,
              adjust="bonferroni")

其他注意事项 #

你后续计划通过AIC选优协方差结构的思路可行，注意保持固定效应完全一致的前提下对比不同协方差结构的AIC即可，REML拟合的模型仅可用于比较协方差结构，不能用于比较不同固定效应的设定，刚好符合你的需求。

内容的提问来源于stack exchange，提问作者Bill Perry