用R語(yǔ)言擬合籽粒生長(zhǎng)模型

1、用R語(yǔ)言擬合籽粒生長(zhǎng)模型簡(jiǎn)介：生物的生長(zhǎng)特征在多少情況下可以用s型曲線進(jìn)行描述，即：起始階段生長(zhǎng)慢，然后 進(jìn)入一個(gè)快速增長(zhǎng)期，最后生長(zhǎng)速度再次變慢并逐漸接近生長(zhǎng)極限。谷物籽粒的發(fā)育過(guò)程也 可以用S型曲線進(jìn)行描述，比如玉米、水稻和小麥的籽粒生長(zhǎng)特征一般用logistic方程進(jìn)行 描述。logistic方程的特點(diǎn)是其曲線的拐點(diǎn)就是其中點(diǎn)，這樣有一個(gè)好處，就是能夠簡(jiǎn)化數(shù) 據(jù)擬合過(guò)程。但是，在生長(zhǎng)過(guò)程中遭受脅迫時(shí)，生長(zhǎng)曲線的拐點(diǎn)就不是曲線的中點(diǎn)，此時(shí)不 再適用logistic方程。除了 logistic方程之外，還有兩個(gè)常用方程即richards方程和gompertz 方程。本文用小麥籽粒的生長(zhǎng)數(shù)據(jù)講

2、解如何用R語(yǔ)言進(jìn)行非線性模型的擬合，并比較3個(gè) 生長(zhǎng)曲線的優(yōu)劣。籽粒生長(zhǎng)過(guò)程主要是籽粒內(nèi)含物的填充過(guò)程，因此也稱為籽粒灌漿。在作物開(kāi)花以后， 籽粒開(kāi)始形成粒重逐漸增加。粒重與時(shí)間的關(guān)系不是線性關(guān)系，因此，用數(shù)學(xué)公式描述粒重 隨時(shí)間的變化動(dòng)態(tài)也稱為非線性模型擬合。擬合的目的在于用嚴(yán)謹(jǐn)?shù)臄?shù)學(xué)公式描述粒重?cái)?shù)據(jù) 的變化規(guī)律，并給出公式的描述誤差。有了這個(gè)公式，就可以在籽粒生長(zhǎng)結(jié)束之前預(yù)測(cè)籽粒 的最終重量。觀測(cè)、擬合與預(yù)測(cè)是研究過(guò)程的3個(gè)階段。非線性模型分為簡(jiǎn)單非線性模型、一般非線性模型和非線性混合模型。非線性混合模型 用群體效應(yīng)、群體和個(gè)體的差異描述一組試驗(yàn)數(shù)據(jù)，結(jié)構(gòu)更加緊湊，模型參數(shù)更少，擬合過(guò) 程

3、也更加復(fù)雜。所以，本文先介紹簡(jiǎn)單的非線性模型為非線性混合模型奠定基礎(chǔ)。關(guān)鍵詞：R語(yǔ)言，籽粒生長(zhǎng)，籽粒灌漿，非線性模型，模型擬合，生長(zhǎng)模擬以小麥籽粒生長(zhǎng)為例，在花后測(cè)量籽粒的干重，以粒重為因變量y），以時(shí)間為自變量 （t），研究粒重與時(shí)間的關(guān)系。數(shù)據(jù)如下: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為開(kāi)花后天數(shù)，grainWeight為籽粒干重，單位是mg。3個(gè)方程分別為:logis

6、tic = function(t, k, a, b)k/(1 + exp(a - b*t)其中 K 為最大粒質(zhì)量,t 為花后時(shí)間，a、 b為待定系數(shù)。gompertz = function(t, k, a, b)k*exp(-exp(a - b*t)其中，t 為自變量，y 為因變量；k、 a、b是未知數(shù)(k，a0，b手1。richards = function(t, k, a, b, c)k*(1 + exp(a-b*t)A(1/(1-c) richards 方程是上述兩個(gè)方 程的更一般的形式，該方程通過(guò)增加參數(shù)c來(lái)調(diào)整曲線的性狀，使擬合值與觀測(cè)值的殘差更 小。1、在R中定義上述3個(gè)方程，運(yùn)行

7、下列代碼：#定義3參數(shù)logistic函數(shù)logistic = function(t, k, a, b)k/(1 + exp(a - b*t)#定義3參數(shù)gompertz函數(shù)gompertz = function(t, k, a, b)k*exp(-exp(a - b*t)#定義函數(shù)richardsrichards = function(t, k, a, b, c)k*(1 + exp(a-b*t)A(1/(1-c)2、做散點(diǎn)圖，查看數(shù)據(jù)趨勢(shì)，運(yùn)行下列代碼：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#從散點(diǎn)圖中可以看出粒重的方差隨花后天數(shù)而增大#因此，擬合異方差模型3、擬合模型并比較3個(gè)模型的優(yōu)劣，運(yùn)行下列代碼：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)#比較模型優(yōu)劣anova(logis_mdl, gomp_mdl)#logistic 優(yōu)于 gompanova(logis_mdl, rich_mdl)#p = 0.03,rich 更好anova(gomp_mdl, rich_mdl)#p 0.0001,rich 更好#異方差richards模型最好下圖是3條曲線的比較：紅色表示richards模型，藍(lán)色表示logistic模型，黑色表示gompertz模型。從圖中可以 直觀地看出紅色曲線最優(yōu)。4、計(jì)算模型參數(shù)以及二級(jí)參數(shù)（灌漿持續(xù)期、最大速率）summary(rich_mdl)#44.47015#求解高階導(dǎo)數(shù)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)#定義最大速率函數(shù)，model=gnls()的返回模型，t0=求解最大速率出現(xiàn)時(shí)間的初始值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)#二 階導(dǎo)數(shù)值dy - deriv(DD(expr, t, 3), t, func = TRUE)#三階導(dǎo)數(shù)值t_max - newton(y, dy, t0)rate - deriv(DD(expr, t, 1), t, func = TRUE)r - rate(t_max)# 一 階導(dǎo)數(shù)值cat(maximum rate =)return(r)#定義持續(xù)期函數(shù)，model=gnls()的返回模型，d0=求解持續(xù)期的初始值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#所需的參數(shù)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