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

用R语言拟合籽粒生长模型简介：生物的生长特征在多少情况下可以用s型曲线进行描述，即：起始阶段生长慢，然后 进入一个快速增长期，最后生长速度再次变慢并逐渐接近生长极限谷物籽粒的发育过程也 可以用S型曲线进行描述，比如玉米、水稻和小麦的籽粒生长特征一般用logistic方程进行 描述logistic方程的特点是其曲线的拐点就是其中点，这样有一个好处，就是能够简化数 据拟合过程但是，在生长过程中遭受胁迫时，生长曲线的拐点就不是曲线的中点，此时不 再适用logistic方程除了 logistic方程之外，还有两个常用方程即richards方程和gompertz 方程本文用小麦籽粒的生长数据讲解如何用R语言进行非线性模型的拟合，并比较3个 生长曲线的优劣籽粒生长过程主要是籽粒内含物的填充过程，因此也称为籽粒灌浆在作物开花以后， 籽粒开始形成粒重逐渐增加粒重与时间的关系不是线性关系，因此，用数学公式描述粒重 随时间的变化动态也称为非线性模型拟合拟合的目的在于用严谨的数学公式描述粒重数据 的变化规律，并给出公式的描述误差有了这个公式，就可以在籽粒生长结束之前预测籽粒 的最终重量观测、拟合与预测是研究过程的3个阶段。

非线性模型分为简单非线性模型、一般非线性模型和非线性混合模型非线性混合模型 用群体效应、群体和个体的差异描述一组试验数据，结构更加紧凑，模型参数更少，拟合过 程也更加复杂所以，本文先介绍简单的非线性模型为非线性混合模型奠定基础关键词：R语言，籽粒生长，籽粒灌浆，非线性模型，模型拟合，生长模拟以小麦籽粒生长为例，在花后测量籽粒的干重，以粒重为因变量y），以时间为自变量 （t），研究粒重与时间的关系数据如下:daagrainWeightdaagrainWeightdaagrainWeightdaagrainWeightdaagrainWeight53.0885.49119.381413.431722.0652.8486.29119.271413.311721.7752.7685.56119.731412.991722.8652.9285.60119.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.923344.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个方程分别为:logistic = 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，a>0，b手1richards = function(t, k, a, b, c){k*(1 + exp(a-b*t))A(1/(1-c))} richards 方程是上述两个方 程的更一般的形式，该方程通过增加参数c来调整曲线的性状，使拟合值与观测值的残差更 小1、 在R中定义上述3个方程，运行下列代码：#定义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 = 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", fill = "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()40-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(daa,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(rich_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#求解高阶导数DD <- function(expr, name, order = 1) (if(order < 1) stop("'order' must be >= 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 <- 0while ((abs(fx) > tol) && (iter < max.iter)) (x <- x - fx/dfxiter <- iter + 1#cat("At iteration", iter, "value of x is :", x, "\n")}if (abs(fx) > tol) (cat("Algorithm failed to converged\n")return(NULL)} else(cat("Algorithm converged\n")return(x)}}#定义最大速率函数，model=gnls()的返回模型，t0=求解最大速率出现时间的初始值r_max <- function(model, t0)(k0 <- coefficients(model)[1]a0 <- coefficients(model)[2]b0 <- coefficients(model)[3]c0 <- coefficients(model)[4]expr <- substit。