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

---

更新：经@木䬕的提示和参考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()的模型形式与正文中给出的不同，但结果是一样的。

，，分别对应Asym，xmid和scal。
SSasymp()的模型形式为

，，分别对应Asym，R0和lrc。
另外，R还提供SSasympOff, SSasympOrig, SSbiexp, SSfol, SSfpl, SSgompertz, SSmicmen, SSweibull这些函数供使用。

---

介绍

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

1. 理查德
   

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

2. 逻辑斯蒂
   

• 参数意义与理查德近似。

数据

# 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))))
}

拟合

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

# 理查德初始值
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

2. 逻辑斯蒂
   逻辑斯蒂的初始值设置相对容易些，先拟合一个对响应变量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. 拟合效果。红-理查德；蓝-逻辑斯蒂

模型评估

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

理查德略优。

参考：

1. Modeling Logistic Growth Data in R
2. Fitting Non-Linear Growth Curves in R
3. Non-linear regression
4. quasi.rsq.nls