引言

采用Richards、Weibull、Logistic、Korf、Mitscherlich、Schumacher6种理论生长方程对某一树种单木树高生长过程进行模拟，筛选出拟合精度最佳的生长方程。 理论生长方程如下

#建模数据

以小兴安岭地区一株265年的天然红松林解析木树高生长数据为例：

age<-c(seq(0,250,10),256)
height<-c(0,0.80,2.20,4.00,5.20,7.70,8.95,10.40,12.00,13.90,16.40,18.10,19.40,20.80,
          22.05,23.02,24.00,24.95,25.90,26.80,27.62,28.27,28.80,29.30,29.80,30.30,30.60)
ht.growth<-data.frame(age,height)

#绘制散点图

plot(height~age,data=ht.growth,xlab="Age (year)",ylab="Tree height (m)")

#构建模型评价函数

首先构建mpc函数，通过MSE,RMSE,Rsquare,MAE,MAPE,RASE,AIC,BIC8个指标来用来评价和筛选模型。

library(modelr)
mpc<-function(model,data){
  MSE<-modelr::mse(model=model,data=data)
  RMSE<-modelr::rmse(model=model,data=data)
  Rsquare<-modelr::rsquare(model=model,data=data)
  MAE<-modelr::mae(model=model,data=data)
  MAPE<-modelr::mape(model=model,data=data)
  RASE<-modelr::rsae(model=model,data=data)
  AIC=AIC(object=model)
  BIC=AIC(object=model)
  return(data.frame(MSE=MSE,RMSE=RMSE,Rsquare=Rsquare, MAE=MAE,MAPE=MAPE,RASE=RASE,
                    AIC=AIC,BIC=BIC))
}

##直线回归

利用R语言自带的lm函数拟合线性模型

mod.lm<-lm(height~age,data=ht.growth)
mod.lm
## 
## Call:
## lm(formula = height ~ age, data = ht.growth)
## 
## Coefficients:
## (Intercept)          age  
##      1.7736       0.1265
mpc(mod.lm,data=ht.growth)
##        MSE     RMSE   Rsquare      MAE MAPE       RASE      AIC      BIC
## 1 3.460883 1.860345 0.9653478 1.609177  Inf 0.08844152 116.1438 116.1438

##非线性回归

利用mipack.lm 程序包中的nlsLM函数，进行非线性模型拟合

library(minpack.lm)
##Richard functions
richard.nls<-nlsLM(height ~a*(1-exp(-b*age))^c, data = ht.growth,
      start = list(a= 2, b= 0.1,c=0.5))
richard.nls
## Nonlinear regression model
##   model: height ~ a * (1 - exp(-b * age))^c
##    data: ht.growth
##        a        b        c 
## 34.93944  0.01023  1.72942 
##  residual sum-of-squares: 2.404
## 
## Number of iterations to convergence: 15 
## Achieved convergence tolerance: 1.49e-08
mpc(richard.nls,data=ht.growth)
##          MSE      RMSE   Rsquare       MAE       MAPE       RASE      AIC
## 1 0.08903215 0.2983826 0.9991126 0.2154854 0.03279708 0.01184323 19.31622
##        BIC
## 1 19.31622
##weibull functions
weibull.nls<-nlsLM(height ~a*(1-exp(-b*age^c)), data = ht.growth,
                   start = list(a= 200, b= 0.1,c=0.5))
weibull.nls
## Nonlinear regression model
##   model: height ~ a * (1 - exp(-b * age^c))
##    data: ht.growth
##         a         b         c 
## 3.289e+01 8.583e-04 1.446e+00 
##  residual sum-of-squares: 1.902
## 
## Number of iterations to convergence: 33 
## Achieved convergence tolerance: 1.49e-08
mpc(weibull.nls,data=ht.growth)
##          MSE      RMSE   Rsquare       MAE       MAPE       RASE      AIC
## 1 0.07043905 0.2654036 0.9992954 0.1939539 0.02006131 0.01065984 12.99148
##        BIC
## 1 12.99148
##logistic functions
logistic.nls<-nlsLM(height~a/(1+exp(b-c*age)),data = ht.growth,
                    start = list(a= 3, b= 10,c=1))
logistic.nls
## Nonlinear regression model
##   model: height ~ a/(1 + exp(b - c * age))
##    data: ht.growth
##        a        b        c 
## 29.94490  2.52614  0.02586 
##  residual sum-of-squares: 21.27
## 
## Number of iterations to convergence: 13 
## Achieved convergence tolerance: 1.49e-08
mpc(logistic.nls,data=ht.growth)
##         MSE      RMSE   Rsquare      MAE MAPE       RASE      AIC      BIC
## 1 0.7879212 0.8876493 0.9922644 0.722891  Inf 0.03973061 78.18704 78.18704
##korf functions
korf.nls<-nlsLM(height~a*exp(-b/(age^c)),data = ht.growth,
                start =list(a= 2, b= 10,c=0.5))
korf.nls
## Nonlinear regression model
##   model: height ~ a * exp(-b/(age^c))
##    data: ht.growth
##         a         b         c 
## 4.540e+00 3.890e-07 7.499e-01 
##  residual sum-of-squares: 7710
## 
## Number of iterations to convergence: 25 
## Achieved convergence tolerance: 1.49e-08
mpc(korf.nls,data=ht.growth)
##        MSE     RMSE    Rsquare     MAE      MAPE      RASE     AIC     BIC
## 1 285.5715 16.89886 0.05390254 14.3135 0.8682723 0.7866799 237.294 237.294
##mitscherlich functions
mitscherlich.nls<-nlsLM(height~a*(1-exp(-b*age)),data = ht.growth,
                start =list(a= 2, b= 0.1))
mitscherlich.nls
## Nonlinear regression model
##   model: height ~ a * (1 - exp(-b * age))
##    data: ht.growth
##         a         b 
## 52.821729  0.003597 
##  residual sum-of-squares: 26.33
## 
## Number of iterations to convergence: 15 
## Achieved convergence tolerance: 1.49e-08
mpc(mitscherlich.nls,data=ht.growth)
##         MSE      RMSE   Rsquare       MAE     MAPE       RASE      AIC      BIC
## 1 0.9750553 0.9874489 0.9906811 0.8913647 0.146115 0.04899004 81.94063 81.94063
##schumacher functions
schumacher.nls<-nlsLM(height~a*exp(-b/age),data = ht.growth,
                        start =list(a= 2, b= 0.1))
schumacher.nls
## Nonlinear regression model
##   model: height ~ a * exp(-b/age)
##    data: ht.growth
##     a     b 
## 44.54 97.77 
##  residual sum-of-squares: 16.86
## 
## Number of iterations to convergence: 9 
## Achieved convergence tolerance: 1.49e-08
mpc(schumacher.nls,data=ht.growth)
##         MSE      RMSE   Rsquare       MAE     MAPE       RASE      AIC      BIC
## 1 0.6243672 0.7901691 0.9941203 0.5265876 0.125357 0.02894163 69.90523 69.90523

##绘制非线性回归图

####regression figure####
plot(height~age,data=ht.growth,xlab="林木年龄 Age (year)",ylab="树高 Tree height (m)")
lines(age,predict(mod.lm),col=1)
lines(age,predict(richard.nls),col=2)
lines(age,predict(weibull.nls),col=3)
lines(age,predict(logistic.nls),col=4)
lines(age,predict(korf.nls),col=5)
lines(age,predict(mitscherlich.nls),col=6)
lines(age,predict(schumacher.nls),col=7)