R语言growthrates包批量绘制fit_easylinear拟合结果遇报错

我用R语言的growthrates包中fit_easylinear函数处理96孔板实验数据，这个函数要求输入数据无重复。我的数据里每个孔（A1、A2等）对应重复的time列，已经借助Stack Overflow用户的帮助完成批量拟合，但绘图时遇到两个问题：

• 用plot(slot(fit_list[[i]], "obs"))只能画出观测值，看不到拟合曲线；
• 用plot(fit_list[[i]])会报错：Error in seq.default(min(obs[, "time"] + lag), max(obs[, "time"]), length = 200) : 'from' must be a finite number。

现有数据与拟合代码 #

library(growthrates)
library(tidyverse)

set.seed(11)

df <- data.frame(time = seq(0, 865),
                 A1 = runif(866, 0, 1),
                 A2 = runif(866, 0, 1),
                 A3 = runif(866, 0, 1),
                 A4 = runif(866, 0, 1),
                 A5 = runif(866, 0, 1),
                 A6 = runif(866, 0, 1),
                 A7 = runif(866, 0, 1),
                 A8 = runif(866, 0, 1),
                 A9 = runif(866, 0, 1),
                 A10 = runif(866, 0, 1),
                 A11 = runif(866, 0, 1),
                 A12 = runif(866, 0, 1),
                 B1 = runif(866, 0, 1),
                 B2 = runif(866, 0, 1),
                 B3 = runif(866, 0, 1),
                 B4 = runif(866, 0, 1),
                 B5 = runif(866, 0, 1),
                 B6 = runif(866, 0, 1),
                 B7 = runif(866, 0, 1),
                 B8 = runif(866, 0, 1),
                 B9 = runif(866, 0, 1),
                 B10 = runif(866, 0, 1),
                 B11 = runif(866, 0, 1),
                 B12 = runif(866, 0, 1),
                 C1 = runif(866, 0, 1),
                 C2 = runif(866, 0, 1),
                 C3 = runif(866, 0, 1),
                 C4 = runif(866, 0, 1),
                 C5 = runif(866, 0, 1),
                 C6 = runif(866, 0, 1),
                 C7 = runif(866, 0, 1),
                 C8 = runif(866, 0, 1),
                 C9 = runif(866, 0, 1),
                 C10 = runif(866, 0, 1),
                 C11 = runif(866, 0, 1),
                 C12 = runif(866, 0, 1))


# 批量拟合easylinear模型
fit_list <- lapply(2:length(colnames(df)), function(x) fit_easylinear(df$time, df[[x]]))

尝试过的绘图代码 #

# 仅绘制观测值，无拟合曲线
for (i in seq_along(fit_list)) {
  jpeg(paste0("C:/Users/Desktop/images/", colnames(df)[i+1], ".jpg"))
  plot(slot(fit_list[[i]], "obs"))
  dev.off() 
}

# 报错的绘图代码
for (i in seq_along(fit_list)) {
  jpeg(paste0("C:/Users/Rahul/Desktop/images", colnames(df)[i+1], ".jpg"))
  plot(fit_list[[i]]) 
  dev.off()
}

问题分析与解决方案 #

问题1：仅显示观测值，无拟合曲线 #

plot(slot(fit_list[[i]], "obs"))只调用了拟合结果中的观测数据部分，自然不会显示拟合曲线。要同时展示观测值和拟合曲线，需要手动提取拟合参数，生成拟合曲线数据后叠加绘图。

问题2：plot(fit_list[[i]])报错 #

报错原因是拟合结果中lag参数为NA或非有限值，导致绘图函数计算序列时出错。这是因为fit_easylinear默认会尝试估计滞后阶段，但随机生成的测试数据没有明显的生长趋势，导致模型无法有效估计lag。

修复后的绘图代码 #

以下代码可以解决两个问题，既显示观测值又叠加拟合曲线，同时避免报错：

# 创建输出目录（如果不存在）
dir.create("C:/Users/Desktop/images/", recursive = TRUE, showWarnings = FALSE)

for (i in seq_along(fit_list)) {
  fit <- fit_list[[i]]
  well_name <- colnames(df)[i+1]
  
  # 提取观测数据
  obs_data <- slot(fit, "obs")
  # 提取拟合参数
  params <- coef(fit)
  # 生成拟合曲线的时间序列
  fit_time <- seq(min(obs_data$time), max(obs_data$time), length.out = 200)
  # 根据easylinear模型计算拟合值：y = y0 + mu*(time - lag)，当time > lag时；否则为y0
  # 处理lag为NA的情况，默认设为0
  lag_val <- ifelse(is.finite(params["lag"]), params["lag"], 0)
  fit_y <- ifelse(fit_time > lag_val, 
                  params["y0"] + params["mu"]*(fit_time - lag_val),
                  params["y0"])
  
  # 绘图并保存
  jpeg(paste0("C:/Users/Desktop/images/", well_name, ".jpg"), width = 800, height = 600)
  plot(value ~ time, data = obs_data, main = well_name, xlab = "Time", ylab = "Value", pch = 16, cex = 0.5)
  lines(fit_time, fit_y, col = "red", lwd = 2)
  legend("topleft", legend = c("Observed", "Fitted"), col = c("black", "red"), pch = c(16, NA), lty = c(NA, 1))
  dev.off()
}

额外说明 #

如果你的真实实验数据有明显生长趋势，fit_easylinear会正确估计lag参数，此时也可以直接使用plot(fit)，但建议先检查拟合结果的参数是否为有限值：

# 检查拟合参数是否有效
lapply(fit_list, function(f) all(is.finite(coef(f))))

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