《回归分析与相关分析精选课件.pptx》由会员分享，可在线阅读，更多相关《回归分析与相关分析精选课件.pptx（52页珍藏版）》请在三一办公上搜索。
1、课程安排回归模型与相关分析方差分析属性数据分析生存数据分析实验设计全基因组关联分析,课程安排,第1,2,3,5次Experimental Design and Data Analysis for Biologists Gerry P. Quinn and Michael J. Keough Cambridge University Press, Cambridge, 2019.2.Biostatistical Design and Analysis Using R: A Practical Guide Murray LoganWiley-Blackwell, Chichester, West
2、Sussex, 2019.第4次3. Survival Analysis: A Self-Learning Text, 2nd edition David G. Kleinbaum and Mitchel Klein, Springer, 2019第6次ABEL tutorialYurii Aulchenko, 2019,参考书,第1,2,3,5次参考书,回归分析和相关分析生物统计学研究所 张洪线性回归1.1回归直线1.2参数估计方法：最小二乘估计1.3 参数推断：F-检验、t-检验、区间估计1.4 回归模型的诊断：残差分析1.5 多重线性回归1.6 Box-Cox变换2. 相关分析2.1 P
3、earson相关系数 2.2 相关系数的统计推断：z-检验、区间估计2.3 偏相关系数2.4 秩相关与多重相关3 实例分析,回归分析和相关分析,与函数关系的区别：同一身高可以体重不同，体重随身高增加的关系不是严格成立，只是有这种趋势,总的趋势：随着身高增加，体重也跟着增加,为什么?因为身高只能解释体重的一部分原因,还有其他未被考虑的因素,如饮食、地域、人种等,与函数关系的区别：总的趋势：为什么?,线性回归目标：建立一个连续型因变量X(身高)与自变量Y(体重)之间的关系,因变量,自变量,随机误差: 未被身高解释的部分： 饮食、种族、地域,截距,斜率: (1) 0 正相关 (2) 0 负相关 (3
4、) = 0 不相关,X 每增加一个单位，Y 增加 个单位,均值为0:,最简单的关系-线性关系：,自变量和因变量的选择：由实际问题本身决定。,线性回归因变量自变量随机误差: 未被身高解释的部分：截距斜率,问题：如何估计参数？ X能解释Y的比例是多少？ 假设检验问题：线性模型是否适用？如果模型不适用，如何修正？,1.1 回归直线 量化X与Y的关系 对于新的个体，如果知道 与 ，可用X来预测Y。,问题：1.1 回归直线,1.2 参数估计方法：最小二乘估计 (Least Squares 估计),观测值与预测值（蓝线）的误差,观测值与预测值（黑线）的误差,总体上看，蓝线比黑线拟合效果要好（从拟合误差看）
5、,两条回归直线,1.2 参数估计方法：最小二乘估计 (Least Squar,观测值：残差：, fit = lm(weightheight, data=dat); print(fit);Call:lm(formula = weight height, data = dat)Coefficients:(Intercept) height -88.6774 0.8902,最小二乘估计,dat 是一个data frame，有两个变量：height 和 weight,R中线性模型拟合,(2) 最小一乘估计,估计方法,R代码,输出的结果, print(fit$coef); Coefficients:(I
6、ntercept) x -88.6774 0.8902,几个有用的函数：summary、names,观测值： fit = lm(weightheight,1.3 参数推断预测值：,总平方和SStot自由度=n-1,回归平方和SSreg自由度=1,残差平方和SSerr自由度=n-2,回归平方和占总平方和的比重越大,则X能解释Y的部分越大！回归方程越好！, summary(fit)$r.squared1 0.693347,方差分解,度量方法,1.3 参数推断总平方和SStot回归平方和SSreg残差平,F检验,适用条件：误差e1,en独立同分布，服从正态分布。, anova(fit)Analysi
7、s of Variance TableResponse: y Df Sum Sq Mean Sq F value Pr(F) x 1 523.71 523.71 29.393 0.0001168 *Residuals 13 231.63 17.82 -Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.,自由度,平方和,平均平方,F,P-值,如果误差独立同分布服从正态分布，则零假设成立时，F服从自由度为1和n-2的F分布,方差分析表,回归部分,残差部分,F统计量,F检验适用条件：误差e1,en独立同分布，服从正态分布。,t 检验: 回归系数的显著性检验零假设：
8、回归系数等于0, summary(fit)$coef Estimate Std. Error t value Pr(|t|)(Intercept) -88.6774259 28.3051787 -3.132905 0.0079283231x 0.8901553 0.1641884 5.421548 0.0001167570,参数估计,标准差,T 统计量,-值,零假设成立且误差独立同分布且服从正态分布时，T 服从自由度为n-2的t分布tn-2,适用条件：误差 e 独立同分布服从正态分布,样本量不太小时，误差分布偏离正态分布不是很大时仍适用，即对正态性假设不是很敏感,当样本量大时总是适用的，因为此
9、时 T 依分布收敛于标准正态分布与tn-2接近(不管正态假设成立与否)。,t 检验: 回归系数的显著性检验 summary(fit),区间估计, confint(fit) 2.5% 97.5% (Intercept) -149.8270467 -27.527805x 0.5354479 1.244863,置信系数 的双侧置信区间， 常取为0.05,结论： 的95%置信区间是 (0.5354479, 1.244863),落在(0.5354479, 1.244863)的概率是95%？错！,样本的两重性：抽样前是随机的，抽样后是非随机的。,区间估计 confint(fit)置信系数,1.4 回归模型
10、的诊断：残差分析,误差项独立同分布且服从正态分布吗?有没有需要剔除的 强影响点?,par(mfrow=c(2,2);plot(fit);,满足下述条件则 回归模型适用：残差没有明显的趋势残差的QQ图表明正态性成立Cook距离1 且leverage2p/n=2/15,1.4 回归模型的诊断：残差分析误差项独立同分布且服从正态,强影响点,强影响点：对回归方程有较大影响的点如位于散点图的两端且对于斜率的影响大，则应被剔出；位于中上、中下的点对截距影响大，但对斜率影响小，关系不大。,有强影响点后的回归直线,没有强影响点的回归直线,强影响点强影响点：对回归方程有较大影响的点有强影响点没有强影,残差散点图
11、残差Q-Q图位置(拟合值) -尺度图4.Cook距离与杠杆图 Cook距离1或 leverage2p/n为强影响点,残差散点图,发现并剔除残差大的强影响点有些强影响点残差不大，不应被提出有些残差大的点不是强影响点，剔除与否关系不大, dat1 = rbind(dat,c(180,200); fit1 = lm(weightheight,data=dat1); fit1.hat = lm.influence(fit1)$hat; influence1 = cbind(dat1,lev=fit1.hat); print(influence1fit1.hat2/n,);,杠杆(leverage)，由
12、“帽子矩阵”计算得到，其值越大对回归方程的影响力越大。,height weight lev1 180 68 0.10842272 160 51 0.249673011 165 57 0.135888012 185 76 0.203897514 164 58 0.154459816 190 200 0.3516872,第16个观测值杠杆很大,发现并剔除残差大的强影响点 dat1 = rbind(da,稳健回归法：最小一乘估计 (最小绝对偏差，LAD),黑线为LAD回归直线, require(quantreg); fit1.lad = rq(weightheight,data=dat1);警告信息
13、：In rq.fit.br(x, y, tau = tau, .) : Solution may be nonunique lines(dat1$height, fit1.lad$coef1+ fit1.lad$coef2*dat1$height, type=l,lwd=2,col=black);,稳健回归法：最小一乘估计 (最小绝对偏差，LAD)黑线为LA,没有一个方法可以完全打败其他的方法，需要在实践中根据具体情况选择最合适的方法,最小二乘估计最小一乘估计优点计算简单对模型假定稳健缺点对正态,1.5 多重线性回归,多个自变量 X1,Xp,dat = data.frame(y,x1,x2,x
14、3,x4);fit = lm(y.4, data=dat);fit.aic = step (fit);,变量过多：降低估计的效率,可能使得那些真正对因变量有贡献的 自变量的效应不显著，过度拟合也会降低预测的准确性(overfit).,变量选择：将真正有对因变量有贡献的自变 量选出来。准则 Akaike Information Criterion (AIC)准则；Bayesian Information Criterion (BIC)准则；筛选方法：逐步法,理论上可以选出任意阶交互效应，阶数小于等于自变量的个数。,自变量越多，拟合的残差越小,1.5 多重线性回归 多个自变量 X1,Xpdat =
15、,1.6 Box-Cox变换,通过残差诊断发现残差不太正态时，可以考虑进行Box-Cox变换（power transformation），使得因变量变换后进行回归分析的残差更正态，统计分析结果更可信。,require(MASS);bc = boxcox(weight., data=dat, lambda=seq(-1,2,0.1);lambda = bc$xwhich.max(bc$y);# install.packages(TeachingDemos);require(TeachingDemos);weight.bc = bct(dat$weight,lambda); ,用变换后的因变量代替
16、原有因变量,给出不同 lambda值对应的似然值,最大似然值对应的lambda,Box-Cox变换,1.6 Box-Cox变换通过残差诊断发现残差不太正态时，可,2. 相关分析,2.1 Pearson相关系数2.2 相关系数的统计推断2.3 偏相关系数2.4 秩相关、多重相关,2. 相关分析2.1 Pearson相关系数,X 与 Y 的Pearson 相关系数,Pearson相关系数与回归系数的关系,2.1 Pearson 相关系数,样本Pearson 相关系数,度量两个变量的相关程度,X 与 Y 的Pearson 相关系数Pearson相关系数,2.2 相关系数的统计推断,零假设：,检验统计
17、量(Fisher 变换)：,z 渐近正态, cortest = cor.test( dat$height,dat$weight, conf.level=0.95, alternative=two.sided); print(cortest$estimate); # rcor 0.8326746 print(cortest$statistic); # z t 5.421548,alternative可选项：”two.sided”, “less”, “greater”, print(cortest$p.value); # p-值1 0.0001167570 print(cortest$conf.i
18、nt); # 置信区间1 0.5587423 0.9427914attr(,conf.level)1 0.95,2.2 相关系数的统计推断零假设：检验统计量(Fisher,2.3 偏相关系数,控制分层等混杂因素的影响,R代码下载：yilab.gatech.edu/pcor.html,pcor.test(x, y, z),控制Z(可以是多维的)的影响下估计/检验X与Y的相关系数,2.3 偏相关系数控制分层等混杂因素的影响R代码下载：yil,2.4 秩相关、多重相关,秩相关系数：X与Y的Spearman秩(大小次序)的Pearson相关系数 不需要正态性假设对奇异值不敏感,Pearson相关系数的
19、统计推断基于正态性假设，且对奇异值敏感, cortest.spearman = cor.test(dat$height,dat$weight, method=spearman);警告信息：In cor.test.default(dat$height, dat$weight, method = spearman) : 无法给连结計算精確p值, print(cortest.spearman$estimate); rho 0.7791741 print(cortest.spearman$statistic); S 123.6625 print(cortest.spearman$p.value); 1
20、 0.0006173036,秩相关,2.4 秩相关、多重相关秩相关系数：X与Y的Spearman,mul.fit = lm(ca+b);cor.mul = cor.test(c, fitted(mul.fit);,多重相关,Y (一元的) 与 X (可以是多元的)的多重相关系数 (复相关系数)定义为因变量 Y 与其基于 X 的预测值 的相关系数，,X 是一元时退化为 X 与 Y 的相关系数乘以一个符号。,c与(a,b)的复相关系数及检验,X与Y是向量,mul.fit = lm(ca+b);多重相关 Y (一元,3 实例分析,Paruelo and Lauenroth (2019) invest
21、igated the geographic (latitude and longitude) and climatic (mean annual temperature, means annual precipitation and the proportion of the mean annual precipitation that fall in the periods June-August and December-February) patterns in the relative abundance of C3 plants throughout 73 sites across
22、North America.,自变量： MAP, MAT, JJAMAP, DJFMAP （气候）, LONG , LAT （地理）因变量：C3,例1,3 实例分析Paruelo and Lauenroth (2,大部分自变量基本服从正态分布，但是因变量正（右）偏,scatter plot+ box plot,C3与MAP及LAT相关性大,大部分自变量基本服从正态分布，但是因变量正（右）偏scatt,dat - read.csv(paruelo.csv);require(car);scatterplot.matrix(C3 + MAP + MAT + JJAMAP + DJFMAP + LON
23、G + LAT, data = dat, diag=boxplot);,Indescriptive statistics, abox plotorboxplot(also known as abox-and-whisker diagramorplot) is a convenient way of graphically depicting groups of numerical data through theirfive-number summaries. A boxplot may also indicate which observations, if any, might be co
24、nsideredoutliers. The spacings between the different parts of the box help indicate the degree ofdispersion(spread) andskewnessin the data.,图例,dat - read.csv(paruelo.csv),scatterplot.matrix(log10(C3 + 0.1) + MAP + MAT + JJAMAP + DJFMAP + LONG + LAT, data = dat, diag=boxplot);,C3可能为0，故加一个小的正数,scatter
25、plot.matrix(log10(C3 +,变换过后的因变量分布 更接近正态分布自变量之间有大的 相关性-共线性性,变换过后的因变量分布, print(cor(dat, 3:8); MAP MAT JJAMAP DJFMAP LONG LATMAP 1.0000000 0.355090766 0.11225905 -0.404512409 -0.73368703 -0.24650582MAT 0.3550908 1.000000000 -0.08077131 0.001478037 -0.21310910 -0.83859041JJAMAP 0.1122590 -0.080771307 1.
26、00000000 -0.791540381 -0.49155774 0.07417497DJFMAP -0.4045124 0.001478037 -0.79154038 1.000000000 0.77074399 -0.06512485LONG -0.7336870 -0.213109100 -0.49155774 0.770743994 1.00000000 0.09655281LAT -0.2465058 -0.838590413 0.07417497 -0.065124848 0.09655281 1.00000000, print(cor(dat, 3:8);,Paruelo &
27、Lauenroth (2019) separated the predictors into two groups for their analyses. One group included LAT and LONG and the other included MAP, MAT, JJAMAP and DJFMAP.,情形一：自变量只包含LAT, LONG模型1：C3 LAT + LONG 模型2：C3 LAT + LONG + LAT * LONG,情形二：自变量包含所有变量（变量选择）,Paruelo & Lauenroth (2019) sep,回归分析与相关分析精选课件,dat$c
28、LAT = scale(dat$LAT);dat$cLONG = scale(dat$LONG);fit1 = lm(log10(C3 + 0.1) cLAT + cLONG, data = dat);fit2 = lm(log10(C3 + 0.1) cLAT * cLONG, data = dat);require(MASS);par(mfrow=c(1,2);qqnorm(stdres(fit1),main=不含交互效应的残差QQ图,xlab=标准正态分布分位数, ylab=残差分位数);qqline(stdres(fit1), col = 2);qqnorm(stdres(fit2),
29、main=含交互效应的残差QQ图,xlab=标准正态分布分位数, ylab=残差分位数);qqline(stdres(fit2), col = 2);,以上为个性化作图，可以控制颜色、x轴与y轴等。也可以用 plot(fit1,which=2); plot(fit2,which=2)代替，但是无法控制作图参数。,dat$cLAT = scale(dat$LAT);以上为个,正态性基本成立方差齐次性不成立没有强影响点,par(mfrow=c(2,2);plot(fit2);,正态性基本成立par(mfrow=c(2,2);,为什么要用 scale（标准化：减去平均值后再除以标准差）？, fit0
30、=lm(log10(C3 + 0.1) LONG * LAT, data = dat); print(1/vif(fit0); # 容忍度：小于0.2认为共线性性大 LONG LAT LONG:LAT 0.014973575 0.003249445 0.002494144, fit2 = lm(log10(C3 + 0.1) cLAT * cLONG, data = dat); print(1/vif(fit2); # 容忍度：现在容忍度变大多了 cLAT cLONG cLAT:cLONG 0.8268942 0.9799097 0.8195915, tmp1 = data.frame(LAT
31、=dat$LAT,LONG=dat$LONG,LAT.LONG=dat$LAT*dat$LONG); tmp2 = data.frame(cLAT=dat$cLAT,cLONG=dat$cLONG,cLAT.cLONG=dat$cLAT*dat$cLONG); print(cor(tmp1); LAT LONG LAT.LONGLAT 1.00000000 0.09655281 0.9137741LONG 0.09655281 1.00000000 0.4894850LAT.LONG 0.91377412 0.48948498 1.0000000 print(cor(tmp2); cLAT c
32、LONG cLAT.cLONGcLAT 1.00000000 0.09655281 -0.4140077cLONG 0.09655281 1.00000000 -0.1344310cLAT.cLONG -0.41400770 -0.13443102 1.0000000,vif: variance inflation factor,为什么要用 scale（标准化：减去平均值后再除以标准差）？, print(summary(fit2)$coef); Estimate Std. Error t value Pr(|t|)(Intercept) -0.55294158 0.02746785 -20.1
33、304997 1.356751e-30cLAT 0.25664114 0.03025181 8.4834978 2.610338e-12cLONG -0.01659533 0.02778966 -0.5971763 5.523443e-01cLAT:cLONG 0.07686131 0.02988419 2.5719723 1.226881e-02, print(summary(fit2)$coef);,dat1 = dat,c(1,3:8); # 1: 因变量C3, 3-8: 协变量 dat1$C3 = log10(dat1$C3+0.1); fit3 = lm(C3., data=dat1
34、); aic.both = step(fit3,scale=0.1, direction=both, trace = FALSE); aic.backward = step (fit3,scale=0.1,direction=backward, trace = FALSE); aic.forward = step (fit3,scale=0.1,direction=forward, trace = FALSE);,情形二：自变量包含所有协变量（变量选择）,scale越大，选到的变量个数越多。,dat1 = dat,c(1,3:8); # 1, print(summary(aic.both)$c
35、oef); Estimate Std. Error t value Pr(|t|)(Intercept) -1.72472260 0.305824641 -5.639580 3.461079e-07JJAMAP -1.00198375 0.433027213 -2.313905 2.365729e-02DJFMAP -1.00534103 0.485632099 -2.070170 4.218135e-02LAT 0.04230919 0.005266712 8.033321 1.737289e-11, print(summary(aic.forward)$coef); Estimate St
36、d. Error t value Pr(|t|)(Intercept) -2.6886133013 1.2391496271 -2.1697245 3.363296e-02MAP 0.0002743328 0.0002175062 1.2612639 2.116537e-01MAT -0.0008468391 0.0116261123 -0.0728394 9.421542e-01JJAMAP -0.8338535725 0.4750796169 -1.7551870 8.386774e-02DJFMAP -0.9618361455 0.7163072633 -1.3427703 1.8394
37、63e-01LONG 0.0069235895 0.0100330182 0.6900804 4.925631e-01LAT 0.0434331291 0.0099284778 4.3746010 4.424396e-05, print(summary(aic.backward)$coef); Estimate Std. Error t value Pr(|t|)(Intercept) -1.72472260 0.305824641 -5.639580 3.461079e-07JJAMAP -1.00198375 0.433027213 -2.313905 2.365729e-02DJFMAP
38、 -1.00534103 0.485632099 -2.070170 4.218135e-02LAT 0.04230919 0.005266712 8.033321 1.737289e-11, print(summary(aic.both)$coef,Loyn (1987) selected 56 forest patches in southeastern Victoria, Australia, and related the abundance of forest birds in each patch (因变量) to six predictor variables: 1) patch
39、 area (ha), 2) distance to nearest patch (km), 3) distance to nearest larger patch (km), 4) grazing stock (1 to 5 indicating light to heavy), 5) altitude (m) and 6) years since isolation (years). Three of the predictor variables (patch area, distance to nearest patch or dist, distance to nearest lar
40、ger patch or ldist) were highly skewed, producing observations with high leverage, so these variables were transformed to log10. A correlation matrix indicated some moderate correlations between predictors, especially between log10(dist) and log10(ldist), log10(area) and graze, and graze and years.,
41、例 2,Loyn (1987) selected 56 forest,回归分析与相关分析精选课件, print(cor(dat,-1); YR.ISOL GRAZE ALT L10DIST L10LDIST L10AREAYR.ISOL 1.00000000 -0.63556710 0.2327154 -0.01957223 -0.16111611 0.2784145GRAZE -0.63556710 1.00000000 -0.4071671 -0.14263922 -0.03399082 -0.5590886ALT 0.23271541 -0.40716705 1.0000000 -0.2
42、1900701 -0.27404380 0.2751428L10DIST -0.01957223 -0.14263922 -0.2190070 1.00000000 0.60386637 0.3021666L10LDIST -0.16111611 -0.03399082 -0.2740438 0.60386637 1.00000000 0.3824795L10AREA 0.27841452 -0.55908864 0.2751428 0.30216662 0.38247952 1.0000000, print(1/vif(fit); # 容忍度 YR.ISOL GRAZE ALT L10DIS
43、T L10LDIST L10AREA 0.5540876 0.3960688 0.6812282 0.6043930 0.4975746 0.5231454, print(cor(dat,-1); print,回归分析与相关分析精选课件,负偏（左偏）：红色曲线左边多出一块。,负偏（左偏）：红色曲线,Box-Cox变换,require(MASS);dat1 = dat;lambda = 0.4;dat1$ABUND = bct( dat$ABUND, lambda); fit1=lm(ABUND ., data = dat1);par(mfrow=c(2,2);plot(fit1);,Box-C
44、ox变换require(MASS);, print(summary(fit1)$coef); Estimate Std. Error t value Pr(|t|)(Intercept) -187.99449889 99.88060702 -1.8821922 6.575550e-02YR.ISOL 0.10835940 0.04924116 2.2019860 3.250884e-02GRAZE -1.46916033 1.01297157 -1.4503471 1.533332e-01ALT 0.02519445 0.02609860 0.9653563 3.391068e-01L10DI
45、ST -0.79187819 2.91466964 -0.2716871 7.870040e-01L10LDIST -0.49433899 2.31225972 -0.2137904 8.315976e-01L10AREA 7.41795158 1.59571239 4.6486771 2.550921e-05, print(summary(fit1)$coef);,考虑2阶交互效应, dat2 = dat1; dat2,-1=apply(dat2,-1,2,scale); fit2=lm(ABUND (YR.ISOL+GRAZE+L10AREA)2, data = dat2); print(
46、1/vif(fit2); YR.ISOL GRAZE L10AREA YR.ISOL:GRAZE YR.ISOL:L10AREA GRAZE:L10AREA 0.3371914 0.3874801 0.5306881 0.3621688 0.4133117 0.5117608, print(c(summary(fit1)$r.squared,summary(fit2)$r.squared);1 0.6675484 0.7689699, print(summary(fit2)$coef); Estimate Std. Error t value Pr(|t|)(Intercept) 30.529
47、0605 1.173432 26.016907 2.442374e-30YR.ISOL -0.9326085 1.346090 -0.692828 4.916883e-01GRAZE -4.3720196 1.255705 -3.481723 1.057595e-03L10AREA 6.9622399 1.072982 6.488682 4.112168e-08YR.ISOL:GRAZE 3.8649862 1.345259 2.873043 5.994109e-03YR.ISOL:L10AREA 3.0593008 1.332670 2.295618 2.601664e-02GRAZE:L1
48、0AREA 4.4983000 1.067310 4.214613 1.071019e-04,考虑2阶交互效应 print(c(summary(fit,回归分析与相关分析精选课件, fit3=lm(ABUND (YR.ISOL+GRAZE+L10AREA)3, data = dat2); print(summary(fit3)$coef); Estimate Std. Error t value Pr(|t|)(Intercept) 30.5394183 1.189513 25.6738765 1.084511e-29YR.ISOL -0.9976824 1.493243 -0.6681314
49、 5.072499e-01GRAZE -4.3808491 1.271331 -3.4458762 1.192371e-03L10AREA 6.8879568 1.292601 5.3287573 2.608908e-06YR.ISOL:GRAZE 3.8684438 1.359437 2.8456213 6.499732e-03YR.ISOL:L10AREA 3.0487675 1.350021 2.2583105 2.850907e-02GRAZE:L10AREA 4.4561959 1.149741 3.8758269 3.221055e-04YR.ISOL:GRAZE:L10AREA
50、-0.1128350 1.069579 -0.1054948 9.164228e-01,3阶交互效应很不显著，不需要考虑,step(fit3)结果包含主效应和 2阶交互效应（fit2）, fit3=lm(ABUND (YR.ISOL+GRA, require(quantreg); fit4 = rq(ABUND (YR.ISOL+GRAZE+L10AREA)2, data = dat2); print(summary(fit4)$coef); coefficients lower bd upper bd(Intercept) 31.16729 29.52578 33.27361YR.ISOL
