基本统计-(3)

1、第第3 3讲讲 常见统计分析常见统计分析内容提要内容提要 描述统计描述统计 频数表分析频数表分析 方差分析方差分析 t检验检验 卡方检验卡方检验 线性回归线性回归 相关分析相关分析3.1 描述分析描述分析(Descriptive statistics)描述统计就是把数据集所包含的信息加以简要地概况，描述统计就是把数据集所包含的信息加以简要地概况，如计算数据的数字特征、制作频数表和频数图等等，如计算数据的数字特征、制作频数表和频数图等等，用所获得的统计量和图表来描述数据集所反映的特用所获得的统计量和图表来描述数据集所反映的特征和规律，使得研究的问题更加简单、直观。征和规律，使得研究的问题更加简单

2、、直观。 描述性统计主要包括反映数据集中趋势的特征值描述性统计主要包括反映数据集中趋势的特征值( (比如比如平均数、中位数、众数、分位数平均数、中位数、众数、分位数) )、数据离散程度的、数据离散程度的特征值特征值( (比如方差、标准差、值域、变异系数比如方差、标准差、值域、变异系数) )和数和数据分布形态的特征值据分布形态的特征值( (比如偏度、峰度比如偏度、峰度) )。 3.1.1使用使用summary()函数函数oats0 比正态分布的高峰更加陡峭尖顶 Kurtosis0 正偏差数值较大，为正偏或右偏。长尾巴拖在右边 Skewness describe(oats$yield) descr

3、ibe(oats$yield) var n mean sd median trimmed mad min max range skew kurtosis se var n mean sd median trimmed mad min max range skew kurtosis se1 1 72 104 27 102 103 27 53 174 121 0.26 -0.46 3.21 1 72 104 27 102 103 27 53 174 121 0.26 -0.46 3.23.2 频数表频数表(Frequency table)分析分析 频数表分析是对数据集按数据范围分成若干区间，频数表

4、分析是对数据集按数据范围分成若干区间，即分成若干组，求出每组组中值，各组数据用组中值代即分成若干组，求出每组组中值，各组数据用组中值代替，计算各组数据的频数，并作出频数表。替，计算各组数据的频数，并作出频数表。 频数表分析例子频数表分析例子summary(oats$yield)# 计算频数计算频数A round(prop.table(A) round(prop.table(A) * * 100,2) # 100,2) # 计算频数比例计算频数比例 (40,60 (60,80 (80,100 (100,120 (120,140 (140,160 (160,180 (40,60 (60,80 (8

5、0,100 (100,120 (120,140 (140,160 (160,180 2.8 16.7 29.2 23.6 18.1 6.9 2.8 2.8 16.7 29.2 23.6 18.1 6.9 2.8Frequency chart of yieldyeildFrequency4060801001201401601800510 15 203.3 方差分析方差分析ANOVA方差分析是一种在若干组能相互比较的试验数据中，把方差分析是一种在若干组能相互比较的试验数据中，把产生变异的原因加以区分的方法与技术，其主要用途是产生变异的原因加以区分的方法与技术，其主要用途是研究外界因素或试验条件的改

6、变对试验结果影响是否显研究外界因素或试验条件的改变对试验结果影响是否显著。著。类型：单因素方差分析类型：单因素方差分析(One-way ANOVA)、双因素方差、双因素方差分析分析(Two-way ANOVA)或多元方差分析或多元方差分析(MANOVA)。方差分析的基本模型是线性模型，并假设随机变量是独方差分析的基本模型是线性模型，并假设随机变量是独立、正态和等方差的。立、正态和等方差的。方差分析是根据平方和的加和原理，利用方差分析是根据平方和的加和原理，利用 F 检验，进而检验，进而判断试验因素对试验结果的影响是否显著。判断试验因素对试验结果的影响是否显著。3.3.1 单因素方差分析单因素方

7、差分析# # 建立数据集建立数据集dfdf # #yield-scan()yield-scan()24 30 28 2624 30 28 2627 24 21 2627 24 21 2631 28 25 3031 28 25 3032 33 33 2832 33 33 2821 22 16 2121 22 16 21Treat - rep(paste(A, 1:5, Treat - rep(paste(A, 1:5, sepsep = ), rep(4, 5) = ), rep(4, 5)dfdf - - data.framedata.frame( Treat, yield)( Treat,

8、yield)# # 方差分析方差分析 #fit - fit TukeyHSDTukeyHSD(fit)(fit) TukeyTukey multiple comparisons of means multiple comparisons of means 95% family-wise confidence level 95% family-wise confidence levelFit: Fit: aovaov(formula = yield (formula = yield Treat, data = Treat, data = dfdf) )$Treat$Treat diff diff

9、 lwrlwr uprupr p p adjadjA2-A1 -2.5 -8.2 3.2 0.66A2-A1 -2.5 -8.2 3.2 0.66A3-A1 1.5 -4.2 7.2 0.92A3-A1 1.5 -4.2 7.2 0.92A4-A1 4.5 -1.2 10.2 0.15A4-A1 4.5 -1.2 10.2 0.15A5-A1 -7.0 -12.7 -1.3 0.01A5-A1 -7.0 -12.7 -1.3 0.01A3-A2 4.0 -1.7 9.7 0.24A3-A2 4.0 -1.7 9.7 0.24A4-A2 7.0 1.3 12.7 0.01A4-A2 7.0 1.

10、3 12.7 0.01A5-A2 -4.5 -10.2 1.2 0.15A5-A2 -4.5 -10.2 1.2 0.15A4-A3 3.0 -2.7 8.7 0.50A4-A3 3.0 -2.7 8.7 0.50A5-A3 -8.5 -14.2 -2.8 0.00A5-A3 -8.5 -14.2 -2.8 0.00A5-A4 -11.5 -17.2 -5.8 0.00A5-A4 -11.5 -17.2 -5.8 0.00-15-10-50510A5-A4A5-A3A4-A3A5-A2A4-A2A3-A2A5-A1A4-A1A3-A1A2-A195% family-wise confidenc

11、e levelDifferences in mean levels of Treat print(duncan.test(fit,Treat,alpha=0.05) print(duncan.test(fit,Treat,alpha=0.05)$statistics$statistics Mean CV MSerror Mean CV MSerror 26 9.9 6.7 26 9.9 6.7$means$means yield std.err r Min. Max. yield std.err r Min. Max.A1 27 1.3 4 24 30A1 27 1.3 4 24 30A2 2

12、4 1.3 4 21 27A2 24 1.3 4 21 27A3 28 1.3 4 25 31A3 28 1.3 4 25 31A4 32 1.3 4 28 33A4 32 1.3 4 28 33A5 20 1.3 4 16 22A5 20 1.3 4 16 22$groups$groups trt means M trt means M1 A4 32 a1 A4 32 a2 A3 28 ab2 A3 28 ab3 A1 27 b3 A1 27 b4 A2 24 b4 A2 24 b5 A5 20 c5 A5 20 c数据正态性、等方差的检验数据正态性、等方差的检验#数据正态性检验数据正态性检

13、验library( car )fit.2 shapiro.test(resid(lm(yield Treat, data = df ) Shapiro-Wilk normality testdata: resid(lm(yield Treat, data = df)W = 0.87, p-value = 0.0126 bartlett.test(yield Treat, data = df)Bartlett test of homogeneity of variancesdata: yield by TreatBartletts K-squared = 0.051, df = 4, p-val

14、ue = 0.99973.3.2 双因素方差分析双因素方差分析 df-read.csv(file=d4.3.2.csv, header=T) #读入数据读入数据library(reshape)df.2-melt(df, id=c(A) #进行数据重构进行数据重构colnames(df.2)2:3-c(B, yield) #变量重命名变量重命名# 方差分析方差分析 #fit (duncan.test(fit, A, alpha = 0.05) (duncan.test(fit, A, alpha = 0.05)$means$means yield std.err r Min. Max. yiel

15、d std.err r Min. Max.A1 74 2.3 3 71 77A1 74 2.3 3 71 77A2 91 2.3 3 90 92A2 91 2.3 3 90 92A3 70 2.3 3 59 80A3 70 2.3 3 59 80A4 79 2.3 3 75 82A4 79 2.3 3 75 82A5 64 2.3 3 60 67A5 64 2.3 3 60 67A6 84 2.3 3 82 86A6 84 2.3 3 82 86$groups$groups trt means M trt means M1 A2 91 a1 A2 91 a2 A6 84 ab2 A6 84 a

16、b3 A4 79 bc3 A4 79 bc4 A1 74 cd4 A1 74 cd5 A3 70 de5 A3 70 de6 A5 64 e6 A5 64 e (duncan.test(fit, B, alpha = 0.05) (duncan.test(fit, B, alpha = 0.05)$means$means yield std.err r Min. Max. yield std.err r Min. Max.B1 74 1.6 6 59 90B1 74 1.6 6 59 90B2 76 1.6 6 60 90B2 76 1.6 6 60 90B3 80 1.6 6 67 92B3

17、 80 1.6 6 67 92$groups$groups trt means M trt means M1 B3 80 a1 B3 80 a2 B2 76 ab2 B2 76 ab3 B1 74 b3 B1 74 bdf-read.csv(file=d4.3.3.csv, header=T)df.2-melt(df,id=c(A)colnames(df.2)2:3-c(B,y)fit duncan.testduncan.test( (fit,Afit,A, alpha=0.05)$groups, alpha=0.05)$groups trttrt means M means M1 A3 33

18、 a1 A3 33 a2 A2 30 b2 A2 30 b3 A1 28 c3 A1 28 c duncan.testduncan.test( fit, B, alpha=0.05)$groups( fit, B, alpha=0.05)$groups trttrt means M means M1 b4 32 a1 b4 32 a2 b3 32 a2 b3 32 a3 b1 30 b3 b1 30 b4 b2 29 c4 b2 29 c5 5b5 28 db5 28 d with(df.2,duncan.test(y,A:B,DFerror=45,MSerror=1.22)$groups)

19、with(df.2,duncan.test(y,A:B,DFerror=45,MSerror=1.22)$groups) trt means M trt means M1 A3:b3 35 a1 A3:b3 35 a2 A3:b4 34 ab2 A3:b4 34 ab3 A3:b2 34 ab3 A3:b2 34 ab4 A3:b1 33 b4 A3:b1 33 b5 A2:b4 33 b5 A2:b4 33 b6 A2:b3 31 c6 A2:b3 31 c7 A1:b3 30 cd7 A1:b3 30 cd8 A1:b4 30 cd8 A1:b4 30 cd9 A3:b5 30 cd9 A

20、3:b5 30 cd10 A2:b1 29 de10 A2:b1 29 de11 A2:b5 28 ef11 A2:b5 28 ef12 A1:b1 27 fg12 A1:b1 27 fg13 A2:b2 26 fgh13 A2:b2 26 fgh14 A1:b2 26 gh14 A1:b2 26 gh15 A1:b5 25 h15 A1:b5 25 h交互作用图形交互作用图形library(HH)interaction2wt( y A * B, data = df.2)y: main effects and 2-way interactionsA1 A2 A3A Ay A | Ab1b2b3

21、b4b5B B24262830323436y yy B | Ay A | B24262830323436y yy B | BBb1b2b3b4b5AA1A2A33.4协方差分析协方差分析(analysis of covariance) 协方差分析是关于如何调节协变量对因变量的影响协方差分析是关于如何调节协变量对因变量的影响效应，从而更加有效地分析实验处理效应的一种统计技效应，从而更加有效地分析实验处理效应的一种统计技术，也是对实验进行统计控制的一种综合方差分析和回术，也是对实验进行统计控制的一种综合方差分析和回归分析的方法。归分析的方法。 当研究者知道有些协变量会影响因变量，却不能够当研究者知

22、道有些协变量会影响因变量，却不能够控制和不感兴趣时，则可以在实验处理前予以观测，然控制和不感兴趣时，则可以在实验处理前予以观测，然后在统计时运用协方差分析来处理。将协变量对因变量后在统计时运用协方差分析来处理。将协变量对因变量的影响从自变量中分离出去，可以进一步提高实验精确的影响从自变量中分离出去，可以进一步提高实验精确度和统计检验灵敏度。例如林木生长量与肥料的关系，度和统计检验灵敏度。例如林木生长量与肥料的关系，施肥条件可以人工控制，但林木初始苗高施肥条件可以人工控制，但林木初始苗高(协变量协变量)是难是难以控制的，通过协方差分析，消除初始苗高的影响，使以控制的，通过协方差分析，消除初始苗高

23、的影响，使得生长量在一致的基础上进行方差分析。得生长量在一致的基础上进行方差分析。library(lsmeans)df-read.csv(file=d4.4.1.csv)fit-lm(y x+A, data=df)#fit2 summary(fit)Coefficients:Coefficients: Estimate Std. Error t value Pr(|t|) Estimate Std. Error t value Pr(|t|) (Intercept) 0.8516 0.1853 4.59 0.00177 (Intercept) 0.8516 0.1853 4.59 0.0017

24、7 * * * x 0.2226 0.0344 6.47 0.00019 x 0.2226 0.0344 6.47 0.00019 * * * *AA2 0.0110 0.1201 0.09 0.92950 AA2 0.0110 0.1201 0.09 0.92950 AA3 0.6468 0.1582 4.09 0.00349 AA3 0.6468 0.1582 4.09 0.00349 * * * -Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Residual standard error: 0.19 on 8 degrees of free

25、domMultiple R-squared: 0.875,Adjusted R-squared: 0.828 F-statistic: 18.6 on 3 and 8 DF, p-value: 0.000578 anova(fit) anova(fit)Analysis of Variance TableAnalysis of Variance TableResponse: yResponse: y Df Sum Sq Mean Sq F value Pr(F) Df Sum Sq Mean Sq F value Pr(F) x 1 1.308 1.308 36.75 0.0003 x 1 1

26、.308 1.308 36.75 0.0003 * * * *A 2 0.677 0.338 9.51 0.0077 A 2 0.677 0.338 9.51 0.0077 * * * Residuals 8 0.285 0.036 Residuals 8 0.285 0.036 -Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1 lsmeanslsmeans(fit, pairwise (fit, pairwise A) A)$A $A lsmeanslsmeans A A lsmeanlsmean SE SE dfdf lower.CL uppe

27、r.CL lower.CL upper.CL A1 1.9 0.084 8 1.8 2.1 A1 1.9 0.084 8 1.8 2.1 A2 2.0 0.085 8 1.8 2.2 A2 2.0 0.085 8 1.8 2.2 A3 2.6 0.134 8 2.3 2.9 A3 2.6 0.134 8 2.3 2.9$A pairwise differences$A pairwise differences estimate SE estimate SE dfdf t.ratiot.ratio p.valuep.valueA1 - A2 -0.011 0.12 8 -0.091 0.9954

28、A1 - A2 -0.011 0.12 8 -0.091 0.9954A1 - A3 -0.647 0.16 8 -4.090 0.0087A1 - A3 -0.647 0.16 8 -4.090 0.0087A2 - A3 -0.636 0.16 8 -3.983 0.0100A2 - A3 -0.636 0.16 8 -3.983 0.0100 p values are adjusted using the p values are adjusted using the tukeytukey method for 3 meansmethod for 3 means双因素协方差分析双因素协方

29、差分析【例子例子 4.4.2】为研究某杨树一年生生长与为研究某杨树一年生生长与 N 肥、肥、K 肥及初始苗高肥及初始苗高的关系，采用正交试验设计，共设置了的关系，采用正交试验设计，共设置了 18 个样地的栽培试验，试个样地的栽培试验，试验因子与水平及测量结果如表验因子与水平及测量结果如表 4-13所示。试分析所示。试分析 N 肥、肥、K 肥及初肥及初始苗高对生长量的影响。始苗高对生长量的影响。dfdf-read.csv(file=d4.4.2.csv, header=T)-read.csv(file=d4.4.2.csv, header=T)dfdf,3-,3-as.factoras.fact

30、or( (dfdf,3),3)fit - lm(mass height + N + K, data = fit - lm(mass height + N + K, data = dfdf ) )summary(fit)summary(fit)anovaanova(fit)(fit)lsmeanslsmeans(fit, pairwise N)(fit, pairwise N)lsmeanslsmeans(fit, pairwise K)(fit, pairwise K)$N pairwise differences$N pairwise differences estimate SE esti

31、mate SE dfdf t.ratiot.ratio p.valuep.valuehigh - low 0.097 0.035 13 2.8 0.016 high - low 0.097 0.035 13 2.8 0.016 $K pairwise differences$K pairwise differences estimate SE estimate SE dfdf t.ratiot.ratio p.valuep.value0 - 12.5 -0.071 0.043 13 -1.7 0.258500 - 12.5 -0.071 0.043 13 -1.7 0.258500 - 25

32、-0.308 0.043 13 -7.2 0.000020 - 25 -0.308 0.043 13 -7.2 0.0000212.5 - 25 -0.237 0.043 13 -5.5 0.0002812.5 - 25 -0.237 0.043 13 -5.5 0.000283.5 t检验检验 在实际工作中，经常需要判断两个样本平均数在实际工作中，经常需要判断两个样本平均数是否差异，以了解两样本所属的两个总体平均数是是否差异，以了解两样本所属的两个总体平均数是否相同。检验方法可以使用否相同。检验方法可以使用 t 检验。检验。 对于两个样本平均数差异显著性检验，可分为对于两个样本平均数差异显著

33、性检验，可分为非配对设计和配对设计。非配对设计和配对设计。height shapiro.testshapiro.test(height)(height)Shapiro-Shapiro-WilkWilk normality test normality testdata: heightdata: heightW = 0.94, p-value = 0.06344W = 0.94, p-value = 0.06344-2-10127.58.08.59.0Normal Q-Q PlotTheoretical QuantilesSample Quantiles t.test(height, mu =

34、8, alternative = two.sided)One Sample t-testdata: heightt = 4.5, df = 31, p-value = 8.492e-05alternative hypothesis: true mean is not equal to 895 percent confidence interval: 8.2 8.5sample estimates:mean of x 8.33.5.2 双样本双样本 t 检验检验weight-scan()16.68 20.67 18.42 1817.44 15.9518.68 23.22 21.42 1918.9

35、2 NAVariety-rep(c(LY1,DXY), rep(6,2)df-data.frame(Variety, weight)a-subset(df$weight, Variety = LY1)b var.test(a,b)F test to compare two variancesdata: a and bF = 0.67, num df = 5, denom df = 4, p-value = 0.6653alternative hypothesis: true ratio of variances is not equal to 195 percent confidence in

36、terval: 0.072 4.971sample estimates:ratio of variances 0.67 t.testt.test( (a,ba,b, paired = F), paired = F)Welch Two Sample t-testWelch Two Sample t-testdata: a and bdata: a and bt = -2.1, t = -2.1, dfdf = 7.8, p-value = 0.06591 = 7.8, p-value = 0.06591alternative hypothesis: true difference in mean

37、s is alternative hypothesis: true difference in means is not equal to 0not equal to 095 percent confidence interval:95 percent confidence interval: -5.0 0.2 -5.0 0.2sample estimates:sample estimates:mean of x mean of y mean of x mean of y 18 20 18 20 3.5.3 成对双样本成对双样本 t 检验检验height - scan() 7.127.264.

38、787.694.254.966.284.82 6.526.074.287.304.175.665.734.52 Variety - rep( c(A, B), rep(8, 2)df - data.frame(Variety, height) a - subset( df$height, Variety = A ) # 选取品种选取品种A的树高的树高b t.testt.test( (a,ba,b, paired=T), paired=T) Paired t-test Paired t-testdata: a and bdata: a and bt = 1.9, t = 1.9, dfdf =

39、7, p-value = 0.09624 = 7, p-value = 0.09624alternative hypothesis: true difference in means is alternative hypothesis: true difference in means is not equal to 0not equal to 095 percent confidence interval:95 percent confidence interval: -0.084 0.812 -0.084 0.812sample estimates:sample estimates:mea

40、n of the differences mean of the differences 0.36 0.36 32 种类放牧不放牧Asclepias syriaca0.0340.247Aster laevis 0.2440.096Aster lateriflorus 0.0410.146Aster novae-angliae 0.3100.365Aster simplex0.0620.088Dactylis glomerata0.0010.055Fragaria virginiana0.4410.385Hieracium pratense0.5920.626Phleum pratense0.3

41、870.911Picris hieracoides1.3691.510Plantago lanceolata 0.2600.208Poa compressa0.6100.773Poa pratensis0.0540.116Solidago altissima0.8431.967Solidago graminifolia0.2010.097Solidago juncea0.2780.148Solidago rugosa0.1560.197Taraxacum officinale0.1000.151N=1818个草地种在放牧和不放牧样方中的生物量（个草地种在放牧和不放牧样方中的生物量（kg/m2）

42、放牧对所研究草地物种生物量的影响是否显著？放牧对所研究草地物种生物量的影响是否显著？t=read.csv(t_test.csv,header=T)head(t)t.test(t$Graze,t$Control,paired=T) t.testt.test( (t$Graze,t$Control,pairedt$Graze,t$Control,paired=T)=T) Paired t-test Paired t-testdata: data: t$Grazet$Graze and and t$Controlt$Controlt = -1.7, t = -1.7, dfdf = 17, p-va

43、lue = 0.1097 = 17, p-value = 0.1097alternative hypothesis: true difference in means is alternative hypothesis: true difference in means is not equal to 0not equal to 095 percent confidence interval:95 percent confidence interval: -0.263 0.029 -0.263 0.029sample estimates:sample estimates:mean of the

44、 differences mean of the differences -0.12 -0.12 3.5 卡方检验卡方检验( 2 test)卡方检验是参照卡方分配来计算概率和临界值的统计检验，是用途很广卡方检验是参照卡方分配来计算概率和临界值的统计检验，是用途很广的一种假设检验方法。的一种假设检验方法。 分析原理：分析原理： (1) 建立零假说（建立零假说（Null Hypothesis），即认为观测值与理论值的差异是），即认为观测值与理论值的差异是由于随机误差所致；由于随机误差所致； (2) 确定数据间的实际差异，即求出确定数据间的实际差异，即求出 2 值；值； (3) 如卡方值大于某特定概

45、率标准（即显著性差异）下的理论值，则拒如卡方值大于某特定概率标准（即显著性差异）下的理论值，则拒绝零假说，即实测值与理论值的差异在该显著性水平下是显著的。绝零假说，即实测值与理论值的差异在该显著性水平下是显著的。卡方检验35freq = c(22,21,22,27,22,36)probs = c(1,1,1,1,1,1)/6 chisq.test(freq,p=probs) chisq.test(freq,p=probs) Chi-squared test for given probabilitiesdata: freqX-squared = 6.7, df = 5, p-value = 0

46、.2423卡方检验36x = c(100,110,80,55,14)probs = c(29, 21, 17, 17, 16)/100 chisq.test(x,p=probs) chisq.test(x,p=probs) Chi-squared test for given probabilitiesdata: xX-squared = 55, df = 4, p-value = 2.685e-11卡方检验（列联表）37yesbelt = c(12813,647,359,42)nobelt = c(65963,4000,2642,303)chisq.test(data.frame(yesbe

47、lt,nobelt) chisq.test(data.frame(yesbelt,nobelt) Pearsons Chi-squared testdata: data.frame(yesbelt, nobelt)X-squared = 59, df = 3, p-value = 8.61e-13练习四练习四以数据以数据stu.data.csv为例，为例，试对体重做频数分析。试对体重做频数分析。请分析身高是否符合正态分布？请分析身高是否符合正态分布？试分析性别对体重有无影响。试分析性别对体重有无影响。问题问题4：请检验总体平均体重与：请检验总体平均体重与60kg有无显著差有无显著差异？男生和女

48、生的平均体重有无显著差异？异？男生和女生的平均体重有无显著差异？问题问题5：男女生比例是否符合：男女生比例是否符合 1.2 : 1.0？练习四练习四 答案答案df-read.csv(file=stu.data.csv,header=T)#问题问题1A - table(cut(df$weight, breaks = 40 + 15 * (0:7) round(prop.table(A) * 100,2) # 计算频数比例计算频数比例hist(df$weight, breaks = 7, xlim = c(40,140), xlab = weight, main = Frequency chart

49、 of weight)#问题问题2shapiro.test(df$height)#问题问题3fit-aov(weight Sex,data=df)summary(fit)library(agricolae)duncan.test(fit, Sex, alpha=0.05)$groups#问题问题4t.test(df$weight, mu = 60, alternative = two.sided)wt.m-subset(df$weight,df$Sex2=1)wt.f-subset(df$weight,df$Sex2=2)var.test(wt.m,wt.f) #等方差检验等方差检验t.tes

50、t(wt.m,wt.f, paired=F) #问题问题5summary(df$Sex)ct-c(87, 33)pt - c(1.2/2.2, 1.0/2.2) chisq.test(ct, p = pt)3.6 线性回归线性回归 比如产量与施肥量有关，病虫害发生时期与气温有比如产量与施肥量有关，病虫害发生时期与气温有关，小麦单位面积产量与单位面积穗数、千粒重有关，关，小麦单位面积产量与单位面积穗数、千粒重有关，等等。因此，还需要研究两个或多个变量之间的关系。等等。因此，还需要研究两个或多个变量之间的关系。 一个变量的变化受另一个或几个变量的影响，称为一个变量的变化受另一个或几个变量的影响，称

51、为因果关系。因果关系。利用回归分析利用回归分析(regression analysis)来研究呈因果关系来研究呈因果关系的变量之间的关系。表示原因的变量为自变量，表示结的变量之间的关系。表示原因的变量为自变量，表示结果的变量为因变量。回归分析有一元和多元回归分析。果的变量为因变量。回归分析有一元和多元回归分析。3.6.1 简单线性回归简单线性回归dfdf - read.csv( file = d4.7.1.csv, header = T) # - read.csv( file = d4.7.1.csv, header = T) # 读入数据读入数据 fit - lm( weight N, da

52、ta = fit summary(fit)Coefficients: Estimate Std. Error t value Pr(|t|) Estimate Std. Error t value Pr(|t|) (Intercept) -87.5167 5.9369 -14.7 1.7e-09 (Intercept) -87.5167 5.9369 -14.7 1.7e-09 * * * *N 3.4500 0.0911 37.9 1.1e-14 N 3.4500 0.0911 37.9 1.1e-14 * * * *-Signif. codes: 0 * 0.001 * 0.01 * 0.

53、05 . 0.1 1Residual standard error: 1.5 on 13 degrees of freedomMultiple R-squared: 0.991, Adjusted R-squared: 0.99 F-statistic: 1.43e+03 on 1 and 13 DF, p-value: 1.09e-14 df$weight 1 115 117 120 123 126 129 132 135 139 142 146 150 154 159 164 fitted(fit) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 2 3 4 5

54、 6 7 8 9 10 11 12 13 14 15 113 116 119 123 126 130 133 137 140 144 147 151 154 157 161 113 116 119 123 126 130 133 137 140 144 147 151 154 157 161 residuals(fit) 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 2.417 0.967 0.517 0.067 -0.383 -0.833 -1.283 -1.733 -1.183 2.417 0.967 0.517 0.067 -0.383 -0.833 -1.28

55、3 -1.733 -1.183 10 11 12 13 14 15 10 11 12 13 14 15 -1.633 -1.083 -0.533 0.017 1.567 3.117 -1.633 -1.083 -0.533 0.017 1.567 3.117 5860626466687072120140160df$Ndf$weight计算校正R243这里这里n是对象的数量（样点数量），是对象的数量（样点数量），m是解释变量的是解释变量的数量（或更准确说，是模型的自由度）。只有模型数量（或更准确说，是模型的自由度）。只有模型自由度（自由度（m）不要比观测值的数量（）不要比观测值的数量（n）大（保

56、证）大（保证n-m-10）,公式是有效的。公式是有效的。3.6.2 多项式回归多项式回归fit summary(fit2)Coefficients: Estimate Std. Error t value Pr(|t|) Estimate Std. Error t value Pr(|t|) (Intercept) 261.87818 25.19677 10.39 2.4e-07 (Intercept) 261.87818 25.19677 10.39 2.4e-07 * * * *N -7.34832 0.77769 -9.45 6.6e-07 N -7.34832 0.77769 -9.4

57、5 6.6e-07 * * * *I(N2) 0.08306 0.00598 13.89 9.3e-09 I(N2) 0.08306 0.00598 13.89 9.3e-09 * * * *-Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Residual standard error: 0.38 on 12 degrees of freedomMultiple R-squared: 0.999,Adjusted R-squared: 0.999 F-statistic: 1.14e+04 on 2 and 12 DF, p-value: 2e-1

58、65860626466687072120140160df$Ndf$weight3.6.3 多元线性回归多元线性回归dfdf- read.csv(file=d4.7.3.csv, header=T)- read.csv(file=d4.7.3.csv, header=T)lmfitlmfit-lm(yx1+x2+x3+x4,data=-lm(yx1+x2+x3+x4,data=dfdf) )step(lm(yx1+x2+x3+x4,data=step(lm(yx1+x2+x3+x4,data=dfdf)lmfit2lmfit2-lm(yx1+x2+x3,data= summary(lmfit)C

59、oefficients: Estimate Std. Error t value Pr(|t|) Estimate Std. Error t value Pr(|t|) (Intercept) -625.358 114.378 -5.47 6.5e-05 (Intercept) -625.358 114.378 -5.47 6.5e-05 * * * *x1 15.196 2.127 7.15 3.4e-06 x1 15.196 2.127 7.15 3.4e-06 * * * *x2 7.378 1.889 3.91 0.0014 x2 7.378 1.889 3.91 0.0014 * *

60、 * x3 9.503 1.342 7.08 3.7e-06 x3 9.503 1.342 7.08 3.7e-06 * * * *x4 -0.847 1.493 -0.57 0.5790 x4 -0.847 1.493 -0.57 0.5790 -Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Residual standard error: 36 on 15 degrees of freedomMultiple R-squared: 0.894,Adjusted R-squared: 0.866 F-statistic: 31.8 on 4 an