R语言DRC包拟合多分组DRM剂量反应模型出现下标越界错误如何解决

错误原因 #

drm() 函数的curveid参数（公式后第一个位置参数）仅支持传入单个分组变量，你直接在该位置传入Trt,Timing两个独立变量，函数无法识别多分组结构，因此抛出subscript out of bounds错误。

解决方法 #

将两个分组变量合并为一个唯一标识每个曲线组的组合变量后再传入curveid参数即可。

修正后完整代码 #

1. 数据预处理阶段新增组合分组变量

library(pacman) 
library(dplyr) 
p_load(MASS, drc, ggplot2, reshape2) 

t1= read.csv('TimingBB.csv',sep = ",",header=TRUE,strip.white=TRUE) 
t2 <- melt(t1,id.vars=c('Dose', 'Rep', 'Timing'), variable.name = "Trt", value.name = 'Germination') 

# 新增代码：创建Trt和Timing的组合分组变量
t2$group <- interaction(t2$Trt, t2$Timing, sep = "_")

2. 修正模型拟合代码

# 修正后的全模型调用
modl_Full <- drm(Germination ~ Dose, curveid = group, data = t2, fct = LL.3())

3. 多分组模型可视化代码（可选）
   如果需要按原需求分面+按Timing配色出图，可通过预测值生成拟合线：

# 生成预测网格数据
newdata <- expand.grid(Dose = seq(min(t2$Dose), max(t2$Dose), length.out = 100),
                       group = unique(t2$group))
newdata$pred <- predict(modl_Full, newdata = newdata)
# 拆分组合变量得到原分组字段
newdata <- cbind(newdata, colsplit(newdata$group, "_", c("Trt", "Timing")))

# 出图
ggplot(t2, aes(x=Dose, y=Germination, color=Timing))+
  geom_point(alpha = 0.7)+
  geom_line(data = newdata, aes(y = pred), linewidth = 1)+
  facet_wrap(~Trt,nrow = 3)

额外提示 #

如果需要对不同组的模型参数做约束检验（比如检验不同Timing的曲线斜率是否一致），可以通过pmodels参数指定参数结构，LL.3()的三个参数顺序为斜率b、下限c、上限d，pmodels传入的列表元素分别对应三个参数的约束规则：

# 示例：约束所有组的斜率b相同，下限和上限随分组变化
modl_constrained <- drm(Germination ~ Dose, curveid = group, data = t2, fct = LL.3(),
                  pmodels = list(~1, ~group, ~group))

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