R語言基于多元回歸分析的大豆植株性狀與產量的回歸分析

1、基于多元回歸分析的大豆植株性狀與產量的回歸分析1. 背景與問題1.1背景選育高產優(yōu)質的作物品種、采用高效的栽培技術，一直是農業(yè)科研工作者研究的主題。多年以來，人們?yōu)榱私庾魑镏仓晷誀钆c產量的關系，以便做到更有效的性狀選擇和采取相應的栽培措施，進行了一系列的研究，關于作物產量與施肥、土壤條件、氣候因子等的關系研究報道也較多。但影響大豆產量的植株自身因素還有許多種，有待于繼續(xù)探究。1.2問題本文主要研究在作物植株諸多性狀如生育日數x1、株高x2、有效分枝數x3、主莖節(jié)數x4、單株莢數x5、單株粒數x6、每莢粒數x7、百粒重x8、單株粒重x9，九種大豆植株性狀與小區(qū)產量y之間的關系，從中找出產量的主導

2、性狀因子，為大豆高產育種中各農藝性狀的選擇和高產栽培技術提供理論依據。2. 程序設計與數據分析2.1程序設計本文將大豆植株的九種性狀與小區(qū)產量的多元回歸分析基礎上剔除自相關性強的自變量、建立較高精度和穩(wěn)定產量回歸方程分析，本文利用多元回歸分析，首先完成初步的回歸分析；利用殘差分析對回歸模型的假設條件即隨機誤差項是否獨立同分布進行檢驗；通過共線性診斷、逐步回歸的方法，建立產量回歸模型，進行產量分析，從而的出對大豆產量影響較大的植株性狀，對各回歸變量的作用大小作出評價。2.2分析方法2.2.1多元線性回歸模型的建立多元線性回歸分析是研究一個響應變量與多個自變量間呈線性相關關系的問題，這種關系可以用

3、多元線性回歸方程來描述：式中 為回歸常數項， （i=1，2，k）稱為偏回歸系數，其意義為當其它自變量對響應變量的影響固定時，對應的第 i個自變量對 的線性影響程度。2.2.2殘差分析殘差是指由回歸方程計算所得的預測值與實際樣本值之間的差距，定義為，它是回歸模型的估計值，由多個形成的序列稱為殘差序列，如果回歸方程能夠很好的反映被解釋變量的特征和變化規(guī)律，那么殘差序列中不應包含明顯的規(guī)律性和趨勢性。2.2.3多重共線性檢驗與修正逐步回歸法逐步回歸的基本思想是：對全部因子按其對影響程度大?。ㄆ貧w平方的大?。?，從大到小地依次逐個地引入回歸方程，并隨時對回歸方程當時所含的全部變量進行檢驗，看其是否仍然

4、顯著，如不顯著就將其剔除，知道回歸方程中所含的所有變量對的作用都顯著是，才考慮引入新的變量。再在剩下的未選因子中，選出對作用最大者，檢驗其顯著性，顯著著，引入方程，不顯著，則不引入。直到最后再沒有顯著因子可以引入，也沒有不顯著的變量需要剔除為止。逐步回歸分析時在考慮的全部自變量中按其對的貢獻程度大小，由大到小地逐個引入回歸方程，而對那些對作用不顯著的變量可能是中不被引入回歸方程。另外，已被引入回歸方程的變量在引入新變量進行檢驗后失去重要性時，需要從回歸方程中剔除出去。step 1 計算變量均值和差平方和記各自的標準化變量為step 2 計算的相關系數矩陣。step 3 設已經選上了個變量：且互

5、不相同，經過變換后為對逐一計算標準化變量的偏回歸平方和，記，作檢驗，對給定的顯著性水平，拒絕域為。step 4 最step 3 循環(huán)，直至最終選上了個變量，且互不相同，經過變換后為，則對應的回歸方程為：，通過代數運算可得。2.3 數據來源及分析2.3.1數據來源本文數據采用2010年吉林省大豆種植研究數據中的274個大豆品種中的8個植株性狀和生育日數及小區(qū)產量進行回歸分析，其中植株性狀選用：株高x2、有效分枝數x3、主莖節(jié)數x4、單株莢數x5、單株粒數x6、每莢粒數x7、百粒重x8、單株粒重x9，生育日數x1及小區(qū)產量y。2.3.2數據分析本文利用r語言對以上數據進行分析，分析過程及結果如下：

6、1、 數據選用及處理整理已選用好的9個自變量及一個因變量，剔除缺失值，進行線性回歸分析：residual standard error: 214.3 on 247 degrees of freedommultiple r-squared: 0.4711, adjusted r-squared: 0.4518 f-statistic: 24.44 on 9 and 247 df, p-value: 2.2e-16 從輸出結果可以看出，修正決定系數為0.4518，剩余方差估計值=，f統(tǒng)計量估計值為24.44，對應p值2.2e-16比顯著水平0.05小，說明回歸方程是顯著的?？蓻Q系數為0.4711，

7、修正的可決系數為0.4518。2、 殘差分析 左上圖是擬合值與殘差的散點圖，從圖上可以發(fā)現(xiàn)，所有點基本上是隨機地分散在縱坐標值為-3 和+3的兩條平行線之間，這說明隨機誤差項具有同方差性；左下圖是擬合值與殘差的標準差的散點圖，其意義與上面類似；右上圖表明隨機誤差項是服從正態(tài)分布的，其原因是正態(tài)q-q 圖近似地可以看成一條直線；右下圖的cook 距離圖進一步證實第6 個觀測值是一個離群點，它對回歸方程的影響是比較大的，要根據具體問題，討論出現(xiàn)這一觀測值的實際背景。3、多重共線性檢驗利用r語言計算解釋變量相關系數矩陣的條件數k，k100多重共線性程度很小，100k1000嚴重，計算結果為：186.

8、93968,k|t|) (intercept) 231.4579 88.2335 2.623 0.00925 * x1 2.0592 1.0755 1.915 0.05667 . x2 -1.2598 0.8086 -1.558 0.12051 x3 -63.9781 20.1098 -3.181 0.00165 * x5 -9.2091 3.2953 -2.795 0.00560 * x6 6.7922 1.3255 5.124 5.98e-07 *x9 32.1232 3.9792 8.073 2.90e-14 *-signif. codes: 0 * 0.001 * 0.01 * 0.0

9、5 . 0.1 1 residual standard error: 214 on 250 degrees of freedommultiple r-squared: 0.4658, adjusted r-squared: 0.453 f-statistic: 36.33 on 6 and 250 df, p-value: 生育日數x1單株莢數x5有效分枝數x3單株粒數x6單株粒重x9。其中，影響最大的是大豆植株的株高，最小的是單株粒重，生育日數、單株粒數與單株粒重與產量呈正相關，而株高、有效分枝數、單株莢數與產量呈現(xiàn)負相關。說明如要提高大豆的產量，應選用株高較低、有效分枝數較少、單株莢數較少

10、、生育日數較多、單株粒數較多、單株粒重較多的大豆植株。4、討論對于大豆的育種，應該提高大豆種子的生育日數，延長生長時間，保證種子獲取充足的養(yǎng)料，同時增加大豆植株的單株粒重和單株粒數，在植株數量不變的條件下，增加每株大豆植株的產量，同時應降低植株的株高，越高的植株產量越低，減少有效分枝數和單株莢數，對提高大豆產量都起著重要的作用。程序：1、 導入數據m-read.table(jilin1.txt,header=true)y-m,10;x1-m,1;x2-m,2;x3-m,3;x4-m,4;x5-m,5;x6-m,6;x7-m,7;x8-m,8;x9-m,92、 做回歸mul_re|t|) (in

11、tercept) -64.408 233.272 -0.276 0.78270 x1 1.970 1.108 1.779 0.07651 . x2 -1.391 1.217 -1.143 0.25410 x3 -64.458 20.164 -3.197 0.00157 * x4 1.639 10.195 0.161 0.87241 x5 -7.692 3.603 -2.135 0.03377 * x6 6.688 1.615 4.142 4.73e-05 *x7 58.408 77.856 0.750 0.45384 x8 10.079 8.812 1.144 0.25382 x9 26.76

12、1 6.299 4.249 3.05e-05 *-signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1 residual standard error: 214.3 on 247 degrees of freedommultiple r-squared: 0.4711, adjusted r-squared: 0.4518 f-statistic: 24.44 on 9 and 247 df, p-value: 2.2e-16 3、回歸診斷：殘差圖plot(mul_re,which=1:4)4、多重共線性診斷x=cor(m1:9)kappa（x）$val

13、ues1 4.37328955 1.76790710 1.04413257 0.84241371 0.46836239 0.26071313 7 0.12949880 0.08029930 0.03338346 $vectors ,1 ,2 ,3 ,4 ,5 1, -0.356736427 0.284860137 0.081557776 -0.36565897 -0.23125149 2, -0.419532972 0.085193652 0.039641097 -0.39838195 0.08393449 3, -0.331105336 -0.301417546 0.159383136 0.

14、09604982 0.80962961 4, -0.420123249 0.114524926 -0.021857366 -0.37200572 -0.02930968 5, -0.420740277 -0.163355262 -0.009297722 0.40645079 -0.10301261 6, -0.412505231 -0.151811974 -0.182899521 0.38224782 -0.25492800 7, -0.009195756 -0.000999126 -0.965555211 -0.10305260 0.16625081 8, 0.114550048 0.675

15、470526 -0.003843284 0.10749279 0.41987990 9, -0.224070741 0.549213074 -0.006772412 0.47590055 -0.05983518 ,6 ,7 ,8 ,9 1, 0.77174628 -0.009563462 0.01131618 0.04412015 2, -0.38018274 0.110268675 -0.70154728 -0.02329927 3, 0.22305831 0.177712313 0.11576309 0.12103842 4, -0.42801145 -0.113018364 0.6873

16、8855 -0.04674896 5, 0.06483988 -0.359950292 -0.08059018 -0.69291355 6, -0.05265000 -0.270465844 -0.06668541 0.69416706 7, 0.10617363 0.088150347 -0.01110029 -0.10095476 8, -0.02962256 -0.574897504 -0.07434453 0.07640361 9, -0.08828021 0.634350797 0.07216092 -0.051461915、多重共線性修正逐步回歸mul_step-step(mul_

17、re)summary(mul_step)start: aic=2768.53y x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 df sum of sq rss aic- x4 1 1187 11340059 2766.6- x7 1 25837 11364710 2767.1- x2 1 59986 11398859 2767.9- x8 1 60056 11398928 2767.9 11338873 2768.5- x1 1 145242 11484114 2769.8- x5 1 209190 11548063 2771.2- x3 1 46910

18、3 11807976 2776.9- x6 1 787576 12126448 2783.8- x9 1 828618 12167491 2784.7step: aic=2766.56y x1 + x2 + x3 + x5 + x6 + x7 + x8 + x9 df sum of sq rss aic- x7 1 26433 11366492 2765.2- x8 1 61003 11401062 2765.9 11340059 2766.6- x2 1 108567 11448626 2767.0- x1 1 159779 11499838 2768.2- x5 1 210458 1155

19、0518 2769.3- x3 1 472766 11812825 2775.1- x6 1 808970 12149029 2782.3- x9 1 827551 12167610 2782.7step: aic=2765.16y x1 + x2 + x3 + x5 + x6 + x8 + x9 df sum of sq rss aic- x8 1 84791 11451283 2765.1 11366492 2765.2- x2 1 109547 11476039 2765.6- x1 1 157130 11523622 2766.7- x5 1 322295 11688787 2770.

20、3- x3 1 468587 11835079 2773.5- x9 1 801868 12168360 2780.7- x6 1 1286258 12652750 2790.7step: aic=2765.07y x1 + x2 + x3 + x5 + x6 + x9 df sum of sq rss aic 11451283 2765.1- x2 1 111176 11562460 2765.6- x1 1 167933 11619216 2766.8- x5 1 357727 11809011 2771.0- x3 1 463617 11914901 2773.3- x6 1 1202684 12653968 2788.7- x9 1 2985034 14436318 2822.6call:lm(formula = y x1 + x2 + x3 + x5 + x6 + x9)residuals: min