R语言基于多元回归分析的大豆植株性状与产量的回归分析.doc

基于多元回归分析的大豆植株性状与产量的回归分析1. 背景与问题1.1背景选育高产优质的作物品种、采用高效的栽培技术，一直是农业科研工作者研究的主题。多年以来，人们为了解作物植株性状与产量的关系，以便做到更有效的性状选择和采取相应的栽培措施，进行了一系列的研究，关于作物产量与施肥、土壤条件、气候因子等的关系研究报道也较多。但影响大豆产量的植株自身因素还有许多种，有待于继续探究。1.2问题本文主要研究在作物植株诸多性状如生育日数x1、株高x2、有效分枝数x3、主茎节数x4、单株荚数x5、单株粒数x6、每荚粒数x7、百粒重x8、单株粒重x9，九种大豆植株性状与小区产量y之间的关系，从中找出产量的主导性状因子，为大豆高产育种中各农艺性状的选择和高产栽培技术提供理论依据。2. 程序设计与数据分析2.1程序设计本文将大豆植株的九种性状与小区产量的多元回归分析基础上剔除自相关性强的自变量、建立较高精度和稳定产量回归方程分析，本文利用多元回归分析，首先完成初步的回归分析；利用残差分析对回归模型的假设条件即随机误差项是否独立同分布进行检验；通过共线性诊断、逐步回归的方法，建立产量回归模型，进行产量分析，从而的出对大豆产量影响较大的植株性状，对各回归变量的作用大小作出评价。2.2分析方法2.2.1多元线性回归模型的建立多元线性回归分析是研究一个响应变量与多个自变量间呈线性相关关系的问题，这种关系可以用多元线性回归方程来描述：式中 为回归常数项， （i=1，2，k）称为偏回归系数，其意义为当其它自变量对响应变量的影响固定时，对应的第 i个自变量对 的线性影响程度。2.2.2残差分析残差是指由回归方程计算所得的预测值与实际样本值之间的差距，定义为，它是回归模型的估计值，由多个形成的序列称为残差序列，如果回归方程能够很好的反映被解释变量的特征和变化规律，那么残差序列中不应包含明显的规律性和趋势性。2.2.3多重共线性检验与修正逐步回归法逐步回归的基本思想是：对全部因子按其对影响程度大小（偏回归平方的大小），从大到小地依次逐个地引入回归方程，并随时对回归方程当时所含的全部变量进行检验，看其是否仍然显著，如不显著就将其剔除，知道回归方程中所含的所有变量对的作用都显著是，才考虑引入新的变量。再在剩下的未选因子中，选出对作用最大者，检验其显著性，显著着，引入方程，不显著，则不引入。直到最后再没有显著因子可以引入，也没有不显著的变量需要剔除为止。逐步回归分析时在考虑的全部自变量中按其对的贡献程度大小，由大到小地逐个引入回归方程，而对那些对作用不显著的变量可能是中不被引入回归方程。另外，已被引入回归方程的变量在引入新变量进行检验后失去重要性时，需要从回归方程中剔除出去。Step 1 计算变量均值和差平方和记各自的标准化变量为Step 2 计算的相关系数矩阵。Step 3 设已经选上了个变量：且互不相同，经过变换后为对逐一计算标准化变量的偏回归平方和，记，作检验，对给定的显著性水平，拒绝域为。Step 4 最Step 3 循环，直至最终选上了个变量，且互不相同，经过变换后为，则对应的回归方程为：，通过代数运算可得。2.3 数据来源及分析2.3.1数据来源本文数据采用2010年吉林省大豆种植研究数据中的274个大豆品种中的8个植株性状和生育日数及小区产量进行回归分析，其中植株性状选用：株高x2、有效分枝数x3、主茎节数x4、单株荚数x5、单株粒数x6、每荚粒数x7、百粒重x8、单株粒重x9，生育日数x1及小区产量y。2.3.2数据分析本文利用R语言对以上数据进行分析，分析过程及结果如下：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统计量估计值为24.44，对应p值2.2e-16比显著水平0.05小，说明回归方程是显著的。可决系数为0.4711，修正的可决系数为0.4518。2、 残差分析 左上图是拟合值与残差的散点图，从图上可以发现，所有点基本上是随机地分散在纵坐标值为-3 和+3的两条平行线之间，这说明随机误差项具有同方差性；左下图是拟合值与残差的标准差的散点图，其意义与上面类似；右上图表明随机误差项是服从正态分布的，其原因是正态Q-Q 图近似地可以看成一条直线；右下图的CooK 距离图进一步证实第6 个观测值是一个离群点，它对回归方程的影响是比较大的，要根据具体问题，讨论出现这一观测值的实际背景。3、多重共线性检验利用R语言计算解释变量相关系数矩阵的条件数k，k<100多重共线性程度很小，100<k<1000较强，>1000严重，计算结果为：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 * 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，可得出回归方程为：3、结果及结论由以上回归方程可得出：对大豆产量影响较大的因素从高到底分别为：株高x2>生育日数x1>单株荚数x5>有效分枝数x3>单株粒数x6>单株粒重x9。其中，影响最大的是大豆植株的株高，最小的是单株粒重，生育日数、单株粒数与单株粒重与产量呈正相关，而株高、有效分枝数、单株荚数与产量呈现负相关。说明如要提高大豆的产量，应选用株高较低、有效分枝数较少、单株荚数较少、生育日数较多、单株粒数较多、单株粒重较多的大豆植株。4、讨论对于大豆的育种，应该提高大豆种子的生育日数，延长生长时间，保证种子获取充足的养料，同时增加大豆植株的单株粒重和单株粒数，在植株数量不变的条件下，增加每株大豆植株的产量，同时应降低植株的株高，越高的植株产量越低，减少有效分枝数和单株荚数，对提高大豆产量都起着重要的作用。程序：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<-lm(yx1+x2+x3+x4+x5+x6+x7+x8+x9)summary(mul_re)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 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 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）$values1 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.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.675470526 -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.68738855 -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_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<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.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<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 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<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 1286258 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 = y x1 + x2 + x3 + x5 + x6 + x9)Residuals: Min 1Q Median 3Q Max -572.39 -130.77 -2.68 138.96 574.46 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 * 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 参考文献1韩秉进,潘相文,金剑,王光华,刘长江,刘晓冰 大豆植株性状相关性与产量回归分析J.20082李秀敏,蔡霞,使用统计软件R进行多元回归分析A3闫 昊,王 博,刘宝泉 大豆主茎节数、节间长度遗传分析及与株高关系研究A.2010