使用purrr包执行含交互项的线性回归及ANCOVA分析时的模型结果异常问题求助

问题根源 #

你遇到的核心问题是数据结构不符合模型的要求：看你的df数据，每个n（比如tree1、tree2）对应的分组里，Type的每个水平（a、b、c）都只有1个观测值。当你用nest(data = -n)按tree分组后，每个子数据集里只有3个点（每个Type各1个），这时候拟合PD~SR*Type的交互模型完全没有意义——因为每个Type在单个tree里没有重复的SR-PD观测，模型根本无法估计出每个Type的独立斜率和截距，最终输出的系数是模型在秩亏（rank-deficient）情况下的默认结果，所以看起来每个Type的系数都一样。

你可以用下面的代码验证这一点：

df %>% count(n, Type)

输出会显示每个n-Type组合的观测数都是1，这就是问题所在。

解决思路与代码示例 #

根据你的研究目标，这里提供两种可行的调整方案：

方案1：按Type分组拟合独立回归（适合目标是看每个Type的PD-SR关系） #

如果你的需求是分别估计每个Type下PD与SR的线性关系，那应该按Type而不是n来嵌套数据，这样每个Type下有4个观测值（来自不同的tree），足够拟合线性模型：

library(tidyverse)
library(broom)

# 原始数据
df <- data.frame(PD=c(10,20,30,40,50,10,20,33,12,52,21,43), 
                 SR=c(5,10,20,24,6,21,59,22,1,11,12,3), 
                 n=c("tree1", "tree1", "tree1", "tree2", "tree2","tree2", "tree3", "tree3", "tree3", "tree4", "tree4","tree4"), 
                 Type=c("a", "b",'c', "a", "b",'c', "a", "b",'c', "a", "b",'c'))

# 按Type嵌套数据并拟合模型
type_models <- df %>%
  nest(data = -Type) %>%
  mutate(
    fit = map(data, ~lm(PD ~ SR, data = .x)),  # 每个Type单独拟合PD~SR
    tidied = map(fit, tidy, conf.int = TRUE)  # 提取系数及置信区间
  )

# 整理成可读的结果
tidy_results <- type_models %>%
  unnest(tidied) %>%
  select(Type, term, estimate, conf.low, conf.high, p.value)

print(tidy_results)

这样你就能得到每个Type对应的截距和斜率，以及它们的统计显著性。

方案2：用混合效应模型做ANCOVA（适合考虑tree为随机效应的情况） #

如果你的目标是做ANCOVA分析（检验Type对PD的影响，控制SR，同时检验SR与Type的交互），并且n（tree）是随机分组（比如不同的tree是重复抽样的单位），那不需要拆分数据，直接用混合效应模型更合适：

library(lme4)
library(broom.mixed)

# 拟合混合效应ANCOVA模型：固定效应为SR*Type，随机效应为tree的截距
ancova_mixed <- lmer(PD ~ SR * Type + (1|n), data = df)

# 提取结果
tidy(ancova_mixed, conf.int = TRUE)

这个模型会同时估计：

• 固定效应：SR的主效应、Type的主效应、SR与Type的交互效应
• 随机效应：不同tree之间的截距变异
  这更符合ANCOVA的统计逻辑，也能利用所有数据的信息。

关于purrr的补充说明 #

purrr本身没有问题，它只是按你指定的分组执行了模型拟合。问题出在分组方式和数据量的匹配上。如果未来你需要用purrr批量拟合模型，一定要先检查每个分组内的样本量是否足够支撑你要拟合的模型复杂度——比如带交互项的模型需要每个分组内有足够的观测来估计交互效应。

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