02FZ横截面回归.pptx

1、横截面数据: 因变量为实数轴上的数量变量,回归,回顾简单回归,w=read.table(COfreewy.txt,header=T)#读入数据 a=lm(CO.,w)#利用3个自变量做线性回归 summary(a)#展示结果 b=step(a,direction=backward)#逐步回归 summary(b)#展示逐步回归结果 shapiro.test(b$res)#做残差的正态性检验 qqnorm(b$res);qqline(b$res)#做残差的QQ图.,对这个数据再作进一步探索. 首先把各个变量的数据用散点图表示(图2.2). 从这6个散点图可以看出, CO和Traffic似乎有些线

2、性关系, CO和Hour则有些类似于正弦曲线一样的关系, 而CO和Wind的关系就比较复杂, 很难用线性关系表示. 根据时间序列分析所用的谐波分析(可参看Chatfield, 2004, p126), 可以用有穷Fourier级数来代表时间序列$x_t:$,w=read.table(COfreewy.txt,header=T)#读入数据 attach(w) #把变量名字放入内存 par(mfrow=c(2,3) #建立6个图的摆放模式 plot(COTraffic);plot(COHour);plot(COWind) plot(TrafficHour);plot(WindHour);plot(

3、TrafficWind) par(mfrow=c(1,1) cor(cbind(CO,Traffic,Tsq=Traffic2,Tcub=Traffic3,Hour,Hsq=Hour2, Hcub=Hour3,Wind,Wsq=Wind2,Wcub=Wind3)#计算Pearson线性相关系数 a=lm(COTraffic+Wind+I(Wind2)+I(Wind3)+sin(2*pi/24)*Hour)+ cos(2*pi/24)*Hour)+sin(4*pi/24)*Hour)+cos(4*pi/24)*Hour) b=step(a) #逐步回归, 按照AIC选择变量 summary(b)

4、;anova(b);shapiro.test(b$res),b1=lm(COTraffic+Wind+I(Wind2)+ cos(2*pi/24)*Hour)+cos(4*pi/24)*Hour) summary(b1)#结果汇总 anova(b1)#方差分析表 shapiro.test(b1$res)#对残差的正态性检验 qqnorm(b1$res);qqline(b1$res)#做QQ图(请读者自己做, 这里不展示),讨论,在上面模型选择中既用了基于AIC的逐步回归, 又用了对系数的$t$检验及$F$检验的$p$值, 到底依照什么标准来选择变量呢? 根据不同模型对数据的解释也不同. 这是经

5、典回归的固有问题. 这正如在物理学中, 一个现象有多种假说来解释, 每一种假说都不是真理, 只不过是人们根据自己的标准从不同角度对客观世界的猜想. 另外, 那些显著性检验的依据是正态性分布, 因此考察正态性是必不可少的, 如果没有了正态性, 这些检验就没有多大道理了. 如果不想评价模型, (对于简单的一元回归)在散点图上信手画一条曲线, 也是完全是一种回归, 只很难说清楚你这条曲线比另一人画的曲线要优越而已. 这个数据中的自变量都做了不同的变换, 但这还是线性模型.,简单线性模型不易处理的横截面数据,u=read.csv(pharynx1.csv)#读入数据 x=1:11;(x=x-c(5,1

6、1)#定性变量的下标 for(i in x)u,i=factor(u,i)#把定性变量从数值型转换成因子型 a=lm(TIME.,data=u);summary(a) #变换 a=lm(TIME0.3 INST + SEX + TX + AGE + COND + T.STAGE + N.STAGE + STATUS, data = u);b=step(a) summary(b);anova(b) shapiro.test(b$res) par(mfrow=c(1,2) plot(b$res, main=Residuals);abline(h=0)#残差图 qqnorm(b$res);qqlin

7、e(b$res) #正态QQ图,hei 讨论. 即使一个回归的检验全部显著而且$R2$也很接近1也不能一定说明该回归就有意义. 请试着运行下面的语句, 并且查看各种输出: set.seed(44) x=c(rnorm(100),50);y=c(rnorm(100),-50) a=lm(yx);summary(a) shapiro.test(a$res),生存分析数据的Cox回归模型,library(survival) fit - survfit(Surv(TIME, as.numeric(STATUS) TX, data=u) plot(fit,lty=1:2, ylab=S(t),xlab=

8、t, main=Survival Functions) legend(1500,1,c(TX=1,TX=2),lty=1:2)#图例,library(survival) fit - coxph( Surv(TIME,as.numeric(STATUS).,data=u) #Cox回归模型 plot(survfit( fit) #拟合的生存函数 summary(fit)#回归结果,Cox比例危险模型的多重分数多项式模型(Multiple Fractional Polynomial Model),library(mfp) f=mfp(Surv(TIME,as.numeric(STATUS)fp(A

9、GE, df=4,select=0.05) +INST+SEX+TX+GRADE+COND+SITE+T.STAGE+N.STAGE,family= cox,data=u) print(f)#输出结果 (rsq=1-sum(f$residuals)2)/sum(u$TIME-mean(u$TIME)2)#R2,数据出现多重共线性情况: 岭回归, lasso回归, 偏最小二乘回归,有一些关于多重共线的度量, 其中之一是容忍度(tolerance)或(等价地)方差膨胀因子(variance inflation factor, 简写成VIF), 而另一个是条件数(Condition Number),

10、 常用$kappa$表示. 其中容忍度与VIF定义为,这里$R2_j$是第$j$个变量在所有其它变量上回归时的确定系数 容忍度太小(按照一些文献, 比如小于0.2或0.1)或VIF太大(比如大于5或10)则有多重共线性问题.,条件数,或曰: 当$kappa 15$则有共线性问题, 而$kappa30$时说明共线性的问题严重.,library(mfp) w=read.csv(diabetes.csv)#注:w,1:10为x,w,11为y,而w,12:75为x2 kappa(w,12:75)#x2的条件数 library(car)#包含vif的软件包 vif(lm(y.,w,11:75),岭回归,

11、OLS,岭回归,library(MASS) a-lm.ridge(y.,lambda=seq(0,150,length=151), data=w,11:75,model =TRUE) names(a)#变量名字 a$lambdawhich.min(a$GCV) #找到GCV最小时的lambdaGCV=81.1 a$coef,which.min(a$GCV) #找到GCV最小时对应的系数 par(mfrow=c(1,2) #画出图形，并作出lambdaGCV 取0.01 时的那条竖直线 matplot(a$lambda,t(a$coef),xlab=expression(lambda), yla

12、b=Coefficients,type=l,lty=1:20) abline(v=a$lambdawhich.min(a$GCV) #下面语句绘出lamda 同GCV 之间关系的图形 plot(a$lambda,a$GCV,type=l,xlab=expression(lambda), ylab=expression(beta) abline(v=a$lambdawhich.min(a$GCV),Lasso回归,library(lars) x=as.matrix(w,1:10);y=as.matrix(w,11);x2=as.matrix(w,12:75) laa=lars(x2,y) #la

13、rs函数只用于矩阵型数据 plot(laa) #绘出图2.5 summary(laa)#给出Cp值(表2.6) cva=cv.lars(x2,y,K=10) #10折交叉验证 best=cva$indexwhich.min(cva$cv)#选适合的值(随机性使得结果不同) coef=coef.lars(laa,mode=fraction,s=best)#使得CV最小步时的系数 names(laa$Cpwhich=min(laa$Cp)#给出17 coef1=coef.lars(laa,mode=step,s=17)#使laa$Cp最小的step时的系数,偏最小二乘回归,library(lars

14、) library(pls) ap=plsr(y x2, 64, validation = CV) #求出所有可能的64个因子 ap$loadings #看代表性, 前28个因子可以代表76.4%的方差 ap$coef #看各个因子作为原变量的线性组合的系数 validationplot(ap)#此图可用来根据CV所用准则最优的原则挑选因子数量 RMSEP(ap);MSEP(ap);R2(ap)#不同准则(MSEP,R2)在不同因子数量时的值,可以看出5个因子时RMSEP最小. 用其他准则, 比如MSEP及R2都选择了5个因子,无法做任何假定的数据: 机器学习回归方法,hei 例2.4 mg数

15、据(mg.csv). 这个数据来自Flake and Lawrence(2002), 可从网站下载footnote.tw/cjlin/libsvmtools/datasets/regression/mg. 只知道该数据有6个自变量和一个因变量, 一共有1385个观测值, 没有关于该数据的任何其他信息.,w=read.csv(mg.csv)#读入数据 a=lm(y.,w)#简单线性回归 cor(w)#相关系数表,下面我们对各种方法都用5折交叉验证的方法来判断其结果的可靠性. 计算中通过随机建立的五个训练集建立5个模型, 得到5个测试集的标准化均方误差(

16、NMSE)及均方误差(MSE), 再得出5次平均的NMSE及MSE.,随机选下标集的代码为,n=1385;zz1=1:n #zz1为所有观测值(行)的下标 zz2=rep(1:5,ceiling(1385/5)1:n set.seed(100);zz2=sample(zz2,n)#zz2为1:5的随机排列,#下面建立一些都是0的5元素向量以存取结果 NMSE=rep(0,z$K);MSE=NMSE;NMSE0=MSE;MSE0=MSE for(i in 1:5) #对每一组训练集和测试集做一次, 共5次 m=zz1zz2=i #m为测试集下标集合 a=lm(y.,data=w-m,) #简单线

17、性回归, 这里-m为训练集下标集合 y0=predict(a,w-m,) #对训练集预测 y1=predict(a,wm,) #对测试集预测 #训练集的MSE: NMSE0i=mean(w$y-m-y0)2)/mean(w$y-m-mean(w$y-m)2) #训练集的NMSE: MSE0i=mean(w$y-m-y0)2) #测试集的MSE: NMSEi=mean(w$ym-y1)2)/mean(w$ym-mean(w$ym)2) #测试集的NMSE: MSEi=mean(w$ym-y1)2) #下面输出训练集及测试集的平均MSE及NMSE: (MMSE0=mean(MSE0);(MNMSE

18、0=mean(NMSE0) (MMSE=mean(MSE);(MNMSE=mean(NMSE),简单回归,决策树回归(回归树),library(rpart.plot)#同时自动打开rpart a=rpart(y.,w);a #计算决策树并输出决策树的细节 rpart.plot(a,type=2,faclen=T) #画出决策树的图,library(rpart); NMSE=rep(0,5);MSE=NMSE;NMSE0=MSE;MSE0=MSE; for(i in 1:5) m=zz1zz2=i a=rpart(y.,w-m,) #决策树回归(回归树) y0=predict(a,w-m,) y

19、1=predict(a,wm,) NMSE0i=mean(w$y-m-y0)2)/mean(w$y-m-mean(w$y-m)2) MSE0i=mean(w$y-m-y0)2) NMSEi=mean(w$ym-y1)2)/mean(w$ym-mean(w$ym)2) MSEi=mean(w$ym-y1)2) (MMSE0=mean(MSE0);(MNMSE0=mean(NMSE0) (MMSE=mean(MSE);(MNMSE=mean(NMSE),Boosting回归,比如$f_j(x(j)$可以是由第$j$个自助法样本$x(j)$得到的决策树. 如果定义了损失函数$rho$, 则模型需要按

20、照下面公式拟合:,这里$n$为训练样本量, $w_i$为某种权重, $x_i$为高维自变量观测值, $y_i$为相应的因变量观测值.,Boosting回归,library(mboost) NMSE=rep(0,5);MSE=NMSE;NMSE0=MSE;MSE0=MSE; for(i in 1:5) m=zz1zz2=i B1=mboost(y btree(x1)+btree(x2)+btree(x3) +btree(x4)+btree(x5)+btree(x6),data =w-m,) y0=predict(B1,w-m,) y1=predict(B1,wm,) NMSE0i=mean(w$

21、y-m-y0)2)/mean(w$y-m-mean(w$y-m)2) MSE0i=mean(w$y-m-y0)2) NMSEi=mean(w$ym-y1)2)/mean(w$ym-mean(w$ym)2) MSEi=mean(w$ym-y1)2) (MMSE0=mean(MSE0);(MNMSE0=mean(NMSE0) (MMSE=mean(MSE);(MNMSE=mean(NMSE),随机森林回归,library(randomForest) set.seed(101) NMSE=rep(0,5);MSE=NMSE;NMSE0=MSE;MSE0=MSE; for(i in 1:5) m=zz

22、1zz2=i A=randomForest(y.,data=w-m,importance=TRUE,proximity=TRUE) y0=predict(A,w-m,) y1=predict(A,wm,) NMSE0i=mean(w$y-m-y0)2)/mean(w$y-m-mean(w$y-m)2) MSE0i=mean(w$y-m-y0)2) NMSEi=mean(w$ym-y1)2)/mean(w$ym-mean(w$ym)2) MSEi=mean(w$ym-y1)2) (MMSE0=mean(MSE0);(MNMSE0=mean(NMSE0) (MMSE=mean(MSE);(MNMSE=mean(NMSE),人工神经网络回归,library(nnet) set.seed(444) NMSE=rep(0,5);MSE=NMSE;NMSE0=