拓端tecdat：R語言因子實驗設計nlme擬合非線性混合模型分析有機農(nóng)業(yè)施氮水平
原文鏈接：http://tecdat.cn/?p=24134?
原文出處：拓端數(shù)據(jù)部落公眾號
測試非線性回歸中的交互作用
因子實驗在農(nóng)業(yè)中非常普遍，它們通常用于測試實驗因素之間相互作用的重要性。例如，可以在兩種不同的施氮水平（例如高和低）下進行基因型評估，以了解基因型的排名是否取決于養(yǎng)分的可用性。對于那些不太了解農(nóng)業(yè)的人，我只會說這樣的評估是相關的，因為我們需要知道我們是否可以推薦相同的基因型，例如，在傳統(tǒng)農(nóng)業(yè)（高氮可用性）和有機農(nóng)業(yè)中農(nóng)業(yè)氮的可用性。
讓我們考慮一個實驗，在該實驗中，我們在完整的區(qū)組因子設計中以兩種氮含量（“高”和“低”）測試了三種基因型（為了簡便起見，我們稱它們?yōu)?A、B 和 C），并進行四次重復。在八個不同的時間（播種后天數(shù)：DAS）從 24 個地塊中的每一個中取出生物量子樣本，以評估生物量隨時間的增長。
加載數(shù)據(jù)并將“Block”變量轉(zhuǎn)換為一個因子。
head(dataset)
數(shù)據(jù)集由以下變量組成：
'Id'：觀察的數(shù)字代碼
'DAS'：即播種后的天數(shù)。這是采集樣本的時刻
'Block', 'Plot', 'GEN' 和 'N' 分別代表每個觀察的塊、圖、基因型和氮水平
“產(chǎn)量”代表收獲的生物量。
查看觀察到的增長數(shù)據(jù)，如下圖所示。
我們看到增長是對稱的（大概是邏輯的）并且觀察的方差隨著時間的推移而增加，即方差與期望因變量成正比。
問題是：我們?nèi)绾畏治鲞@些數(shù)據(jù)？
模型
我們可以憑經(jīng)驗假設生物量和時間之間的關系是邏輯的：
其中Y是第i個基因型、第j個氮水平、第k個區(qū)塊和第l個小區(qū)在X時間觀察到的生物量產(chǎn)量，d是時間進入無窮大時的最大漸進生物量水平，b是拐點處的斜率，而e是生物量產(chǎn)量等于d/2時的時間。我們主要對參數(shù)d和e感興趣：第一個參數(shù)描述基因型的產(chǎn)量潛力，而第二個參數(shù)給出生長速度的測量。
?
每個小區(qū)都有重復測量，因此，模型參數(shù)可能顯示出一些變化，取決于基因型、氮水平、區(qū)塊和小區(qū)。特別是，假設b是相當恒定的，并且獨立于上述因素，而d和e可能根據(jù)以下公式發(fā)生變化，這是可以接受的。
其中，對于每個參數(shù)，μμ是截距，gg是第i個基因型的固定效應，NN是第j個氮水平的固定效應，gNgN是固定交互效應，θ是區(qū)塊的隨機效應，而ζζ是區(qū)塊內(nèi)地塊的隨機效應。這兩個方程完全等同于通常用于線性混合模型的方程，在雙因素因子區(qū)塊設計的情況下，其中ζζ是殘差誤差項。事實上，原則上，我們也可以考慮兩步法的擬合程序，即我們。
將邏輯模型擬合到每個圖的數(shù)據(jù)并獲得 d?和 e?的估計值
使用這些估計來擬合線性混合模型
我們不會在這里追求這種兩步法，我們將專注于一步擬合。
錯誤的方法
如果觀察是獨立的（即沒有塊和沒有重復測量），這個模型可以通過使用傳統(tǒng)的非線性回歸來擬合。
編碼報告如下。產(chǎn)量 "是(～)DAS的函數(shù)，通過一個三參數(shù)的Logistic函數(shù)。對于基因型和氮水平的不同組合必須擬合不同的曲線（id = N:GEN），盡管這些曲線應該部分地基于共同的參數(shù)值（'models = ...）。model"參數(shù)需要一些補充說明。它必須是一個矢量，其元素數(shù)與模型中的參數(shù)數(shù)一樣多（在本例中是三個：B、D和E）。每個元素代表一個變量的線性函數(shù)，并按字母順序指向參數(shù)，即第一個元素指b，第二個指d，第三個指e。參數(shù)b不依賴于任何變量（'~1'），因此在不同的曲線上擬合出一個常數(shù)；d和e依賴于基因型和氮水平的完全因子組合（~N*GEN = ~N + GEN + N:GEN）。最后，我們使用參數(shù)'bcVal = 0.5'來指定我們打算使用轉(zhuǎn)換兩邊方法，即對方程的兩邊進行對數(shù)轉(zhuǎn)換。這對于考慮異方差是必要的，但它不影響參數(shù)估計。
rm(Yield ~ DAS, dta =daas,id = GEN：N,model = c( ~ 1, ?~ N*GEN, ?~ N*GEN))
這個模型對于其他情況（無區(qū)塊和無重復測量）可能是有用的，但在我們的例子中是錯誤的。事實上，觀測值在區(qū)塊和地塊內(nèi)是聚在一起的；如果忽略這一點，我們就違反了模型殘差獨立的假設。殘差與擬合值的圖顯示，不存在異方差的問題。
?
考慮到上述情況，我們必須在這里使用不同的模型，盡管我將證明這種擬合可能會很有用。
非線性混合模型擬合
為了解釋觀察的類，我們切換到非線性混合效應模型（NLME）。一個不錯的選擇是'nlme()' 函數(shù)（Pinheiro 和 Bates，2000），盡管有時語法可能很麻煩。我們需要指定以下內(nèi)容：
模型參數(shù)的線性函數(shù)。nlme'函數(shù)中的'fixed'參數(shù)與上面函數(shù)中的'models'參數(shù)非常相似，即它需要一個列表，其中每個元素都是變量的線性函數(shù)。唯一的區(qū)別是，參數(shù)名稱需要在函數(shù)的左側(cè)指定。
模型參數(shù)的隨機效應。這些是通過使用 "隨機 "參數(shù)來指定的。在這種情況下，參數(shù)d和e有望在一個區(qū)塊內(nèi)的不同區(qū)塊和不同地塊之間顯示隨機變化。為了簡單起見，由于參數(shù)b不受基因型和氮水平的影響，我們也希望它在區(qū)塊和地塊之間不顯示任何隨機變化。
模型參數(shù)的起始值。我們需要指定模型參數(shù)的初始值。在這種情況下，我決定使用上面非線性回歸的輸出。
方程兩邊的變換。
nlme(sqtYld ~ srtLS.L3(DAS, b, d, e)tTable
從上圖中，我們看到整體擬合良好。隨機效應的固定效應和方差分量按如下方式獲得：
summary(mdnle1)
現(xiàn)在，讓我們回到我們最初的目的：測試 "基因型x氮 "交互作用的顯著性。事實上，我們有兩個可用的測試：一個是參數(shù)d，一個是參數(shù)e。首先，我們對兩個 "簡化 "模型進行編碼。為此，我們把固定效應從'~ N*GEN'改為'~ N + GEN'。同樣在這種情況下，我們使用非線性回歸擬合來獲得模型參數(shù)的起始值，用于下面的NLME模型擬合。
?
mdNave2 <- Yied ~ DASmodes = c( ~ 1, ?~ N + GEN, ?~ N * GENmdnle2 <-sqrt(Yeld) ~ sqrt(NS.3(DAS, b, d, emdaie3 <- Yied ~ DAScrid = N:GENmodls = c( ~ 1, ?~ N*GEN, ?~ N + GENmdne3 <- sqrt(Yild) ~ sqrt(NS.L3(DAS, b, d, eranom = d + e~ 1|Blck/Potfixd = lit(b ~ 1 d ~ N*GN, e ~ N + GN)
讓我們考慮第一個縮小的模型'modnlme2'。在這個模型中，"基因型x氮 "的交互作用已經(jīng)被移除，用于d參數(shù)。我們可以通過使用似然比檢驗來比較這個模型和完整的模型 "modnlme1"。
aova(mode1, mdne2)
該檢驗是顯著的，但兩個模型的AIC值非常接近。考慮到混合模型中的LRT通常比較寬松，應該可以得出結論，"基因型x氮素 "的交互作用不顯著，因此，用d參數(shù)衡量的基因型在產(chǎn)量潛力方面的排名應該與氮素水平有關。
現(xiàn)在讓我們考慮第二個簡化模型“modnlme3”。在第二個模型中，“e”參數(shù)的“基因型 x 氮”交互作用已被刪除。我們還可以使用似然比檢驗將此簡化模型與完整模型“modnlme1”進行比較：
anva(mole1, mdle3)
在這第二次檢驗中，"基因型x氮素 "交互作用的不顯著性似乎比第一次檢驗更可信。
最受歡迎的見解
1.基于R語言的lmer混合線性回歸模型
2.R語言用Rshiny探索lme4廣義線性混合模型（GLMM）和線性混合模型（LMM）
3.R語言線性混合效應模型實戰(zhàn)案例
4.R語言線性混合效應模型實戰(zhàn)案例2
5.R語言線性混合效應模型實戰(zhàn)案例
6.線性混合效應模型Linear Mixed-Effects Models的部分折疊Gibbs采樣
7.R語言LME4混合效應模型研究教師的受歡迎程度
8.R語言中基于混合數(shù)據(jù)抽樣(MIDAS)回歸的HAR-RV模型預測GDP增長
9.使用SAS，Stata，HLM，R，SPSS和Mplus的分層線性模型HLM
