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

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

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

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

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

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

6、1、 數(shù)據(jù)選用及處理整理已選用好的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 從輸出結(jié)果可以看出，修正決定系數(shù)為0.4518，剩余方差估計值=，f統(tǒng)計量估計值為24.44，對應(yīng)p值2.2e-16比顯著水平0.05小，說明回歸方程是顯著的?？蓻Q系數(shù)為0.

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

8、1000較強，>1000嚴重，計算結(jié)果為：186.93968,k<100多重共線性程度很小；4、修正多重共線性逐步回歸coefficients: estimate std. error t value pr(>|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

9、* 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.05 . 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: < 2.2e-16 p值遠小于0.05，可得出回歸方程

10、為：3、結(jié)果及結(jié)論由以上回歸方程可得出：對大豆產(chǎn)量影響較大的因素從高到底分別為：株高x2>生育日數(shù)x1>單株莢數(shù)x5>有效分枝數(shù)x3>單株粒數(shù)x6>單株粒重x9。其中，影響最大的是大豆植株的株高，最小的是單株粒重，生育日數(shù)、單株粒數(shù)與單株粒重與產(chǎn)量呈正相關(guān)，而株高、有效分枝數(shù)、單株莢數(shù)與產(chǎn)量呈現(xiàn)負相關(guān)。說明如要提高大豆的產(chǎn)量，應(yīng)選用株高較低、有效分枝數(shù)較少、單株莢數(shù)較少、生育日數(shù)較多、單株粒數(shù)較多、單株粒重較多的大豆植株。4、討論對于大豆的育種，應(yīng)該提高大豆種子的生育日數(shù)，延長生長時間，保證種子獲取充足的養(yǎng)料，同時增加大豆植株的單株粒重和單株粒數(shù)，在植株數(shù)量不變的

11、條件下，增加每株大豆植株的產(chǎn)量，同時應(yīng)降低植株的株高，越高的植株產(chǎn)量越低，減少有效分枝數(shù)和單株莢數(shù)，對提高大豆產(chǎn)量都起著重要的作用。程序：1、 導(dǎo)入數(shù)據(jù)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<-lm(yx1+x2+x3+x4+x5+x6+x7+x8+x9)summary(mul_re)

12、call:lm(formula = y x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9)residuals: min 1q median 3q max -585.88 -130.43 -9.72 141.78 595.01 coefficients: estimate std. error t value pr(>|t|) (intercept) -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

13、 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.761 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 de

14、grees 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）$values1 4.37328955 1.76790710 1.04413257 0.84241371 0.46836239 0.26071313 7 0.12949880 0.08029930 0.03338346 $

15、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.09604982 0.80962961 4, -0.420123249 0.114524926 -0.021857366 -0.37200572 -0.02930968 5, -0.420740277 -0.163

16、355262 -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.675470526 -0.003843284 0.10749279 0.41987990 9, -0.224070741 0.549213074 -0.006772412 0.47590055 -0.05983518 ,

17、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.68738855 -0.04674896 5, 0.06483988 -0.359950292 -0.08059018 -0.69291355 6, -0.05265000 -0.270465844 -0.06668541

18、 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_re)summary(mul_step)start: aic=2768.53y x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 df sum of sq rss aic

19、- x4 1 1187 11340059 2766.6- x7 1 25837 11364710 2767.1- x2 1 59986 11398859 2767.9- x8 1 60056 11398928 2767.9<none> 11338873 2768.5- x1 1 145242 11484114 2769.8- x5 1 209190 11548063 2771.2- x3 1 469103 11807976 2776.9- x6 1 787576 12126448 2783.8- x9 1 828618 12167491 2784.7step: aic=2766.5

20、6y 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<none> 11340059 2766.6- x2 1 108567 11448626 2767.0- x1 1 159779 11499838 2768.2- x5 1 210458 11550518 2769.3- x3 1 472766 11812825 2775.1- x6 1 808970 12149029 2782.3- x9 1 827

21、551 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<none> 11366492 2765.2- x2 1 109547 11476039 2765.6- x1 1 157130 11523622 2766.7- x5 1 322295 11688787 2770.3- x3 1 468587 11835079 2773.5- x9 1 801868 12168360 2780.7- x6 1 1

22、286258 12652750 2790.7step: aic=2765.07y x1 + x2 + x3 + x5 + x6 + x9 df sum of sq rss aic<none> 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