如何从R语言lm()模型中推导正确的y关于x的计算公式？

问题描述 #

我用R语言的lm()函数拟合曲线，想从模型里导出计算y的方程应用到更大数据集，但推导的方程和模型预测曲线不符。之前做线性/多元回归能正常提取方程，但这次模型包含常数K、a、p，推导后曲线不对。

数据与代码如下：

# Data
x <- c(18.212128, 23.453468, 30.275368, 18.721986, 22.831441, 12.624084, 13.572420,  9.564936,  4.282807, 13.042168, 12.032649, 13.082956, 13.144139, 11.247467, 10.074794, 14.123067,  8.830741, 10.380709, 10.054400,  6.444605,  3.936104,  2.590079,  2.834810, 16.111513, 12.399747, 10.186963, 15.397712,  3.324274, 13.123745,  8.698177, 27.960613,  5.455481,  6.373225,  3.263091,  6.842294,  7.013732, 16.373396, 34.661337, 35.608483)

y <- c(1.5740,  28.1570, 30.9710, 3.5450, 5.7480, 0.4950,  2.9810,   0.8400,   1.7500,  0.7165,   0.9920, 0.8065, 1.0680, 0.5640, 0.1430,  1.0605, 0.2745, 0.7750,  0.6750,  0.7800,   0.4050,   0.4035,   0.4210,   7.6080,   5.5820, 2.3250, 5.7760,  0.3125,   6.9100,  3.0575,  40.0450,   3.5900,   3.8050,   0.6435,   1.2600,   2.8100,  12.4200, 176.400, 128.8000)


# set variables
K = 1
a = -4
p= .34

# Create the model
model<- lm(K*(y+a)^p ~ x)

# Create predicted values 
model.predict <- exp(predict(model,list(Level=x))) 

# Print the Summary
summary(model)

# Plot 
plot_ly(x = x, y=y)%>%
  add_markers()%>%
  add_lines(x = x, y=model.predict)%>%
#add_lines(x = x.dummy, y=y.dummy)%>%
  layout(
    xaxis = list(title = "Water Level (cm)"),
    yaxis = list(title = "Discharge")
  )

从模型摘要得到系数：(Intercept)=-1.2304，x=0.1712，我推导的方程是y = ((-1.23036 + 0.17123 * x)/K)^(1/p) - a，但生成的曲线和model.predict不一致，请问怎么从模型输出得到正确的y关于x的计算方程？

问题分析与解决 #

1. 核心错误：多余的exp()转换 #

你的代码中用exp(predict(model,list(Level=x)))生成预测值，但你的模型没有做对数变换，完全不需要用exp()反变换，这是导致曲线不符的根本原因。

2. 正确的方程推导 #

你的模型定义是：

lm(K*(y+a)^p ~ x)

对应的线性模型表达式为：
$$K \times (y + a)^p = \beta_0 + \beta_1 \times x$$
其中$\beta_0$是截距项，$\beta_1$是x的系数。

整理得到y的表达式：
$$y = \left( \frac{\beta_0 + \beta_1 \times x}{K} \right)^{\frac{1}{p}} - a$$
你的推导本身是正确的，只是预测值生成时的错误操作导致结果偏差。

3. 修正后的代码 #

正确生成预测值

# 直接获取线性模型的预测值，无需exp转换
linear_pred <- predict(model, list(x = x))
# 按推导方程转换为y的预测值
model.predict <- (linear_pred / K)^(1/p) - a

完整修正代码

# Data
x <- c(18.212128, 23.453468, 30.275368, 18.721986, 22.831441, 12.624084, 13.572420,  9.564936,  4.282807, 13.042168, 12.032649, 13.082956, 13.144139, 11.247467, 10.074794, 14.123067,  8.830741, 10.380709, 10.054400,  6.444605,  3.936104,  2.590079,  2.834810, 16.111513, 12.399747, 10.186963, 15.397712,  3.324274, 13.123745,  8.698177, 27.960613,  5.455481,  6.373225,  3.263091,  6.842294,  7.013732, 16.373396, 34.661337, 35.608483)

y <- c(1.5740,  28.1570, 30.9710, 3.5450, 5.7480, 0.4950,  2.9810,   0.8400,   1.7500,  0.7165,   0.9920, 0.8065, 1.0680, 0.5640, 0.1430,  1.0605, 0.2745, 0.7750,  0.6750,  0.7800,   0.4050,   0.4035,   0.4210,   7.6080,   5.5820, 2.3250, 5.7760,  0.3125,   6.9100,  3.0575,  40.0450,   3.5900,   3.8050,   0.6435,   1.2600,   2.8100,  12.4200, 176.400, 128.8000)

# set variables
K = 1
a = -4
p= .34

# Create the model
model<- lm(K*(y+a)^p ~ x)

# 正确生成预测值
linear_pred <- predict(model, list(x = x))
model.predict <- (linear_pred / K)^(1/p) - a

# Print the Summary
summary(model)

# Plot 
plot_ly(x = x, y=y)%>%
  add_markers()%>%
  add_lines(x = x, y=model.predict)%>%
  layout(
    xaxis = list(title = "Water Level (cm)"),
    yaxis = list(title = "Discharge")
  )

4. 验证说明 #

修正后，用推导方程计算出的model.predict和模型的预测逻辑完全一致，绘制的曲线会与数据趋势匹配。

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