1、课程安排课程安排1.回归模型与相关分析回归模型与相关分析2.方差分析方差分析3.属性数据分析属性数据分析4.生存数据分析生存数据分析5.实验设计实验设计6.全基因组关联分析全基因组关联分析 第第1,2,3,5次次1.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
2、 Murray LoganWiley-Blackwell,Chichester,West Sussex,2019.第第4次次3.Survival Analysis:A Self-Learning Text,2nd edition David G.Kleinbaum and Mitchel Klein,Springer,2019第第6次次4.ABEL tutorialYurii Aulchenko,2019参考书参考书回归分析和相关分析回归分析和相关分析生物统计学研究所生物统计学研究所 张洪张洪1.线性回归线性回归1.1 回归直线回归直线1.2 参数估计方法：最小二乘估计参数估计方法：最小二乘估
3、计1.3 参数推断：参数推断：F-检验、检验、t-检验、区间估计检验、区间估计1.4 回归模型的诊断：残差分析回归模型的诊断：残差分析1.5 多重线性回归多重线性回归1.6 Box-Cox变换变换2.相关分析相关分析2.1 Pearson相关系数相关系数 2.2 相关系数的统计推断：相关系数的统计推断：z-检验、区间估计检验、区间估计2.3 偏相关系数偏相关系数2.4 秩相关与多重相关秩相关与多重相关3 实例分析实例分析与函数关系的区别与函数关系的区别：同一身高可以体重不同，同一身高可以体重不同，体重随身高增加的关系体重随身高增加的关系不是严格成立，不是严格成立，只是有这种趋势只是有这种趋势总
4、的趋势总的趋势：随着身高增加，随着身高增加，体重也跟着增加体重也跟着增加为什么为什么?因为身高只能解释因为身高只能解释体重的一部分原因体重的一部分原因,还有其他未被考虑还有其他未被考虑的因素的因素,如饮食、地域、人种等如饮食、地域、人种等1.线性回归线性回归目标目标：建立一个：建立一个连续型连续型因变量因变量X(身高身高)与自变量与自变量Y(体重体重)之间的关系之间的关系eXY因变量因变量自变量自变量随机误差随机误差:未被身高解释的部未被身高解释的部分：分：饮食、种族、地饮食、种族、地域域截距截距斜率斜率:(1)0 正相关正相关 (2)fit=lm(weightheight,data=dat)
5、;print(fit);Call:lm(formula=weight height,data=dat)Coefficients:(Intercept)height -88.6774 0.8902(1)最小二乘估计最小二乘估计dat 是一个是一个data frame，有两个变量：，有两个变量：height 和和 weightR中线性模型拟合中线性模型拟合(2)最小一乘估计最小一乘估计估计方法估计方法R代码代码输出的结果输出的结果 print(fit$coef);Coefficients:(Intercept)x -88.6774 0.8902几个有用的函数：几个有用的函数：summary、nam
6、es1.3 参数推断参数推断预测值：预测值：iixy)()(yyyyyyiiiiniiiniiniiyyyyyy121212)()()(总平方和总平方和SStot自由度自由度=n-1回归平方和回归平方和SSreg自由度自由度=1残差平方和残差平方和SSerr自由度自由度=n-2回归平方和占总平方和的比重越大回归平方和占总平方和的比重越大,则则X能解释能解释Y的部分越大！的部分越大！回归方程越好！回归方程越好！toterrtotregSSSSSSSSR12 summary(fit)$r.squared1 0.693347niiyny11方差分解方差分解度量方法度量方法F检验检验)2/(1/nSS
7、SSFerrreg适用条件：误差适用条件：误差e1,en独立同分布，服从正态分布。独立同分布，服从正态分布。anova(fit)Analysis 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.自由度自由度平方和平方和平均平方平均平方FP-值值如果误差独立同分布如果误差独立同分布服从正态分布，则服从正态分布，则零假设成立时，零假设成立时，
8、F服从服从自由度为自由度为1和和n-2的的F分布分布方差分析表方差分析表回归部分回归部分残差部分残差部分0:0HF统计量统计量t 检验检验:回归系数的显著性检验回归系数的显著性检验n零假设零假设：回归系数等于：回归系数等于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 统计量统计量p-值值)(seT 零假设成
9、立且误差独立同分布零假设成立且误差独立同分布且服从正态分布时，且服从正态分布时，T 服从服从自由度为自由度为n-2的的t分布分布tn-2适用条件适用条件：n误差误差 e 独立同分布服从正态分布独立同分布服从正态分布n样本量不太小时，误差分布偏离正态分布样本量不太小时，误差分布偏离正态分布不是很大时仍适用，即对正态性假设不是很不是很大时仍适用，即对正态性假设不是很敏感敏感n当样本量大时总是适用的，因为此时当样本量大时总是适用的，因为此时 T 依依分布收敛于标准正态分布与分布收敛于标准正态分布与tn-2接近接近(不管正不管正态假设成立与否态假设成立与否)。区间估计区间估计)(),(2/1,22/1
10、,2setsetnn confint(fit)2.5%97.5%(Intercept)-149.8270467 -27.527805x 0.5354479 1.244863置信系数置信系数 的双侧置信区间，的双侧置信区间，常取为常取为0.051结论：结论：的的95%置信区间是置信区间是(0.5354479,1.244863)落在落在(0.5354479,1.244863)的概率是的概率是95%？错！？错！样本的两重性：抽样前是随机的，抽样后是非随机的。样本的两重性：抽样前是随机的，抽样后是非随机的。1.4 回归模型的诊断：残差分析回归模型的诊断：残差分析n 误差项独立同分布误差项独立同分布且服
11、从正态分布吗且服从正态分布吗?n 有没有需要剔除的有没有需要剔除的 强影响点强影响点?par(mfrow=c(2,2);plot(fit);满足下述条件则满足下述条件则 回归模型适用：回归模型适用：n残差没有明显的趋势残差没有明显的趋势n残差的残差的QQ图表明正态性成立图表明正态性成立nCook距离距离1 且且leverage1或或 leverage2p/n为强影响点为强影响点发现并剔除残差大的强影响点发现并剔除残差大的强影响点1)有些强影响点残差不大，不应被提出有些强影响点残差不大，不应被提出2)有些残差大的点不是强影响点，剔除与否关系不大有些残差大的点不是强影响点，剔除与否关系不大 dat
12、1=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)，由，由“帽帽子矩阵子矩阵”计算得到，其计算得到，其值越大对回归方程的影值越大对回归方程的影响力越大。响力越大。height weight lev1 180 68 0.10842272 160 51 0.249673011 165 57 0.135888012 18
13、5 76 0.203897514 164 58 0.154459816 190 200 0.3516872第第16个观测值杠杆很大个观测值杠杆很大稳健回归法：最小一乘估计稳健回归法：最小一乘估计(最小绝对偏差，最小绝对偏差，LAD)黑线为黑线为LAD回归直线回归直线 require(quantreg);fit1.lad=rq(weightheight,data=dat1);警告信息：警告信息：In rq.fit.br(x,y,tau=tau,.):Solution may be nonunique lines(dat1$height,fit1.lad$coef1+fit1.lad$coef2*
14、dat1$height,type=l,lwd=2,col=black);最小二乘估计最小二乘估计最小一乘估计最小一乘估计优点优点n计算简单计算简单n解唯一解唯一n最有效（模型假设成立时）最有效（模型假设成立时）n可给出方差分析表可给出方差分析表n对模型假定稳健对模型假定稳健n对奇异值不敏感对奇异值不敏感缺点缺点n对正态性假定高对正态性假定高n要求方差一致性要求方差一致性n对于奇异值敏感对于奇异值敏感n计算复杂度较高计算复杂度较高n解可能不唯一解可能不唯一n会损失一些效率会损失一些效率n目前的方差分析目前的方差分析不成熟不成熟没有一个方法可以完全打败其他的方法，没有一个方法可以完全打败其他的方法
15、，需要在实践中根据具体情况选择最合适的方法需要在实践中根据具体情况选择最合适的方法1.5 多重线性回归多重线性回归 多个自变量多个自变量 X1,XpeXXYpp11dat=data.frame(y,x1,x2,x3,x4);fit=lm(y.4,data=dat);fit.aic=step(fit);变量过多变量过多：降低估计的效率：降低估计的效率,可能使得那些真正对因变量有贡献的可能使得那些真正对因变量有贡献的 自变量的效应不显著，过度拟合也会降低预测的准确性自变量的效应不显著，过度拟合也会降低预测的准确性(overfit).变量选择变量选择：将真正有对因变量有贡献的自变：将真正有对因变量有
16、贡献的自变 量选出来。量选出来。准则准则 1)Akaike Information Criterion(AIC)准则；准则；2)Bayesian Information Criterion(BIC)准则；准则；筛选方法筛选方法：逐步法：逐步法理论上可以选出任意阶交互效应，理论上可以选出任意阶交互效应，阶数小于等于自变量的个数。阶数小于等于自变量的个数。eXXXXXXXXYpppppp1,13113211211.自变量越多，拟合的残差越小自变量越多，拟合的残差越小1.6 Box-Cox变换变换通过残差诊断发现通过残差诊断发现残差不太正态时残差不太正态时，可以考虑进行，可以考虑进行Box-Cox变
17、换变换（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);用变换后的因变量代替原有因变量用变换后的因变量代替原有因变量给出
18、不同给出不同 lambda值值对应的似然值对应的似然值最大似然值对应的最大似然值对应的lambda0),log(0,1)(xxxzBox-Cox变换变换2.相关分析相关分析2.1 Pearson相关系数相关系数2.2 相关系数的统计推断相关系数的统计推断2.3 偏相关系数偏相关系数2.4 秩相关、多重相关秩相关、多重相关X 与与 Y 的的Pearson 相关系数相关系数的样本方差的样本方差的样本协方差与YXYXr 不相关与负相关与正相关与YXYXYX,0,0,0Pearson相关系数与回归系数的关系相关系数与回归系数的关系的样本标准差的样本标准差的标准差的标准差YXrYX2.1 Pearson
19、 相关系数相关系数的方差的方差的协方差与YXYX样本样本Pearson 相关系数相关系数度量两个变量的相关程度度量两个变量的相关程度2.2 相关系数的统计推断相关系数的统计推断零假设：零假设：0检验统计量检验统计量(Fisher 变换变换)：rrz11log21z 渐近正态渐近正态*,11log21N 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
20、5.421548alternative可选项：可选项：”two.sided”,“less”,“greater”print(cortest$p.value);#p-值值1 0.0001167570 print(cortest$conf.int);#置信区间置信区间1 0.5587423 0.9427914attr(,conf.level)1 0.952.3 偏相关系数偏相关系数控制分层等混杂因素的影响控制分层等混杂因素的影响R代码下载：代码下载：yilab.gatech.edu/pcor.htmlpcor.test(x,y,z)pcor.test(x,y,z)控制控制Z(Z(可以是多维的可以是多
21、维的)的影响下的影响下估计估计/检验检验X X与与Y Y的的相关系数相关系数YYYXXXeZYeZX的相关系数与YXee2.4 秩相关、多重相关秩相关、多重相关秩相关系数：秩相关系数：X与与Y的的Spearman秩秩(大小次序大小次序)的的Pearson相关系数相关系数n 不需要正态性假设不需要正态性假设n对奇异值不敏感对奇异值不敏感Pearson相关系数的统计推断基于相关系数的统计推断基于正态性假设，正态性假设，且且对奇异值敏感对奇异值敏感 cortest.spearman=cor.test(dat$height,dat$weight,method=spearman);警告信息：警告信息：I
22、n 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 0.0006173036 秩相关秩相关mul.fit=lm(ca+b);cor.mul=cor.test(c,fitted(mul.fit);多重相关多重相关 Y(一元的一元的)与与 X(
23、可以是多元的可以是多元的)的的多重相关系数多重相关系数(复相关系数复相关系数)定义为因变量定义为因变量 Y 与其基于与其基于 X 的预测值的预测值 的相关系数，的相关系数，YX 是一元时退化为是一元时退化为 X 与与 Y 的相关系数乘以一个符号。的相关系数乘以一个符号。c与与(a,b)的复相关系数及检验的复相关系数及检验),(),(XYcorYYcorXY),()(),(),(YXcorsignXYcorYYcorXYX与与Y是向量是向量3 实例分析实例分析Paruelo and Lauenroth(2019)investigated the geographic(latitude and l
24、ongitude)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 North America.自变量：自变量：MAP,MAT,JJAMAP,DJFM
25、AP （气候）（气候）,LONG,LAT（地理）（地理）因变量：因变量：C3 例例1大部分自变量基本服从正态分布，但是因变量正（右）偏大部分自变量基本服从正态分布，但是因变量正（右）偏scatter plot+box plotC3与MAP及LAT相关性大dat 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
26、.21310910-0.83859041JJAMAP 0.1122590-0.080771307 1.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
27、.00000000Paruelo&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情形二：自变量包含所有变量（变量选择）情形二：自变量包含所有变量（变量选择）dat$cLAT=scale(da
28、t$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),main=含交互效应的残差含交互
29、效应的残差QQ图图,xlab=标准正态分布分位数标准正态分布分位数,ylab=残差分位数残差分位数);qqline(stdres(fit2),col=2);以上为个性化作图，可以控制颜色、以上为个性化作图，可以控制颜色、x轴与轴与y轴等。轴等。也可以用也可以用 plot(fit1,which=2);plot(fit2,which=2)代替，但是无法控代替，但是无法控制作图参数。制作图参数。正态性基本成立正态性基本成立方差齐次性不成立方差齐次性不成立没有强影响点没有强影响点par(mfrow=c(2,2);plot(fit2);为什么要用为什么要用 scale（标准化：减去平均值后再除以标准差）
30、？（标准化：减去平均值后再除以标准差）？fit0=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.81959
31、15 tmp1=data.frame(LAT=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(
32、cor(tmp2);cLAT cLONG cLAT.cLONGcLAT 1.00000000 0.09655281 -0.4140077cLONG 0.09655281 1.00000000 -0.1344310cLAT.cLONG-0.41400770-0.13443102 1.0000000vif:variance inflation factor print(summary(fit2)$coef);Estimate Std.Error t value Pr(|t|)(Intercept)-0.55294158 0.02746785-20.1304997 1.356751e-30cLAT
33、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 dat1=dat,c(1,3:8);#1:因变量因变量C3,3-8:协变量协变量 dat1$C3=log10(dat1$C3+0.1);fit3=lm(C3.,data=dat1);aic.both=step(fit3,scale=0.1,direction=both,trace=FALSE);
34、aic.backward=step(fit3,scale=0.1,direction=backward,trace=FALSE);aic.forward=step(fit3,scale=0.1,direction=forward,trace=FALSE);情形二：自变量包含所有协变量（变量选择）情形二：自变量包含所有协变量（变量选择）scale越大，选到的变量个数越多。越大，选到的变量个数越多。print(summary(aic.both)$coef);Estimate Std.Error t value Pr(|t|)(Intercept)-1.72472260 0.305824641-5.
35、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 Std.Error t value Pr(|t|)(Intercept)-2.6886133013 1.2391496271-2.1697245 3.363296e-02M
36、AP 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.839463e-01LONG 0.0069235895 0.0100330182 0.6900804 4.925631e-01LAT 0.0434331291 0.0099284778 4
37、.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 -1.00534103 0.485632099-2.070170 4.218135e-02LAT 0.04230919 0.005266712 8.033321 1.737289e-11Lo
38、yn(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 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)
39、altitude(m)and 6)years since isolation(years).Three of the predictor variables(patch area,distance to nearest patch or dist,distance to nearest larger patch or ldist)were highly skewed,producing observations with high leverage,so these variables were transformed to log10.A correlation matrix indicat
40、ed some moderate correlations between predictors,especially between log10(dist)and log10(ldist),log10(area)and graze,and graze and years.例例 2 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.000000
41、00-0.4071671-0.14263922-0.03399082-0.5590886ALT 0.23271541 -0.40716705 1.0000000-0.21900701-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.
42、30216662 0.38247952 1.0000000 print(1/vif(fit);#容忍度容忍度 YR.ISOL GRAZE ALT L10DIST L10LDIST L10AREA 0.5540876 0.3960688 0.6812282 0.6043930 0.4975746 0.5231454 负偏（左偏）：红色曲线负偏（左偏）：红色曲线左边多出一块。左边多出一块。Box-Cox变换变换require(MASS);dat1=dat;lambda=0.4;dat1$ABUND=bct(dat$ABUND,lambda);fit1=lm(ABUND .,data=dat1);p
43、ar(mfrow=c(2,2);plot(fit1);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-01L10
44、DIST -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考虑考虑2阶交互效应阶交互效应 dat2=dat1;dat2,-1=apply(dat2,-1,2,scale);fit2=lm(ABUND (YR.ISOL+GRAZE+L10AREA)2,data=dat2);print(1/vif(fit2);YR.ISOL GRAZE L10AR
45、EA 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.5290605 1.173432 26.016907 2.442374e-30Y
46、R.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:L10AREA 4.4983000 1.067310 4.214613 1.071
47、019e-04 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 5.072499e-01GRAZE -4.3808491 1.271331-3.4458762 1.192371e-03L10AREA 6.8879568 1.292601 5.3
48、287573 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-0.1128350 1.069579-0.1054948 9.164228e-013阶交互效应很不显著，不需要考虑阶交互效应很不显著，不需要考虑 step(fit3)结果包含主效应和结
49、果包含主效应和 2阶交互效应（阶交互效应（fit2）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 0.08424 -3.72091 0.74918GRAZE -4.58954 -6.94855 -3.46013 *L10AREA 6.23743 5.20394 8.37718 *YR.ISOL:GRAZE
50、3.98259 1.40819 6.97070 *YR.ISOL:L10AREA 2.42730 0.72038 4.06048 *GRAZE:L10AREA 4.42584 3.82665 5.36264 *print(summary(fit2)$coef);Estimate Std.Error t value Pr(|t|)(Intercept)30.5290605 1.173432 26.016907 2.442374e-30YR.ISOL -0.9326085 1.346090-0.692828 4.916883e-01GRAZE -4.3720196 1.255705-3.48172
