R语言分析不同处理组成熟动物数量增长差异的方法咨询

第一步：先把数据转成适合分析的格式 #

你现在的数据是宽格式（每个时间点单独一列），但要做时间相关的增长分析，必须转成长格式（时间作为一个独立变量列）。用tidyr包的pivot_longer就能快速搞定，代码示例：

library(tidyr)
library(dplyr)

# 假设你的原始数据框叫data_wide
data_long <- data_wide %>%
  pivot_longer(
    cols = starts_with("Mature_D"),  # 选中所有记录成熟数量的列
    names_to = "Day",                # 新列名：天数
    values_to = "Mature_Count",      # 新列名：成熟动物数量
    names_prefix = "Mature_D"        # 去掉列名里的"Mature_D"前缀
  ) %>%
  mutate(Day = as.integer(Day))      # 把Day转成数值型，方便后续建模

模型选择：绝对不要单独给每个处理建模 #

单独建模的问题太多：没法用统计检验直接比较组间斜率差异、容易犯多重比较错误、效率极低。正确的做法是把所有处理放进同一个模型，加入处理和时间的交互项，这样直接就能检验不同处理的增长速度（斜率）是否有显著差异。

推荐的模型类型 #

你的数据是计数型、左偏，还存在重复测量（同一个重复在不同时间点的观测是相关的），优先考虑：

1. 负二项混合效应模型：解决计数数据左偏、过度离散的问题，同时控制重复的随机效应（同一个重复的观测不独立），用glmmTMB包实现：

library(glmmTMB)

# 构建模型：Treatment和Day的交互项是核心，(1|Replicate)控制重复的随机效应
model <- glmmTMB(
  Mature_Count ~ Treatment * Day + (1 | Replicate),
  data = data_long,
  family = nbinom2,  # 负二项分布，适配左偏计数数据
  na.action = na.exclude  # 自动忽略NA值
)

2. 如果数据里有大量0值，再考虑零膨胀负二项模型，把family改成zinb2即可。

怎么判断组间斜率差异？ #

模型跑出来后，用summary(model)查看交互项Treatment:Day的显著性：

• 如果交互项显著：说明不同处理的增长速度（斜率）有显著差异；
• 如果交互项不显著：说明所有处理的增长速度一致，只需要关注主效应即可。

如果交互项显著，用emmeans包做事后比较，直接对比各组的斜率：

library(emmeans)

# 提取各组的斜率（Day的效应）并做两两比较
slope_comparison <- emmeans(model, ~ Day | Treatment, type = "response") %>%
  contrast(method = "pairwise")

summary(slope_comparison)

模型假设验证 #

新手不用过度纠结，但至少做这两步：

1. 检查过度离散：负二项模型本身就是用来解决过度离散的，重点看残差分布；
2. 残差可视化：用DHARMa包画残差图，判断是否符合模型假设：

library(DHARMa)

sim_res <- simulateResiduals(model)
plot(sim_res)

如果残差图明显不符合假设，再考虑换模型（比如零膨胀模型）。

绘图建议：直观展示增长趋势 #

用ggplot2画原始数据点+拟合的回归线，清晰展示各组的增长差异：

library(ggplot2)

# 先生成模型的预测值，用来画拟合线
pred_data <- expand.grid(
  Day = 4:12,
  Treatment = unique(data_long$Treatment),
  Replicate = unique(data_long$Replicate)
)
pred_data$Mature_pred <- predict(model, newdata = pred_data, type = "response")

# 绘图
ggplot(data_long, aes(x = Day, y = Mature_Count, color = Treatment)) +
  geom_point(alpha = 0.5, size = 1.5) +  # 原始数据点，透明化避免重叠
  geom_line(data = pred_data, aes(y = Mature_pred), linewidth = 1) +  # 拟合的增长曲线
  labs(x = "天数", y = "成熟动物数量", color = "处理组") +
  theme_bw() +
  theme(legend.position = "bottom")

关于NA值的处理 #

如果NA是随机缺失（比如偶然没记录到），直接用na.exclude忽略即可；如果是非随机缺失（比如某组某时间点全部缺失），先排查实验原因，再考虑是否剔除对应组或时间点。

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