#rstats 生长曲线、生长方程、逻辑斯蒂与nls R_square

我本身不做生长方程拟合，而身边常有。拟合模型时在初始值选择上常遇见困难，时而成功，时而报错。该怎么选择初始值？

更新：经@木䬕的提示和参考3中的例子，发现R可以自动给初始值，如逻辑斯蒂用SSlogis()，渐进回归模型SSasymp()。

# SSlogis
nlsfitSS <- nls(height ~ SSlogis(age, Asym, xmid, scal),
                data=Loblolly)
# SSasymp
nlsfitSS <- nls(height ~ SSasymp(age, Asym, R0, lrc),
                data=Loblolly)
# 0.993144

SSlogis()的模型形式与正文中给出的不同，但结果是一样的。
$Y = \frac{a}{1+e^{(b-time)/c}}$
$a$ ， $b$ ， $c$ 分别对应Asym，xmid和scal。
SSasymp()的模型形式为
$Y = a+(b-a)*e^{-exp(c)*time}$
$a$ ， $b$ ， $c$ 分别对应Asym，R0和lrc。
另外，R还提供SSasympOff, SSasympOrig, SSbiexp, SSfol, SSfpl, SSgompertz, SSmicmen, SSweibull这些函数供使用。

介绍

生长方程是森林经理上生物量估计，或是数量遗传学中功能作图用到的几个可以拟合生物生长过程函数。常见的有：

理查德
$Y = Y_{max}*(1-e^{-b*time})^c$

$Y_{max}$ 【有些写作 $a$ 】：渐近线，即成年稳定值，停止生长，理论上的最大值；
$b$ ：增长参数，描述接近渐近线的速度，越大说明单位时间内生长得越快；
$c$ ：没有特别含义。

逻辑斯蒂
$Y = \frac{a}{1+e^{-(b+c*time)}}$

参数意义与理查德近似。

数据

# R中自带数据Loblolly
data(Loblolly)
# 散点图查看关系
plot(Loblolly$height~Loblolly$age)

图1. height~age.png

构建方程函数

# 理查德
f_richard <- function(x, a, b, c){
  return(a*(1-exp(-b*x))^c)
}
# 逻辑斯蒂
f_logist <- function(x, a, b, c){
  return(a/(1+exp(-(b+c*x))))
}

拟合

理查德
理查德初始值没有好的方法，只能试，但不是随意的试，要按照一定规则：知道 $a$ 是渐近线，通过图1可以看到大致在60的位置，所以设置为60； $b$ 是速率，要大于0，一般设置为0.1， $c$ 没有太大限制。然后可以通过画图来查看你给定的这一组初始值效果：

# 理查德初始值
plot(Loblolly$height~Loblolly$age)
curve(f_richard(x, a=60, b=0.1, c=2), add=T)

图2. 理查德初始值

效果还可以，那么运行一下模型：

m_richard <- nls(height ~ f_richard(age, a, b, c),
                 data=Loblolly,
                 start=list(a=60, b=0.1, c=2))
summary(m_richard)
# Formula: height ~ f_richard(age, a, b, c)
# 
# Parameters:
#   Estimate Std. Error t value Pr(>|t|)    
# a 76.933519   2.580073   29.82   <2e-16 ***
#   b  0.082659   0.006067   13.62   <2e-16 ***
#   c  1.846938   0.093736   19.70   <2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 1.731 on 81 degrees of freedom
# 
# Number of iterations to convergence: 4 
# Achieved convergence tolerance: 9.319e-07

逻辑斯蒂
逻辑斯蒂的初始值设置相对容易些，先拟合一个对响应变量log转换的线性模型，将得到的系数视为初始值：

coef(lm(log(height/60)~age,data=Loblolly))
# (Intercept)         age 
# -2.409599    0.111560

这里60也是渐近线值，b和c分别是-2.4和0.1,。然后拟合模型：

m_logist <- nls(height ~ f_logist(age, a, b, c),
                data=Loblolly,
                start=list(a=60, b=-2.4, c=0.1))
summary(m_logist)
# Formula: height ~ f_logist(age, a, b, c)
# 
# Parameters:
#   Estimate Std. Error t value Pr(>|t|)    
# a 61.344070   0.984761   62.29   <2e-16 ***
#   b -2.731120   0.083262  -32.80   <2e-16 ***
#   c  0.231854   0.009177   25.27   <2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 2.657 on 81 degrees of freedom
# 
# Number of iterations to convergence: 10 
# Achieved convergence tolerance: 1.234e-06

都可以成功运行。
最终的拟合效果间图3。

图3. 拟合效果。红-理查德；蓝-逻辑斯蒂

模型评估 $R^2$

quasi.rsq.nls <- function (mdl, y, param){
  adj <- (sum(!is.na(y)) - 1)/(sum(!is.na(y)) - param)
  sum.sq <- (sum(!is.na(y)) - 1) * var(y, na.rm = TRUE)
  rsq <- 1 - (adj * (deviance(mdl)/sum.sq))
  return(rsq)
}
# sample
quasi.rsq.nls(m_richard, Loblolly$height, param=3)
quasi.rsq.nls(m_logist, Loblolly$height, param=3)
#0.9929895
#0.9834879

理查德略优。

参考：

最后编辑于：2024.05.29 15:22:56