用R语言拟合籽粒生长模型

1、用R语言拟合籽粒生长模型简介：生物的生长特征在多少情况下可以用s型曲线进行描述，即：起始阶段生长慢，然后 进入一个快速增长期，最后生长速度再次变慢并逐渐接近生长极限。谷物籽粒的发育过程也 可以用S型曲线进行描述，比如玉米、水稻和小麦的籽粒生长特征一般用logistic方程进行 描述。logistic方程的特点是其曲线的拐点就是其中点，这样有一个好处，就是能够简化数 据拟合过程。但是，在生长过程中遭受胁迫时，生长曲线的拐点就不是曲线的中点，此时不 再适用logistic方程。除了 logistic方程之外，还有两个常用方程即richards方程和gompertz 方程。本文用小麦籽粒的生长数据讲

2、解如何用R语言进行非线性模型的拟合，并比较3个 生长曲线的优劣。籽粒生长过程主要是籽粒内含物的填充过程，因此也称为籽粒灌浆。在作物开花以后， 籽粒开始形成粒重逐渐增加。粒重与时间的关系不是线性关系，因此，用数学公式描述粒重 随时间的变化动态也称为非线性模型拟合。拟合的目的在于用严谨的数学公式描述粒重数据 的变化规律，并给出公式的描述误差。有了这个公式，就可以在籽粒生长结束之前预测籽粒 的最终重量。观测、拟合与预测是研究过程的3个阶段。非线性模型分为简单非线性模型、一般非线性模型和非线性混合模型。非线性混合模型 用群体效应、群体和个体的差异描述一组试验数据，结构更加紧凑，模型参数更少，拟合过 程

3、也更加复杂。所以，本文先介绍简单的非线性模型为非线性混合模型奠定基础。关键词：R语言，籽粒生长，籽粒灌浆，非线性模型，模型拟合，生长模拟以小麦籽粒生长为例，在花后测量籽粒的干重，以粒重为因变量y），以时间为自变量 （t），研究粒重与时间的关系。数据如下:daagrainWeightdaagrainWeightdaagrainWeightdaagrainWeightdaagrainWeight53.0885.49119.381413.431722.0652.8486.29119.271413.311721.7752.7685.56119.731412.991722.8652.9285.60119

4、.111413.981721.1252.5885.78118.951412.681720.5152.5285.18118.051412.671718.5352.8785.50119.101413.161722.2652.8584.57117.931412.571721.0553.0785.58118.561413.791720.31daagrainWeightdaagrainWeightdaagrainWeightdaagrainWeightdaagrainWeight2028.652333.862739.102939.053342.642028.842334.862741.632941.92

5、3344.492032.322337.012740.622941.393344.912029.522336.232740.722939.863346.402028.512334.692739.872944.693345.722027.612334.282740.992941.603345.552029.662333.672739.782939.453344.762027.832333.642740.212941.223341.122031.292334.362744.182941.753341.07表中daa为开花后天数，grainWeight为籽粒干重，单位是mg。3个方程分别为:logis

6、tic = function(t, k, a, b)k/(1 + exp(a - b*t)其中 K 为最大粒质量,t 为花后时间，a、 b为待定系数。gompertz = function(t, k, a, b)k*exp(-exp(a - b*t)其中，t 为自变量，y 为因变量；k、 a、b是未知数(k，a0，b手1。richards = function(t, k, a, b, c)k*(1 + exp(a-b*t)A(1/(1-c) richards 方程是上述两个方 程的更一般的形式，该方程通过增加参数c来调整曲线的性状，使拟合值与观测值的残差更 小。1、在R中定义上述3个方程，运行

7、下列代码：#定义3参数logistic函数logistic = function(t, k, a, b)k/(1 + exp(a - b*t)#定义3参数gompertz函数gompertz = function(t, k, a, b)k*exp(-exp(a - b*t)#定义函数richardsrichards = function(t, k, a, b, c)k*(1 + exp(a-b*t)A(1/(1-c)2、做散点图，查看数据趋势，运行下列代码：library(ggplot2)sk_theme - theme(panel.background = element_rect(fill

8、 = NA, colour = black)+theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank()+theme(axis.text = element_text(face = italic, size = 8)+theme(strip.text.x = element_text(size = 10)+theme(axis.title = element_text(size = 10)+theme(strip.background = element_rect(colour = black, fil

9、l = white)+ theme(legend.margin = margin(6, 6, 6, 6),legend.text = element_text(size = 8)p1 - ggplot(data = gyt_data, aes(daa, grainWeight) +geom_point()+labs(x = DAA, y = expression(GM*(*frac(mg, grain)*)#作图png(filename = figl.png, width = 8, height = 6, units = cm, res = 600)p1+ sk_themedev.off()4

10、0-II30-2010-3020DAA#从散点图中可以看出粒重的方差随花后天数而增大#因此，拟合异方差模型3、拟合模型并比较3个模型的优劣，运行下列代码：library(nlme)logis_mdl - gnls(grainWeight logistic(daa,k,a,b),data = gyt_data,start = list(k = 45, a = 3.9, b = 0.2),weights = varPower()lines(gyt_data$daa, fitted(logis_mdl), col=blue)gomp_mdl - gnls(grainWeight gompertz(d

11、aa,k,a,b),data = gyt_data,start = list(k = 45, a = 3.9, b = 0.2),weights = varPower()lines(gyt_data$daa, fitted(gomp_mdl), col=black)rich_mdl - gnls(grainWeight richards(daa,k,a,b,c),data = gyt_data,start = list(k = 44.5, a = 2.5, b = 0.1, c = 1.88),weights = varPower()lines(gyt_data$daa, fitted(ric

12、h_mdl),col=red)#比较模型优劣anova(logis_mdl, gomp_mdl)#logistic 优于 gompanova(logis_mdl, rich_mdl)#p = 0.03,rich 更好anova(gomp_mdl, rich_mdl)#p 0.0001,rich 更好#异方差richards模型最好下图是3条曲线的比较：红色表示richards模型，蓝色表示logistic模型，黑色表示gompertz模型。从图中可以 直观地看出红色曲线最优。4、计算模型参数以及二级参数（灌浆持续期、最大速率）summary(rich_mdl)#44.47015#求解高阶导数D

13、D - function(expr, name, order = 1) (if(order = 1)if(order = 1) D(expr, name)else DD(D(expr, name), name, order - 1)#递归#牛顿迭代方法求方程的根newton - function(ftn, dftn, x0, tol = 1e-9, max.iter = 100)(x - x0fx - ftn(x)dfx - dftn(x)iter tol) & (iter max.iter) (x - x - fx/dfxiter tol) (cat(Algorithm failed to

14、convergedn)return(NULL) else(cat(Algorithm convergedn)return(x)#定义最大速率函数，model=gnls()的返回模型，t0=求解最大速率出现时间的初始值r_max - function(model, t0)(k0 - coefficients(model)1a0 - coefficients(model)2b0 - coefficients(model)3c0 - coefficients(model)4expr - substitute( k * (1 + exp(a - b * t) A (1 / (1 - c), list(

15、k=k0,a=a0,b=b0,c=c0)y - deriv(DD(expr, t, 2), t, func = TRUE)#二 阶导数值dy - deriv(DD(expr, t, 3), t, func = TRUE)#三阶导数值t_max - newton(y, dy, t0)rate - deriv(DD(expr, t, 1), t, func = TRUE)r - rate(t_max)# 一 阶导数值cat(maximum rate =)return(r)#定义持续期函数，model=gnls()的返回模型，d0=求解持续期的初始值a0 - coefficients(model)2

16、b0 - coefficients(model)3c0 - coefficients(model)4expr - substitute(1 + exp(a - b * t) A (1 / (1 - c) - 19/20,list(a=a0,b=b0,c=c0)f - deriv(expr, t, func = TRUE)df - deriv(DD(expr, t, 1), t, func = TRUE)dur - newton(f, df, d0)cat(Duration =)return(dur)#plot curverich_plot - function(model, start = 0

17、, end = 35, points = 101)(expr - function(t, k, a, b, c)(k - coefficients(model)1a - coefficients(model)2b - coefficients(model)3c - coefficients(model)4k*(1 + exp(a-b*t)A(1/(1-c)curve(expr, from = start, to = end, n = points,type = l, xlab = DAA, ylab = mg/grain)#拟合的模型rich_mdl#所需的参数coefficients(rich_mdl) #44.4701549 5.0014415 0.2588514 2.3851854summary(rich_mdl) #power = 0.6409366r_max(rich_mdl, 18) #2.576643duration(rich_mdl, 25) #29.39937rich_plot(rich_mdl, start = 0, end = 3