如何针对特定处理及天数子集拟合线性回归并找-2MPa交点

需求与问题 #

• 目标：针对biochar和no biochar两种处理，筛选干旱条件下天数20至64的子集数据，拟合days与水势（mean_wp）的线性回归，绘制回归线并查看其与-2MPa的交点
• 问题：现有代码无法正确筛选子集、完成拟合与绘图，仅能输出单组回归系数

原代码与数据 #

原代码 #

##fitting a linear line on graph to find TLP####
##subset to biochar and no biochar in drought##
dat_line<-subset(dat, days > 20 & days < 64)
line_BD<-subset(dat_line,drought=="drought" & treatment=="biochar")
line_NBD<-subset(dat_line, drought=="drought" & treatment=="no biochar")

lm(days~wp_Mpa, data=line_BD)
abline(a=28.243, b=-1.919)

lm(days~wp_Mpa, data=line_NBD)
abline(a=20.83, b=-2.80)

单处理数据集（biochar干旱组） #

treatment drought days  N    mean_wp        sd        se
1    biochar drought    0 12 -0.5784129 0.6308831 0.1821203
2    biochar drought    1  8 -1.4366956 0.9335251 0.3300510
3    biochar drought    3  7 -2.0759137 1.0529204 0.3979665
4    biochar drought    8  6 -1.7999920 0.5245320 0.2141393
5    biochar drought   11  5 -1.9390823 0.9582212 0.4285295
6    biochar drought   17 10 -1.3570956 0.6604088 0.2088396
7    biochar drought   21  6 -1.3617151 1.8131385 0.7402107
8    biochar drought   23  3 -1.6917443 0.4362822 0.2518876
9    biochar drought   29  3 -1.8935309 0.7635187 0.4408178
10   biochar drought   36  6 -3.0623077 1.2009647 0.4902918
11   biochar drought   43 10 -3.4186978 2.7052291 0.8554686
12   biochar drought   64 15 -4.4986011 1.8798342 0.4853711

原回归输出（biochar组） #

Call:
lm(formula = days ~ wp_Mpa, data = line_BD)

Coefficients:
(Intercept)       wp_Mpa  
     28.243       -1.919  

问题分析与解决方案 #

核心问题 #

1. 子集筛选错误：原代码days > 20 & days < 64排除了days=20和days=64，不符合需求；
2. 变量名不匹配：数据集中水势列是mean_wp，但回归公式用了wp_Mpa，导致模型调用错误；
3. 未保存模型对象：直接调用lm()未保存结果，无法复用系数或计算交点；
4. 绘图逻辑缺失：仅用abline()但未先绘制基础散点图，导致无绘图输出。

完整修正代码 #

# 1. 正确筛选子集：包含20-64天的干旱处理数据
dat_line <- subset(dat, drought == "drought" & days >= 20 & days <= 64)
line_BD <- subset(dat_line, treatment == "biochar")
line_NBD <- subset(dat_line, treatment == "no biochar")

# 2. 拟合线性回归并保存模型
model_BD <- lm(days ~ mean_wp, data = line_BD)
model_NBD <- lm(days ~ mean_wp, data = line_NBD)

# 3. 绘制散点图与回归线
# 先设置绘图窗口，添加基础散点
plot(days ~ mean_wp, data = line_BD, col = "red", pch = 16, 
     xlab = "水势 (MPa)", ylab = "干旱天数", xlim = c(-5, 0), ylim = c(20, 70))
points(days ~ mean_wp, data = line_NBD, col = "blue", pch = 17)

# 添加回归线
abline(model_BD, col = "red", lwd = 2)
abline(model_NBD, col = "blue", lwd = 2)

# 4. 添加-2MPa水平线并计算交点
abline(v = -2, col = "gray", lty = 2)

# 计算biochar组交点天数
bd_intercept_day <- coef(model_BD)[1] + coef(model_BD)[2] * (-2)
# 计算no biochar组交点天数（假设模型已拟合）
nbd_intercept_day <- coef(model_NBD)[1] + coef(model_NBD)[2] * (-2)

# 标注交点
points(x = -2, y = bd_intercept_day, col = "red", pch = 18, cex = 1.5)
points(x = -2, y = nbd_intercept_day, col = "blue", pch = 18, cex = 1.5)
text(x = -2.2, y = bd_intercept_day, labels = round(bd_intercept_day, 1), col = "red")
text(x = -2.2, y = nbd_intercept_day, labels = round(nbd_intercept_day, 1), col = "blue")

# 添加图例
legend("bottomleft", legend = c("Biochar", "No Biochar"), 
       col = c("red", "blue"), pch = c(16,17), lwd = 2)

代码说明 #

• 子集筛选：用days >=20 & days <=64确保包含目标天数范围；
• 变量匹配：回归公式使用数据集中实际存在的mean_wp列；
• 模型保存：将lm()结果保存为对象，方便后续提取系数、计算交点；
• 绘图逻辑：先绘制散点图，再添加回归线、水平线与交点标注，确保可视化完整；
• 交点计算：利用回归方程days = 截距 + 斜率*mean_wp，代入mean_wp=-2计算对应的天数。

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