响应曲线初始值的确定(待补充)

祝介东

2018/09/06

让人烦恼的初始值

又是一年生长季临近结束时,不少小伙伴们是第一次用 LI-6400XT 或 LI-6800,好不容易通过各种渠道掌握了操作,又历尽艰辛,经过一个生长季风吹日晒(做光合的挨雨淋的机会比做其他测量的少,听不懂但又要接着用光合仪的小伙伴们赶紧面壁翻说明书翻资料去。)的洗礼,终于到分析数据的时候了,在面对响应曲线数据时,很多小伙伴有点懵了,怎么拟合(其实也不怪小伙伴们,我的印象里生物统计基本上就讲到简单线性回归,多元线性回归、非线性拟合基本不讲或跳过),好不容易看文献资料明白了光响应曲线、二氧化碳响应曲线模型是怎么回事,选好了使用的模型,照着教程开始拟合时候傻眼了——初始值是什么东西?怎么确定?能不能不要初始值?下面的内容我尽量用最简单的语言解释一下(如果照抄教材的内容,你们也没兴趣看那一堆公式和晦涩语言写的东西,否则也不会看我写的内容了)。

通俗易懂但不严谨的解释

在解释初始值之前我们首先需要了解一个数学上的概念——迭代

“迭代法”也称“辗转法”是一种不断用变量的旧值递推新值的过程。

用通俗但不是特别严谨的说法可解释为:每次执行这种算法是,程序都会从原值(也就是我抄的上面迭代法定义的旧值)推出一个新值。

之所以先介绍这个迭代,原因很简单,非线性拟合就是通过迭代的方法,需要对每一个变量最初的估计值进行不断的迭代,得到一个向一个点收缩或汇聚的值,这个估计值必须在实际值的一定范围内,程序通过不断调整这个值来改善拟合结果。这就解释了上面的问题,初始值是让程序开始运行的前提,不然没法迭代,必须设定。至于怎么确定初始值,这件事情解释起来要麻烦一些,请接着看下面的内容。

初始值(范围)确定的方法

我下面的内容将以 LI-6800 的光响应曲线的测试数据,使用非直角双曲线模型进行拟合来讲解具体的 R 中的一些实现方法,我们首先导入数据,然后再利用这些数据逐个举例不同的确定初始值的方式。

lrc <- read.csv("data/lrc.csv")


# 光响应曲线比较简单,我们将需要的数据直接提取,方便后面操作
lrc_Q <- lrc$Qin
lrc_A <- lrc$A

nlsLM 解决方案

nlsLM 来自于 Elzhov et al. (2016)minpack.lm,利用 C 语言的 MINPACK 库,修改了 Levenberg-Marquardt 算法,在实际操作中,很多时候并不准确的输入初始值,他也能得出比较好的拟合结果。但结果未必完美,出现下面让人烦恼的报错:

singular gradient matrix at initial parameter estimates

的概率会大大降低,而且尽管结果不如意,我们也可以利用他的结果缩小初始值的范围,继续尝试其他初始值。

例如下面的例子中,非直角双曲线的 Rd 的初始值我们可以利用暗呼吸的实测值大致估计,同理最大光合速率也是如此,剩下的分别为非直角双曲线曲率,我们暂定为 1,alpha 也暂定为 0.1,使用 nlsLM 进行拟合,结果如下:

library(minpack.lm)

lrcnls_lm <- nlsLM(lrc_A ~ (1/(2*theta))*
        (alpha*lrc_Q+Am-sqrt((alpha*lrc_Q+Am)^2 - 
        4*alpha*theta*Am*lrc_Q))- 
        Rd, start=list(Am=(max(lrc_A)-min(lrc_A)),
        alpha=0.1,Rd=-min(lrc_A),theta=0.8)) 

结果没有报错,看上去没有问题,那我们观察一下具体的拟合结果:

summary(lrcnls_lm)
## 
## Formula: lrc_A ~ (1/(2 * theta)) * (alpha * lrc_Q + Am - sqrt((alpha * 
##     lrc_Q + Am)^2 - 4 * alpha * theta * Am * lrc_Q)) - Rd
## 
## Parameters:
##        Estimate Std. Error t value Pr(>|t|)    
## Am    12.307570   0.406739  30.259 2.30e-10 ***
## alpha  0.045706   0.003423  13.352 3.09e-07 ***
## Rd     0.656638   0.132646   4.950 0.000791 ***
## theta  0.707522   0.079738   8.873 9.59e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1852 on 9 degrees of freedom
## 
## Number of iterations to convergence: 8 
## Achieved convergence tolerance: 1.49e-08

结果看上去还可以1

作图比对法

模型很多参数可以用已有数据去估计,我们可以只来分析难以判断的参数,流程如下:

  • Rd、Am等我们可以利用测量值来确定一个范围。
  • 剩余的参数,我们也可以根据经验或文献来有一个大致的判断。
  • 然后我们根据数学的方式来判断哪个参数对曲线形状影响最大(例如在分母上的参数,或者是乘以该参数,该参数可以显著改变计算结果,例如整体乘以或除以 0.1 还是 0.01,像 Rd 之类的参数本身就很小,多数公式中都是减去该值,对结果影响很小,我们通常直接使用实测值 )。
  • 将该参数取一系列值带入模型来求解净光合速率。
  • 将计算的A值与光强进行作图,看我们计算的曲线与测量数据点的重合程度,必要时在修改其他参数,使曲线和散点重合度最好,重合程度最高的参数值即为我们需要的初始值。

实现过程

# 我们选择的模型,将其写为一个函数,用于计算净光合速率
expfct <- function(x, Am, alpha, Rd, theta) {(1/(2 * theta)) * (alpha * x + Am - sqrt((alpha * x + Am)^2 - 4 * alpha * theta * Am * x)) - Rd
}

# 我们的数据
test <- data.frame(x = lrc_Q, y = lrc_A)

# 先做实测数据的散点图
plot(y ~ x, data = test)

# 利用上面的函数,假定 alpha 的值为0.8,看计算值与测量值重合程度
curve(expfct(x, Am = (max(lrc_A)-min(lrc_A)),
     alpha=0.8, Rd=-min(lrc_A), theta=0.8), add = TRUE
             )

观察上图的结果可以看到,曲线在 0-600 的范围内,拟合值明显偏大,观察模型的方程式,以及其他起始值的设定方式,我们初步判断 alpha 的值偏大,于是乎我们将其改小观察,观察曲线和测量点的重合仍然不是很好,我们尝试修改 theta 值与 alpha 值(也即曲线高于测量点,则需要减小纵坐标的值,低于测量点,则需要增加该值,该过程省略,我大概设置了五分钟完成),最终得出的结果如下:

plot(y ~ x, data = test)
curve(expfct(x, Am = (max(lrc_A)-min(lrc_A)),
     alpha=0.06, Rd=-min(lrc_A), theta=0.82), add = TRUE)

尽管看上去效果仍然不满意,但我们可试着进行拟合,看能否得到显著差异的结果:

lrcnls_manual <- nls(lrc_A ~ 
        (1/(2*theta))*
        (alpha*lrc_Q+Am-sqrt((alpha*lrc_Q+Am)^2 - 
                               4*alpha*theta*Am*lrc_Q))- 
        Rd, start=list(Am=(max(lrc_A)-min(lrc_A)),
                       alpha=0.03,Rd=-min(lrc_A),theta=0.6))
summary(lrcnls_manual)
## 
## Formula: lrc_A ~ (1/(2 * theta)) * (alpha * lrc_Q + Am - sqrt((alpha * 
##     lrc_Q + Am)^2 - 4 * alpha * theta * Am * lrc_Q)) - Rd
## 
## Parameters:
##        Estimate Std. Error t value Pr(>|t|)    
## Am    12.307585   0.406741  30.259 2.30e-10 ***
## alpha  0.045706   0.003423  13.352 3.09e-07 ***
## Rd     0.656642   0.132646   4.950 0.000791 ***
## theta  0.707518   0.079739   8.873 9.59e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1852 on 9 degrees of freedom
## 
## Number of iterations to convergence: 7 
## Achieved convergence tolerance: 4.601e-06
# 对拟合之后的结果作图,观察使用我们的估计值,
# 迭代的最终值与元数据的重合程度
plot(y ~ x, data = test, ylim = c(-2, 14))
curve(expfct(x, Am = 12.307586,
     alpha=0.045706, Rd= 0.656643, theta=0.707518), add = TRUE)

从最终图形的呈现以及 F 检验的 p 值来讲,图形已经比较完美了。也就是说尽管我们作图的时候看到重合度并不高,但是非线性拟合本来就是一个迭代的过程,只要我们的数据与真实值相差不大,还是能够得到完美结果的

直观展示

上面的表述太啰嗦,直接用下面的图形说明一下,其中 alhpa 的取值在此处选择从 0.01 到 0.07,每次增加 0.05,其他值分别为 Am = 12.31, Rd= 0.66, theta=0.71 (此处为展示效果和方便,将这些值直接按照拟合结果设定了,实际差别不大)

library(ggplot2)
library(purrr)

# 光响应曲线比较简单,我们将需要的数据直接提取,方便后面操作
lrc_Q <- lrc$Qin
lrc_A <- lrc$A
n <- length(lrc_A)


alp <- paste0("a=", seq(0.01, 0.07, by = 0.005))
alpn <- rep(alp, each = n)

expfct <- function(x, Am, alpha, Rd, theta) {(1/(2 * theta)) * (alpha * x + Am - sqrt((alpha * x + Am)^2 - 4 * alpha * theta * Am * x)) - Rd
}

paras <- data.frame(alpha = rep(seq(0.01, 0.07, by = 0.005), each = n), 
           x = rep(lrc_Q, n), Am = rep(12.31, n), Rd = rep(0.66, n), 
           theta = rep(0.71, n))
y = unlist(pmap(paras, expfct))

show <- data.frame(x = rep(lrc_Q, 14),
           y = c(lrc_A, y), 
           a = factor(c(rep("measured", n), alpn),
           level = c("measured", alp)
             ))

ggplot(data = show, aes(x, y, group = a, color=a)) + 
  geom_point() + 
  geom_smooth(se = FALSE) 

从上图我们我们可以看到,实测值在 alpha =0.04 和 alpha = 0.05 两条曲线之间,在 0.045 时最接近测量点,也就是我们把初始值设为 0.04 和 0.05 之间最接近,本例中可认为是0.045,实际这三个值均可。

自动多次尝试法

该方法实际为使用 nls2 来实现,具体方法参考 Bouvier and Huet (1994) 的文章,可简单概括为使用一系列的起始值梯度(例如下面的代码中, alpha 的取值在 0.01 到 0.08 之间 ),然后软件循序使用不同的起始值,即排列组合所有的起始值序列,最终找到合适的值,具体实现如下:

library(nls2)
## Loading required package: proto
grid.test <- expand.grid(list(
  Am=c(12),
  alpha = seq(0.01, 0.08, by =0.01),
  Rd = seq(0, 3),
  theta=seq(0.1, 1, by = 0.1)
  ))

lrcnls2 <- nls2(lrc_A ~ 
        (1/(2*theta))*
        (alpha*lrc_Q+Am-sqrt((alpha*lrc_Q+Am)^2 - 
                               4*alpha*theta*Am*lrc_Q))- 
        Rd, start = grid.test, algorithm = "brute-force")
summary(lrcnls2)
## 
## Formula: lrc_A ~ (1/(2 * theta)) * (alpha * lrc_Q + Am - sqrt((alpha * 
##     lrc_Q + Am)^2 - 4 * alpha * theta * Am * lrc_Q)) - Rd
## 
## Parameters:
##        Estimate Std. Error t value Pr(>|t|)    
## Am    12.000000   0.623023  19.261 1.27e-08 ***
## alpha  0.050000   0.006414   7.795 2.72e-05 ***
## Rd     1.000000   0.260153   3.844  0.00394 ** 
## theta  0.800000   0.102143   7.832 2.62e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3663 on 9 degrees of freedom
## 
## Number of iterations to convergence: 320 
## Achieved convergence tolerance: NA

通过结果可以看到,虽然和之前采用手动方法判定的结果比较接近,但是还是略有差异,可以看一下他们各自的结果同测量值的重合程度:

plot(lrc_Q, lrc_A)
lines(lrc_Q,  predict(lrcnls2), col="red")
lines(lrc_Q,  predict(lrcnls_manual), col="blue")

可以看到,使用 nls2 的拟合结果似乎和测量值更匹配,当然这只是第一印象,后续的判断还要进一步通过 F 检验、 AIC、BIC 等统计方式才能判定。

小结

采用如上三种方式都可以有效的解决起始值的问题,nlsLM 操作上更易实现,对初始值的大小不敏感,但设置不能太离谱,否则仍然会报错。作图比对法操作上更麻烦一些,但是这种方式一定能得出合理的初始值设置。采用 nls2 类似于将手动作图方式自动化,类似于 SPSS 中非线性拟合中需要给出一个初始值的范围,且该范围不能过大。如有一定的经验,操作起来将非常迅速。或者使用作图法确定大致的范围,将该范围输入到 nls2 中,这样会节省时间,也更加方便。

需要注意的是,这三种方法结合起来使用会更好,例如,即使使用 nlsLM 的结果不合理,也可以参考他们参数的范围(部分结果也可能是差异显著),然后将这些结果用于手动作图判定参数或者 nls2 中判定参数范围。

声明

以上内容仅为个人理解,如有错误欢迎指出

参考文献

Bouvier A, Huet S. 1994. Nls2 non-linear regression by s-plus—functions. Computational Statistics & Data Analysis 18, 187–190.

Elzhov TV, Mullen KM, Spiess AN, Bolker B. 2016. Minpack.lm: R interface to the levenberg-marquardt nonlinear least-squares algorithm found in minpack, plus support for bounds.


  1. 有些时候结果并不理想,该方法并不是万无一失