R语言基础培训第二讲常用统计分析

1、基于R的基本统计分析,内容提要,描述统计 频数表分析 方差分析 t检验,卡方检验 线性回归 相关分析,描述分析(Descriptive statistics),描述统计就是把数据集所包含的信息加以简要地概况，如计算数据的数字特征、制作频数表和频数图等等，用所获得的统计量和图表来描述数据集所反映的特征和规律，使得研究的问题更加简单、直观。 描述性统计主要包括反映数据集中趋势的特征值(比如平均数、中位数、众数、分位数)、数据离散程度的特征值(比如方差、标准差、值域、变异系数)和数据分布形态的特征值(比如偏度、峰度)。,标准差（std.dev）和标准误（SE.mean）,标准差（std.dev）,真

2、实均值,SE,样本均值,当样本含量 n 足够大时，标准差趋向稳定；而标准误随n的增大而减小，甚至趋于0 。,标准差（std.dev）和标准误（SE.mean）,比如，某学校共有500名学生，现在要通过抽取样本量为30的一个样本，来推断学生的身高。这时可以依据抽取的样本信息，计算出样本的均值与标准差。如果我们抽取的不是一个样本，而是10个样本，每个样本30人，那么每个样本都可以计算出均值，这样就会有10个均值。也就是形成了一个10个数字的数列，然后计算这10个数字的标准差，此时的标准差就是标准误。但是，在实际抽样中我们不可能抽取10个样本。所以，标准误就由样本标准差除以样本量来表示。当然，这样的

3、结论也不是随心所欲，而是经过了统计学家的严密证明的。 SE.mean= std.dev / n1/2,在实际的应用中，标准差主要有两点作用，一是统计量样本离散程度的表征；二是用来对样本进行标准化处理，即样本观察值减去样本均值，然后除以标准差，这样就变成了标准正态分布。标准误的作用主要是用来做区间估计，常用的估计区间是均值加减n倍的标准误（例如95%的置信区间是：均值+1.96*SE）,标准差（std.dev）和标准误（SE.mean）,95% CI: 假设上面这个随机抽样估计学生身高的例子，抽样100次，每次抽10个学生测量身高，均值估计值及标准误为 152cm12cm。但有时需要表示为估计量

4、的95%的置信区间152cm-1.96*12cm，152cm+1.96*12cm。可以解释为，如果从再从总体中抽样100次（每次抽样10个），产生100个平均值，这100个平均值将有95次落在152cm-1.96*12cm，152cm+1.96*12cm这个范围内，5次落在这个范围外，如果抽样次数越多，这个推断越准确。这个来源于中心极限定理的应用：任何分布(总体)抽样n次，每次抽样的和符合正态分布。通俗一点说, 不管是学校的学生身高是怎么分布, 每次随机抽取10个求和, 抽取n次，这n个身高总和是符合正态分布的。平均身高为身高总和除于10，所以平均身高也是正态分布的。正态分布双尾95%的分界点

5、所对应的值刚好是1.96。,峰度（Kurtosis）,峰度（Kurtosis）是描述某变量所有取值分布形态陡缓程度的统计量。 它是和正态分布相比较的。 Kurtosis=0 与正态分布的陡缓程度相同。 Kurtosis0 比正态分布的高峰更加陡峭尖顶 Kurtosis0 比正态分布的高峰来得平台平顶 计算公式：,偏度（Skewness）,偏度（Skewness）是描述某变量分布偏离正态对称性的统计量。 Skewness=0 分布形态与正态分布偏度相同 Skewness0 正偏差数值较大，为正偏或右偏。长尾巴拖在右边 Skewness0 负偏差数值较大，为负偏或左偏。长尾巴拖在左边 计算公式：

6、Skewness 越大，分布形态偏移程度越大,频数表(Frequency table)分析,频数表分析是对数据集按数据范围分成若干区间，即分成若干组，求出每组组中值，各组数据用组中值代替，计算各组数据的频数，并作出频数表。,频数表分析例子,summary(oats$yield) # 计算频数 A - table(cut(oats$yield, breaks = 40 + 20 * (0:7) round(prop.table(A) * 100,2) # 计算频数比例 # 画频数表 hist(oats$yield, # breaks = 7, xlim = c(40,180), xlab = y

7、ield, main = Frequency chart of yield),方差分析ANOVA,方差分析是一种在若干组能相互比较的试验数据中，把产生变异的原因加以区分的方法与技术，其主要用途是研究外界因素或试验条件的改变对试验结果影响是否显著。 类型：单因素方差分析(One-way ANOVA)、双因素方差分析(Two-way ANOVA) 。 方差分析的基本模型是线性模型，并假设随机变量是独立、正态和等方差的。 方差分析是根据平方和的加和原理，利用 F 检验，进而判断试验因素对试验结果的影响是否显著。,单因素方差分析,# Tukey HSD 方法 # install.packages(mu

8、ltcomp) library(multcomp) tuk - glht(fit, linfct = mcp(Treat = Tukey) summary(tuk) # standard display tuk.cld - cld(tuk) # letter-based display opar - par(mai=c(1,1,1.5,1) plot(tuk.cld) par(opar),双因素（无重复）方差分析,多重比较,library(agricolae) # 对A因素在a = 0.05水平上进行多重比较 (duncan.test(fit, A, alpha = 0.05) # 对B因素进

9、行多重比较 (duncan.test(fit, B, alpha = 0.05),重复试验的双因素方差分析,协方差分析(analysis of covariance),协方差分析是关于如何调节协变量对因变量的影响效应，从而更加有效地分析实验处理效应的一种统计技术，也是对实验进行统计控制的一种综合方差分析和回归分析的方法。 当研究者知道有些协变量会影响因变量，却不能够控制和不感兴趣时，则可以在实验处理前予以观测，然后在统计时运用协方差分析来处理。将协变量对因变量的影响从自变量中分离出去，可以进一步提高实验精确度和统计检验灵敏度。例如林木生长量与肥料的关系，施肥条件可以人工控制，但林木初始苗高(协

10、变量)是难以控制的，通过协方差分析，消除初始苗高的影响，使得生长量在一致的基础上进行方差分析。,单因素协方差分析,双因素协方差分析,【例子 4.4.2】为研究某杨树一年生生长与 N 肥、K 肥及初始苗高的关系，采用正交试验设计，共设置了 18 个样地的栽培试验，试验因子与水平及测量结果如表 4-13所示。试分析 N 肥、K 肥及初始苗高对生长量的影响。,t检验,在实际工作中，经常需要判断两个样本平均数是否差异，以了解两样本所属的两个总体平均数是否相同。检验方法可以使用 t 检验。 对于两个样本平均数差异显著性检验，可分为非配对设计和配对设计。,3.5.1 单样本检验,【例子 4.5.1】 杨树

11、某无性系试验林造林 5 年后，调查树高生长量，随机抽取 32 棵树，调查结果如下表 4-16 所示。有一无性系 B5 的 5 年树龄树高 = 8 m。试分析该试验林的树高与 B5 有无显著差异？,25,18个草地种在放牧和不放牧样方中的生物量（kg/m2）,放牧对所研究草地物种生物量的影响是否显著？,成对双样本 t 检验,卡方检验( 2 test),卡方检验是参照卡方分配来计算概率和临界值的统计检验，是用途很广的一种假设检验方法。 分析原理： (1) 建立零假说（Null Hypothesis），即认为观测值与理论值的差异是由于随机误差所致； (2) 确定数据间的实际差异，即求出 2 值； (

12、3) 如卡方值大于某特定概率标准（即显著性差异）下的理论值，则拒绝零假说，即实测值与理论值的差异在该显著性水平下是显著的。,27,卡方检验,freq = 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 probabilities data: freq X-squared = 6.7, df = 5, p-value = 0.2423,28,卡方检验,x = c(100,110,80,55,14) p

13、robs = c(29, 21, 17, 17, 16)/100 chisq.test(x,p=probs) chisq.test(x,p=probs) Chi-squared test for given probabilities data: x X-squared = 55, df = 4, p-value = 2.685e-11,29,卡方检验（列联表）,yesbelt = c(12813,647,359,42) nobelt = c(65963,4000,2642,303) chisq.test(data.frame(yesbelt,nobelt) chisq.test(data.f

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

15、) #问题1 A - 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 of weight) #问题2 shapiro.test(df$height) #问题3 fit-aov(weight Sex,data=df) summary(fit) library(agricolae) duncan.test(fit,

16、 Sex, alpha=0.05)$groups,#问题4 t.test(df$weight, mu = 60, alternative = two.sided) wt.m-subset(df$weight,df$Sex=男) wt.f-subset(df$weight,df$Sex=女) var.test(wt.m,wt.f) #等方差检验 t.test(wt.m,wt.f, paired=F) #问题5 summary(df$Sex) ct-c(87, 33) pt - c(1.2/2.2, 1.0/2.2) chisq.test(ct, p = pt),线性回归,比如产量与施肥量有关，病

17、虫害发生时期与气温有关，小麦单位面积产量与单位面积穗数、千粒重有关，等等。因此，还需要研究两个或多个变量之间的关系。 一个变量的变化受另一个或几个变量的影响，称为因果关系。 利用回归分析(regression analysis)来研究呈因果关系的变量之间的关系。表示原因的变量为自变量，表示结果的变量为因变量。回归分析有一元和多元回归分析。,3.6.1 简单线性回归,35,计算校正R2,这里n是对象的数量（样点数量），m是解释变量的数量（或更准确说，是模型的自由度）。只有模型自由度（m）不要比观测值的数量（n）大（保证n-m-10）,公式是有效的。,多项式回归,fit - lm( weight

18、N + I(N2), data = df ) # 二次项回归模型 summary( fit ) # 回归分析结果 plot( df$N, df$weight ) # 绘制散点图 lines( df$N, fitted(fit) ) # 添加回归线 添加光滑曲线 require(graphics) plot(cars, main = lowess(cars) lines(cars) lines(lowess(cars), col = 2),多元线性逐步回归,最优线性回归方程为： y= -649.779 + 14.592x1 + 6.841x2 + 9.329x3 回归方程表明：对于川农 16 号

19、小麦而言，当 x2 和 x3 固定时，穗数 x1 每增加 1万/亩，产量 y 将平均增加 14.592Kg/亩；当 x1 和 x3 固定时，每穗粒数 x2 每增加 1粒，产量 y 将平均增加 6.841Kg/亩；当 x1 和 x2 固定时，千粒重 x3 每增加 1g，产量y 将平均增加 9.329 Kg/亩。,自变量的作用主次,第一种方法：计算通径系数(path coefficient, p)。 df.2 round(coef(lmfit),3) (Intercept) x1 x2 x3 0.000 0.777 0.410 0.609 第二种方法：根据 x1、x2、x3 回归系数的 t 值大小

20、，也可判断自变量 x 的主次顺序，凡是 t 值较大者就是较重要的因子。 结论： 自变量 x 的主次顺序为：穗数 x1 千粒重 x3 每穗粒数 x2。,交互作用多元线性回归,练习四,以数据stu.data.csv为例， 试对体重做频数分析。 请分析身高是否符合正态分布？ 试分析性别对体重有无影响。 问题4：请检验总体平均体重与60kg有无显著差异？男生和女生的平均体重有无显著差异？ 问题5：男女生比例是否符合 1.2 : 1.0？,练习四 答案,df-read.csv(file=stu.data.csv,header=T) #问题1 A - table(cut(df$weight, breaks

21、 = 40 + 15 * (0:7) round(prop.table(A) * 100,2) # 计算频数比例 hist(df$weight, breaks = 7, xlim = c(40,140), xlab = weight, main = Frequency chart of weight) #问题2 shapiro.test(df$height) #问题3 fit-aov(weight Sex,data=df) summary(fit) library(agricolae) duncan.test(fit, Sex, alpha=0.05)$groups,#问题4 t.test(d

22、f$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.test(wt.m,wt.f, paired=F) #问题5 summary(df$Sex) ct-c(87, 33) pt - c(1.2/2.2, 1.0/2.2) chisq.test(ct, p = pt),相关分析,相关关系是指两个变量的数值变化存在不完全确定的依存关系，它们之间的数值不能用方程表示出来，但可用

23、某种相关性度量来描述。 按照相关程度: 完全相关、不完全相关和不相关； 按照相关方向: 正相关和负相关； 按照相关形式: 线性相关和非线性相关； 按涉及变量的多少: 一元相关和多元相关； 按影响因素: 单相关和复相关。,Pearson、Spearman和Kendall相关,Pearson 相关是积差相关，衡量两个定量变量之间的线性相关程度，可描述两个正态分布变量间线性相关关系的密切程度。 Spearman 相关是等级相关，衡量分级定序变量之间的相关程度。 Kendall 相关，也是一种非参数的等级相关。 当两变量不符合正态分布时，变量间的关系应通过计算 Spearman 或 Kendall 相

24、关来考察。 cor( )函数可以计算上述三种相关系数，cov( )函数可以计算协方差。 这两个函数的使用模板为： cor( X, use = ”A”, method = ”B” ) X 是数据集；默认情况下，A 是”everything”(所有观测值)， B 是”pearson”。,46,Pearson积矩相关系数,相关系数r的显著性检验可以用t检验,47,Spearman,48,49,【例子 】 在某杉木林分内随机抽取 30 棵树，测定了树高 h (m)、胸径 dbh (cm)、材积 v (m3)、心材比例 cpro (%)、木材基本密度 wd (kg/m3)、木材吸水率 wp

25、ro (%)、管胞长度 tl (m)、管胞宽度 tw (m)和管胞长宽比 lrt，测试结果如表 4-29 所示。试分析各性状之间的相关系数。,偏相关(Partial correlation),偏相关是描述在控制一个或多个定量变量保持不变时，指定的两个变量之间的相关关系。 ggm 包中的 pcor( )函数可以计算偏相关。 函数 pcor( )的使用格式为： pcor(u, S) u是一个数值向量，前两个数值表示要计算相关系数的变量的下标，其余的数值为保持不变的变量的下标；S是变量的协方差矩阵。 函数 pcor .test( )的使用格式为： pcor .test(r , q, n) r 是偏相

26、关系数值，q 是控制变量的数量，n 是样本大小。,仍以例子的数据集为例，进行偏相关系数计算，并检验其显著性。 library(ggm) # source(/biocLite.R) pcor(c(1,3,2),cov(df) pcor.test(0.58,1,30),相关显著性的检验,cor.test( )函数可以对单个 pearson、spearman 和 kendall 相关系数进行检验。 cor.test( )函数的使用方法如下： cor.test(x, y, alternative = , method = ),虽然 cor.te

27、st( )函数可以检验相关的显著性，但每次只能检验一个相关值。如果需要检验多个相关时，可使用 psych 包中的corr.test( )函数。 corr.test( )函数的使用方法如下： corr.test(x, y, use = , method = ), corr.test(df, use = complete) Call:corr.test(x = df, use = complete) Correlation matrix h dbh v cpro wd wpro tl tw lrt h 1.00 0.86 0.90 -0.43 -0.29 0.32 -0.20 -0.04 -0.1

28、7 dbh 0.86 1.00 0.98 -0.41 -0.35 0.39 -0.18 -0.15 -0.07 v 0.90 0.98 1.00 -0.43 -0.38 0.43 -0.20 -0.11 -0.12 cpro -0.43 -0.41 -0.43 1.00 0.22 -0.25 0.18 -0.08 0.20 wd -0.29 -0.35 -0.38 0.22 1.00 -0.98 -0.09 -0.08 0.00 wpro 0.32 0.39 0.43 -0.25 -0.98 1.00 0.02 -0.02 0.01 tl -0.20 -0.18 -0.20 0.18 -0.09 0.02 1.00 0.16 0.79 tw -0.04 -0.15