使用 R 解析非线性模型与洞察的入门指南

本教程由社区成员贡献，已由 DataCamp 为确保清晰与准确进行编辑。

想要分享您的专业见解？我们非常期待！欢迎通过我们的社区投稿表单提交文章或想法。

非线性模型如今不再像过去那样流行。现在有许多书籍要么跳过非线性模型的讨论，要么仅作有限介绍。造成这一现象的原因并不少。首先，通过在线性回归中使用高阶多项式项，或应用基函数（也就是对解释变量进行对数或平方根变换），就能拟合非线性曲线。其次，出于可解释性与表述的便利，很多人会倾向于选择更简单的线性模型。

此外，要应用非线性模型，通常需要先对数据进行可视化并把握其走势，而这在多变量场景下并不容易。最后，非线性模型更适用于机制性问题，即现象主要由物理规律或确定性过程驱动；而在大数据中，用户行为多由机器记录，情况并非如此。

然而，尽管存在这些合理考量，非线性建模仍是一门值得探索并与其线性对应方法加以比较的精妙科学。

非线性回归模型

非线性模型的估计技术本质上是迭代式的。此外，非线性公式的应用方式与线性公式有根本不同。下面这三种不同方程/公式都被认为是线性的，因为变量 y 与 B₀、B₁ 和 B₂ 之间的关系是直线。

而使公式呈现非线性的，是如下这种形式：

是否使用 B₀ 或 B₁ 并非关键。然而，后面两个方程中，变量 y 与 B₀ 的关系不是直线。

不幸的是，这会让计算更加困难，对估计问题的数值解也无法得到保证。因此，非线性模型自带“风险提示”：所得结果未必直观或精确。

特定形式模型

通常在应用非线性模型时，会选用特定类型的模型。这使得拟合更便捷或在数值上更易处理，但需与领域专家达成一致，并且与数据相匹配。所谓特定类型，是指使用预定义的数学方程。以下列出最常用的一些非线性方程（为简洁，部分略去）。

非线性回归模型的策略

1. 可视化数据，观察是否有某个特定数学函数能较好地解释数据。
2. 通过目测图形或使用自启动函数（见下文）选择初始参数值。
3. 拟合模型并对其进行检视。
4. 对拟合模型做统计分析，了解参数置信区间、残差误差以及解释的变异度。
5. 如有助于降低残差误差，尝试为模型添加更多参数。
6. 务必同时拟合一个线性模型，并将非线性与线性模型结果加以比较。

通过示例解释非线性模型

我们将针对一个常用数据集 Puromycin 拟合一个特定的非线性模型。该数据集展示了涉及 Puromycin 药物的细胞酶促反应中，反应速率与底物浓度的信息。数据包含两类细胞：一类接受 Puromycin 处理，另一类未处理。

names(Puromycin)Puromycin'conc''rate''state'技术支持

该数据框中的三个变量分别为底物浓度、反应初速率，以及是否接受处理的指示变量。

将速率对浓度作图，并标注两种处理水平，可以看到一个总体模式：曲线呈渐近趋势，且处理组与未处理组的轨迹明显不同。

plot(Puromycin$conc, Puromycin$rate, type="n", xlab = "conc", ylab = "rate")text(Puromycin$conc, Puromycin$rate,      ifelse(Puromycin$state == "treated", "T", "U"))技术支持

由于该数据在反应速率与浓度之间呈现出明显的 Michaelis-Menten 关系，实验者预计可拟合如下数学方程。

其中 E 为实验误差。借鉴线性模型的思路，通过特定的非线性模型，我们尝试估计参数 Vmax 与 K，使残差平方和最小化（当然也可采用其他方法估计这些参数）：

拟合非线性模型

即便我们知道二者的关系，拟合非线性模型也并非易事，因为模型需要为 V_max 与 K 给出初始猜测。如果初值不佳，模型可能会拟合得很差。通常有两种选择：

• 通过目测图形选取初始值。
• 使用自启动函数（如下所示）。

对于每一种模型类型，都有相应的自启动函数可用于初值猜测。

使用初始猜测拟合模型

通过目测图形，我们选取 Vmax = 160、K = 0.05 作为初值，可用 R 函数 nls() 对数据进行拟合。在线性回归中通常无需显式指定 Vm 或 K，但非线性模型不同。任何迭代算法都需要良好的起点，否则可能无法收敛。下面展示如何在 Puromycin 上使用初值拟合非线性模型。

Purboth_1 <- nls(rate ~ (Vm)*conc/(K+conc), Puromycin, list(Vm=160, K=0.05))summary(Purboth_1)Formula: rate ~ (Vm) * conc/(K + conc)Parameters:Estimate Std. Error t value Pr(>|t|)Vm 190.80624 8.76459 21.770 6.84e-16 ***K 0.06039 0.01077 5.608 1.45e-05 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 18.61 on 21 degrees of freedomNumber of iterations to convergence: 5Achieved convergence tolerance: 5.121e-06技术支持

计算残差误差

残差误差计算如下：

sse <- Purboth_1$m$deviance()sse7276.54697931423技术支持

计算解释变异度

要计算该模型解释的变异度（相当于线性模型术语中的 R 方），需先拟合空模型以估计总体均值，从空模型摘要中提取总平方和 sst，然后据此计算解释变异度的百分比。

sse <- Purboth_1$m$deviance()null <- lm(rate~1, Puromycin)sst <- data.frame(summary.aov(null)[[1]])$Sum.Sqpercent_variation_explained = 100*(sst-sse)/sstpercent_variation_explained85.3488323994681技术支持

该模型解释了 85.34% 的数据变异。

使用自启动函数拟合模型

由于我们使用的是 Michaelis-Menten 模型，可以改用自启动函数 SSmicmen()。不过，结果相同。

Purboth_self <- nls(rate~SSmicmen(conc,Vm,K), Puromycin)summary(Purboth_self)Formula: rate ~ SSmicmen(conc, Vm, K)Parameters:Estimate Std. Error t value Pr(>|t|)Vm 190.80667 8.76462 21.770 6.84e-16 ***K 0.06039 0.01077 5.608 1.45e-05 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 18.61 on 21 degrees of freedomNumber of iterations to convergence: 4Achieved convergence tolerance: 6.745e-06技术支持

计算残差误差

sse <- Purboth_self$m$deviance()sse7276.54697947673技术支持

计算解释变异度

可进一步检查模型解释的变异度。

sse <- Purboth_self$m$deviance()null <- lm(rate~1, Puromycin)sst <- data.frame(summary.aov(null)[[1]])$Sum.Sqpercent_variation_explained = 100*(sst-sse)/sstpercent_variation_explained85.348832399141技术支持

总体而言，拟合得到的速率与浓度之间关系的方程如下：

可视化拟合效果

可在数据上叠加拟合曲线来观察模型表现，该曲线呈现非线性。

plot(Puromycin$conc, Puromycin$rate, xlab="conc", ylab = "rate")lines(Puromycin[Puromycin$state == "treated",]$conc,       predict(Purboth_1, Puromycin[Puromycin$state == "treated",]), col=2)lines(Puromycin[Puromycin$state != "treated",]$conc,       predict(Purboth_1, Puromycin[Puromycin$state != "treated",]), col=2)技术支持

更优模型选择

可以通过增加参数进一步调优模型。观察速率与浓度的散点图可注意到，不同类别（处理 vs 未处理）走向不同，两条轨迹在 y 轴上的差异大约为 40。

plot(Puromycin$conc, Puromycin$rate, type="n", xlab = "conc", ylab = "rate")text(Puromycin$conc, Puromycin$rate,      ifelse(Puromycin$state == "treated", "+", "*"))技术支持

本质上，我们尝试在条件变量 state 上进行分层拟合，得到如下非线性方程。

引入额外参数 delV 后，模型表现更佳。我们将 delV 初值设为 40，用包含三个参数的 nls 进行拟合。

Purboth <- nls(rate ~ (Vm + delV*(state=="treated"))*conc/(K+conc), Puromycin,                list(Vm=160, delV=40, K=0.05))summary(Purboth)Formula: rate ~ (Vm + delV * (state == "treated")) * conc/(K + conc)Parameters:      Estimate Std. Error t value Pr(>|t|)    Vm   166.60397    5.80742  28.688  < 2e-16 ***delV  42.02591    6.27214   6.700 1.61e-06 ***K      0.05797    0.00591   9.809 4.37e-09 ***---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 10.59 on 20 degrees of freedomNumber of iterations to convergence: 5 Achieved convergence tolerance: 8.93e-06技术支持

摘要显示三个参数均具有统计显著性。

计算残差误差

该模型的残差误差显著降低：

sse <- Purboth$m$deviance()sse2240.89143883885技术支持

计算解释变异度

该模型解释的变异度更多。

sse <- Purboth$m$deviance()null <- lm(rate~1, Puromycin)sst <- data.frame(summary.aov(null)[[1]])$Sum.Sqpercent_variation_explained = 100*(sst-sse)/sstpercent_variation_explained95.4880142822744技术支持

AIC 检查

也可以对两个模型做 AIC 检查。尽管如今更常用交叉验证，但那通常需要更多数据。一般而言，AIC 越低模型越好。

AIC(Purboth)AIC(Purboth_1)178.591273238391203.680274951742技术支持

可视化拟合效果

通过在数据上绘制拟合曲线，可以直观看到调优后的模型表现；不出意料，调优模型对不同类别的拟合更好。

plot(Puromycin$conc, Puromycin$rate, xlab="conc", ylab = "rate")lines(Puromycin[Puromycin$state == "treated",]$conc, predict(Purboth, Puromycin[Puromycin$state == "treated",]), col=2)lines(Puromycin[Puromycin$state != "treated",]$conc, predict(Purboth, Puromycin[Puromycin$state != "treated",]), col=3)legend(0, 200, legend = c("treated", "untreated"), fill = c(2,3))技术支持

参数的置信区间

不应仅依赖 p 值。95% 置信区间给出了每个显著参数在 2.5% 与 97.5% 处的系数取值范围，且均不包含 0。

confint(Purboth)技术支持

利用置信区间，可以统计地说明：例如 Vmax 的真值有 95% 的概率落在上述区间内，对 Vmax 而言为 (154.62, 179.25)。

剖面分析

还可通过直接考察目标函数，对模型进行剖面分析，以理解参数估计的不确定性。在本模型中，我们有三个参数 Vmax、delV、K，可对其进行剖面分析：

• 获得置信区间；
• 并使用类似于线性模型 t 统计量的剖面 t 函数。剖面 t 函数的图像给出单个参数的似然区间，并说明估计过程的非线性程度。

绘制置信区间

该图绘制了在不同取值下参数的置信区间（对应不同的置信水平：99%、95%、90%、80% 与 50%）。其与 confint(object, parm, level = 0.95, ...) 函数给出的 95% 置信区间相似。

Purbothpf = profile(Purboth)par( mfrow= c(2,2) )plot(Purbothpf)技术支持

绘制 _T 与参数的关系

将三个参数 Vmax、delV、K 与 _T 作图，可感知拟合的线性或非线性程度。可以观察到，delV 与 Vm 的曲线近似线性，而 K 略有弯曲。这在预期之中，因为 K 与响应之间呈现非线性关系。

dataVm = data.frame(Purbothpf$Vm$par.vals)datadelV = data.frame(Purbothpf$delV$par.vals)dataK = data.frame(Purbothpf$K$par.vals)par( mfrow= c(2,2) )plot(dataVm$Vm,Purbothpf$Vm$tau, type="l", xlab = "Vm", ylab = "Tau")plot(datadelV$delV,Purbothpf$delV$tau, type="l", xlab = "delV", ylab = "Tau")plot(dataK$K,Purbothpf$K$tau, type="l", xlab = "K", ylab = "Tau")技术支持

拟合线性模型

也可使用同一个 R 库函数 nls() 来拟合线性模型，如下所示。笔者认为，为了可解释性与比较，始终尝试线性模型是个好习惯。

PurbothRealLinear <- nls(rate ~ i + m*conc, Puromycin, list(i=160, m=0.05) )summary(PurbothRealLinear)Formula: rate ~ i + m * concParameters:  Estimate Std. Error t value Pr(>|t|)    i    93.92       8.00   11.74 1.09e-10 ***m   105.40      16.92    6.23 3.53e-06 ***---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 28.82 on 21 degrees of freedomNumber of iterations to convergence: 2 Achieved convergence tolerance: 2.571e-09技术支持

可视化拟合效果

线性模型表现不佳，可在数据上绘制拟合直线，并与其非线性对应模型进行对比。

plot(Puromycin$conc, Puromycin$rate, xlab="conc", ylab = "rate")lines(Puromycin$conc, predict(PurbothRealLinear), col=4)技术支持

非线性与线性模型的分析性比较

从上述图中可知，线性模型表现较差，但应如何比较这些模型？这如同“苹果与橘子”的比较，无法直接使用 AIC 结果。

拟合优度

何谓良好的回归拟合——即预测值与真实值尽可能接近。因此，我们可以用 R 的 lm() 函数来检验模型预测的质量并进行比较。更准确的模型，其点将更接近理想参考线。理想线表示 100% 准确，即每个数据点的预测值都等于真实值，因此所有点将落在同一直线上。

非线性拟合

对每个输入数据点，使用模型预测目标值，然后按如下方式应用 lm()。通过绘制真实速率（x 轴）与预测速率（y 轴）的散点图，并叠加理想线，可以观察点到该线的接近程度。

nonLinear_df = data.frame("rate"=Puromycin$rate, "pred_rate"= predict(Purboth))fit_nonLinear = lm(nonLinear_df$rate ~ nonLinear_df$pred_rate)AIC(fit_nonLinear)plot(nonLinear_df$rate, nonLinear_df$pred_rate, xlab="actual_rate",      ylab = "pred_rate",      main = "Comparing Actual vs Predicted for Non Linear Model")lines(Puromycin$rate, Puromycin$rate)174.711198508104技术支持

线性拟合

与下方的线性图相比，上方的非线性图明显拟合更好。

linear_df = data.frame("rate"=Puromycin$rate, "pred_rate"= predict(PurbothRealLinear))fit_linear = lm(linear_df$rate~linear_df$pred_rate)AIC(fit_linear)plot(linear_df$rate, linear_df$pred_rate, xlab="actual_rate", ylab = "pred_rate",      main= "Comparing Actual vs Predicted using lm()")lines(Puromycin$rate, Puromycin$rate)223.783416979551技术支持

改进线性模型

直接用线性模型可能不公平，因为数据显然不是直线关系。我们可尝试在变换后的数据上拟合模型。此处选择对 conc 使用 sqrt() 变换，并拟合线性模型如下。

PurbothLinearSqrt = nls(rate ~ i + m*sqrt(conc), Puromycin, list(i=160, m=0.05) )summary(PurbothLinearSqrt)Formula: rate ~ i + m * sqrt(conc)Parameters:  Estimate Std. Error t value Pr(>|t|)    i   60.961      8.758   6.961 7.11e-07 ***m  139.134     15.675   8.876 1.50e-08 ***---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 22.31 on 21 degrees of freedomNumber of iterations to convergence: 2 Achieved convergence tolerance: 1.838e-09技术支持

残差标准误确实有所下降，但低于预期幅度。

绘制拟合效果

可以清楚地看到，对 conc 使用 sqrt() 变换能改进线性模型。然而，与非线性模型 Purboth 相比仍不占优。PurbothLinearSqrt 的残差标准误几乎是 Purboth 的两倍，因此我们不再深入分析。也可尝试更高阶多项式项（多项式回归），但微调所需工作量更大，而使用非线性模型往往更易取得更好结果。

plot(Puromycin$conc, Puromycin$rate, xlab="conc", ylab = "rate")lines(Puromycin[Puromycin$state != "treated",]$conc, predict(PurbothLinearSqrt, Puromycin[Puromycin$state != "treated",]), col=2)lines(Puromycin[Puromycin$state != "treated",]$conc, predict(PurbothLinearSqrt, Puromycin[Puromycin$state != "treated",]), col=3)技术支持

“解开”非线性方程

还可以用一种技巧将非线性函数“解开”为线性模型。此法在手工计算时代较常见，但如今往往难以获得更好效果。

将非线性转换为线性形式

例如，可对 Michaelis-Menten 方程取倒数，得到如下形式，并据此应用线性模型。

我们仍可使用 nls() 来拟合新的方程。

PurbothPseudoLinear <- nls((1/rate) ~ K /(Vm*conc) + 1/Vm , Puromycin,                            list(Vm=160, K=0.05))summary(PurbothPseudoLinear)Formula: (1/rate) ~ K/(Vm * conc) + 1/VmParameters:    Estimate Std. Error t value Pr(>|t|)    Vm 1.674e+02  1.362e+01  12.287 4.70e-11 ***K  3.900e-02  6.172e-03   6.318 2.89e-06 ***---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.001786 on 21 degrees of freedomNumber of iterations to convergence: 3 Achieved convergence tolerance: 2.997e-09技术支持

虽然残差标准误很低，但所用的是经变换的数据，因而无法与非线性模型直接比较。更可取的方法是可视化两种模型在原始数据集上的表现。

绘制拟合效果

我们将只含两个参数 Vmax、K 的 Purboth_1 与 PurbothPseudoLinear 进行比较。图示清楚地表明，非线性模型效果更好。

plot(Puromycin$conc, Puromycin$rate, xlab="conc", ylab = "rate",      main = "Scatterplot conc vs rate")lines(Puromycin[Puromycin$state == "treated",]$conc,       predict(Purboth_1, Puromycin[Puromycin$state == "treated",]), col=2)lines(Puromycin[Puromycin$state != "treated",]$conc,       predict(Purboth_1, Puromycin[Puromycin$state != "treated",]), col=2)lines(Puromycin[Puromycin$state == "treated",]$conc,       1/predict(PurbothPseudoLinear, Puromycin[Puromycin$state == "treated",])      , col=4)lines(Puromycin[Puromycin$state != "treated",]$conc,       1/predict(PurbothPseudoLinear, Puromycin[Puromycin$state != "treated",])      , col=4)legend(0, 200, legend = c("nls", "nls-turned-linear"), fill = c(2,4))技术支持

拟合优度

可比较模型的拟合优度。

nonLinear_df = data.frame("rate"=Puromycin$rate, "pred_rate"= predict(Purboth_1))fit_nonLinear = lm(nonLinear_df$rate ~ nonLinear_df$pred_rate)AIC(fit_nonLinear)203.017860461858技术支持

linear_df = data.frame("rate"=Puromycin$rate,                        "pred_rate"= 1/predict(PurbothPseudoLinear))fit_linear = lm(linear_df$rate~linear_df$pred_rate)AIC(fit_linear)206.999538212925技术支持

这表明使用非线性函数的模型效果更佳。

par( mfrow= c(1,2) )plot(nonLinear_df$rate, nonLinear_df$pred_rate, xlab="actual_rate",      ylab = "pred_rate",      main = "Non-Linear")lines(Puromycin$rate, Puromycin$rate)plot(linear_df$rate, linear_df$pred_rate, xlab="actual_rate", ylab = "pred_rate",      main = "Turned Linear")lines(Puromycin$rate, Puromycin$rate)技术支持

结论

在机制性数据场景下，非线性模型通常优于线性模型。然而，非线性模型并不能保证估计问题一定有数值解。应通过作图选择最能解释数据关系的数学函数。在这类情况下，非线性模型可能会胜过线性模型。

如果本文对您有所启发，并希望进一步探索回归模型，欢迎查看 DataCamp 的R 中级回归课程。