两变量曲线关系的拟合（代码+可视化+结果解读）

医学和生信笔记

发布于 2026-05-22 21:01:29

560

文章被收录于专栏：医学和生信笔记医学和生信笔记

曲线拟合

有些情况两个变量的关系并不是直线形式的，有可能是曲线形式的，此时可以通过曲线拟合来刻画两变量之间的关系。主要方法就是对变量做转换，比如log、平方、多项式、样条等。

孙振球《医学统计学》第4版例9-11。

以不同剂量的标准促肾上腺皮质激素释放因子CRF（nmol/L）刺激离体培养的大鼠腺垂体细胞，监测其垂体合成分泌肾上腺皮质激素ACTH的量（pmol/L）。根据测得的5对数据建立ACTH-CRF工作曲线。

data9_11 <- foreign::read.spss("datasets/例09-11.sav", to.data.frame = T)
## re-encoding from CP936
str(data9_11)
## 'data.frame':    5 obs. of  3 variables:
##  $ number: num  1 2 3 4 5
##  $ x     : num  0.005 0.05 0.5 5 25
##  $ y     : num  34.1 58 94.5 128.5 170
##  - attr(*, "variable.labels")= Named chr [1:3] "编号" " CRF浓度" "ACTH的合成量"
##   ..- attr(*, "names")= chr [1:3] "number" "x" "y"
##  - attr(*, "codepage")= int 936

先画图查看趋势：

library(ggplot2)

ggplot(data9_11, aes(x,y))+
  geom_point(size=4)

可以发现这个趋势非常像高中学过的对数函数的图像，所以我们选择对自变量X做对数转换，再画图看一看：

ggplot(data9_11, aes(log10(x),y))+
  geom_point(size=4)

果然就基本上呈直线趋势了，所以我们选择对数转换后的X建立直线回归方程：

f9_11 <- lm(y ~ log10(x), data = data9_11)
summary(f9_11)
## 
## Call:
## lm(formula = y ~ log10(x), data = data9_11)
## 
## Residuals:
##      1      2      3      4      5 
##  7.152 -5.083 -4.698 -6.804  9.433 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  110.060      4.095   26.88 0.000113 ***
## log10(x)      36.115      2.968   12.17 0.001195 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.838 on 3 degrees of freedom
## Multiple R-squared:  0.9801, Adjusted R-squared:  0.9735 
## F-statistic: 148.1 on 1 and 3 DF,  p-value: 0.001195

把决定系数R2开方，得到的R值称为相关指数（correlation-index），其值在0到1之间，此数值离1越近，表示两变量间关系越密切。

• 如果两变量X与Y为直线关系，则相关指数在数值上等于X与Y的积差相关系数r的绝对值；
• 如果两变量X与Y为变换X之后可直线化的曲线关系，则相关指数在数值上等于变换后的X与Y的积差相关系数r的绝对值；
• 如果两变量X与Y为变换Y之后可直线化的曲线关系，那么相关指数须通过公式（9-42）（孙振球《医学统计学》第4版P148）定义的决定系数R开方得到，或者等于Y与Y的积差相关系数的绝对值，而不能等于X与Y或者X与Y的积差相关系数r的绝对值。

更一般地说，不论何种情况，Y与Y的相关系数绝对值即相关指数R，可以反映两变量曲线关系的密切程度。前面我们提到，积差相关系数r为0不一定表示两变量没有关系，而只是说没有直线关系，如果是曲线关系，可以用相关指数R来描述这种关系的密切程度。

该例中相关系数R2即为：Multiple R-squared: 0.9801，相关指数为：

sqrt(0.9801)
## [1] 0.99

或者根据定义直接计算积差相关系数：

cor(log10(data9_11$x),data9_11$y)
## [1] 0.9900221

下面是另外一个示例。

孙振球《医学统计学》第4版例9-12。

一位医院管理人员想建立一个回归模型，对重伤患者出院后的长期恢复情况进行预测。自变量为患者住院天数（X），应变量为患者出院后长期恢复的预后指数（Y），指数取值越大表示预后结局越好。

data9_12 <- foreign::read.spss("datasets/例09-12.sav", to.data.frame = T,
                               reencode = "gbk")
## re-encoding from gbk
str(data9_12)
## 'data.frame':    15 obs. of  2 variables:
##  $ x: num  2 5 7 10 14 19 26 31 34 38 ...
##  $ y: num  54 50 45 37 35 25 20 16 18 13 ...
##  - attr(*, "variable.labels")= Named chr [1:2] "住院天数" "预后指数"
##   ..- attr(*, "names")= chr [1:2] "x" "y"

先画图看趋势：

ggplot(data9_12, aes(x,y))+
  geom_point(size=4)

这个图有点像指数函数的图像，我们可以尝试对因变量Y做对数转换，再画图看看：

ggplot(data9_12, aes(x,log(y)))+
  geom_point(size=4)

果然就基本上呈直线趋势了，所以我们选择对数转换后的Y建立直线回归方程：

f9_12 <- lm(log(y) ~ x, data = data9_12)
summary(f9_12)
## 
## Call:
## lm(formula = log(y) ~ x, data = data9_12)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.37241 -0.07073  0.02777  0.05982  0.33539 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.037159   0.084103   48.00 5.08e-16 ***
## x           -0.037974   0.002284  -16.62 3.86e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1794 on 13 degrees of freedom
## Multiple R-squared:  0.9551, Adjusted R-squared:  0.9516 
## F-statistic: 276.4 on 1 and 13 DF,  p-value: 3.858e-10

此时得到的Multiple R-squared:0.9551是x和log10(y)的决定系数，并不是原始的Y和x的决定系数，原始的Y和x的决定系数和相关指数计算如下：

# 计算预测值（在log10尺度）
pred_logy <- predict(f9_12)
# 反变换回原始Y尺度
pred_y <- exp(pred_logy)
# 计算相关指数R = |cor(Y, Ŷ)|
R_index <- abs(cor(data9_12$y, pred_y))
R_index # 相关指数，和书中相差0.001，可忽略
## [1] 0.9933893
R_squared <- R_index^2
R_squared # 决定系数
## [1] 0.9868222

回归方程为：

Ŷ′ = 4.037 − 0.038X

因为此时使用的是log10(Y)，因此：

Ŷ1 = e(4.037 − 0.038X) = 56.66e−0.038X

计算原始尺度上的残差平方和（RSS）：

# 计算原始尺度的Y的预测值
y_pred_original <- exp(4.037-0.038 * data9_12$x)
# 计算原始尺度的 RSS
sum((data9_12$y - y_pred_original)^2)
## [1] 56.14713

以上这种曲线直线化的方法计算出的残差平方和可能并不是最小的，RSS越小说明模型拟合的越准确，因此还可以使用非线性最小二乘法拟合模型，并计算一下此时的RSS。

R语言中通过nls()拟合非线性最小二乘法模型，模型需要一组初始值。通过上面的线性模型我们已得到初步的模型公式为：

Ŷ1 = e(4.037 − 0.038X)

拟合非线性最小二乘模型时需要提供的初始值就是上面公式中的常数项，这个数值可能需要多次尝试才能得到最合适的，如果你不嫌烦，你可以通过机器学习中超参数调优的方式挨个试，用RSS最小的那个，或者你也可以直接使用线性模型的结果，本例中也就是直接使用4.037和-0.038。

拟合模型：

# 最主要的是提供a和b的值，我们就用线性模型的尝试即可
fit_nls <- nls(y ~ exp(a - b * x), # 直接写公式即可
               data = data9_12,
               start = list(a = 4.037, b = 0.038) # 提供初始值
               )
summary(fit_nls)
## 
## Formula: y ~ exp(a - b * x)
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a 4.070846   0.025119  162.06  < 2e-16 ***
## b 0.039586   0.001711   23.13 6.01e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.951 on 13 degrees of freedom
## 
## Number of iterations to convergence: 3 
## Achieved convergence tolerance: 8.78e-06

计算RSS：

sum(residuals(fit_nls)^2)
## [1] 49.4593

非线性最小二乘法的RSS确实更小。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2026-05-19，如有侵权请联系 cloudcommunity@tencent.com 删除

data

本文分享自医学和生信笔记微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

登录后参与评论

0 条评论

热度