02FZ横截面回归

横截面数据:因变量为实数轴上的数量变量回归回顾简单回归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似乎有些线性关系,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(CO~Traffic);plot(CO~Hour);plot(CO~Wind)plot(Traffic~Hour);plot(Wind~Hour);plot(Traffic~Wind)par(mfrow=c(1,1))cor(cbind(CO,Traffic,Tsq=Traffic^2,Tcub=Traffic^3,Hour,Hsq=Hour^2,Hcub=Hour^3,Wind,Wsq=Wind^2,Wcub=Wind^3))#计算Pearson线性相关系数a=lm(CO~Traffic+Wind+I(Wind^2)+I(Wind^3)+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);anova(b);shapiro.test(b$res)b1=lm(CO~Traffic+Wind+I(Wind^2)+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$值,到底依照什么标准来选择变量呢?•根据不同模型对数据的解释也不同.这是经典回归的固有问题.这正如在物理学中,一个现象有多种假说来解释,每一种假说都不是真理,只不过是人们根据自己的标准从不同角度对客观世界的猜想.•另外,那些显著性检验的依据是正态性分布,因此考察正态性是必不可少的,如果没有了正态性,这些检验就没有多大道理了.•如果不想评价模型,(对于简单的一元回归)在散点图上信手画一条曲线,也是完全是一种回归,只很难说清楚你这条曲线比另一人画的曲线要优越而已.•这个数据中的自变量都做了不同的变换,但这还是线性模型.简单线性模型不易处理的横截面数据u=read.csv(pharynx1.csv)#读入数据x=1:11;(x=x[-c(5,11)])#定性变量的下标for(iinx)u[,i]=factor(u[,i])#把定性变量从数值型转换成因子型a=lm(TIME~.,data=u);summary(a)#变换a=lm(TIME^0.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);qqline(b$res)#正态QQ图{\hei讨论.}即使一个回归的检验全部显著而且$R^2$也很接近1也不能一定说明该回归就有意义.请试着运行下面的语句,并且查看各种输出:set.seed(44)x=c(rnorm(100),50);y=c(rnorm(100),-50)a=lm(y~x);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=t,main='SurvivalFunctions')legend(1500,1,c(TX=1,TX=2),lty=1:2)#图例0500100015000.00.20.40.60.81.0SurvivalFunctionstS(t)TX=1TX=2library(survival)fit-coxph(Surv(TIME,as.numeric(STATUS))~.,data=u)#Cox回归模型plot(survfit(fit))#拟合的生存函数summary(fit)#回归结果Cox比例危险模型的多重分数多项式模型(MultipleFractionalPolynomialModel)library(mfp)f=mfp(Surv(TIME,as.numeric(STATUS))~fp(AGE,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))#R^2数据出现多重共线性情况:岭回归,lasso回归,偏最小二乘回归有一些关于多重共线的度量,其中之一是容忍度(tolerance)或(等价地)方差膨胀因子(varianceinflationfactor,简写成VIF),而另一个是条件数(ConditionNumber),常用$\kappa$表示.其中容忍度与VIF定义为这里$R^2_j$是第$j$个变量在所有其它变量上回归时的确定系数容忍度太小(按照一些文献,比如小于0.2或0.1)或VIF太大(比如大于5或10)则有多重共线性问题.条件数或曰:当$\kappa15$则有共线性问题,而$\kappa30$时说明共线性的问题严重.library(mfp)w=read.csv(diabetes.csv)#注:w[,1:10]为x,w[,11]为y,而w[,12:75]为x2kappa(w[,12:75])#x2的条件数library(car)#包含vif的软件包vif(lm(y~.,w[,11:75]))岭回归OLS岭回归library(MASS)a-lm.ridge(y~.,lambda=seq(0,150,length=151),data=w[,11:75],model=TRUE)names(a)#变量名字a$lambda[which.min(a$GCV)]##找到GCV最小时的lambdaGCV=81.1a$coef[,which.min(a$GCV)]##找到GCV最小时对应的系数par(mfrow=c(1,2))#画出图形，并作出lambdaGCV取0.01时的那条竖直线matplot(a$lambda,t(a$coef),xlab=expression(lambda),ylab=Coefficients,type=l,lty=1:20)abline(v=a$lambda[which.min(a$GCV)])#下面语句绘出lamda同GCV之间关系的图形plot(a$lambda,a$GCV,type=l,xlab=expression(lambda),ylab=expression(beta))abline(v=a$lambda[which.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)#lars函数只用于矩阵型数据plot(laa)#绘出图2.5summary(laa)#给出Cp值(表2.6)cva=cv.lars(x2,y,K=10)#10折交叉验证best=cva$index[which.min(cva$cv)]#选适合的值(随机性使得结果不同)coef=coef.lars(laa,mode=fraction,s=best)#使得CV最小步时的系数names(laa$Cp[which=min(laa$Cp)])#给出17coef1=coef.lars(laa,mode=step,s=17)#使laa$Cp最小的step时的系数偏最小二乘回归library(lars)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.4mg数据(mg.csv).}这个数据来自FlakeandLawrence(2002),可从网站下载\footnote{~cjlin/libsvmtools/datasets/regression/mg.}.只知道该数据有6个自变量和一个因变量,一共有1385个观测值,没有关于该数据的任何其他信息.w=read.csv(mg.csv)#读入数据a=lm(y~.,w)#简单线性回归cor(w)#相关系数表下面我们对各种方法都用5折交叉验证的方法来判断其结果的可靠性.计算中通过随机建立的五个训练集建立5个模型,得到5个测试集的标准化均方误差(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=MSEfor(iin1:5){#对每一组训练集和测试集做一次,共5次m=zz1[zz2==i]#m为测试集下标集合a=lm(y~.,data=w[-m,])#简单线性回归,这里[-m]为训练集下标集合y0=predict(a,w[-m,])#对训练集预测y1=predict(a,w[m,])#对测试集预测#训练集的MSE:NMSE0[i]=mean((w$y[-m]-y0)^2)/mean((w$y[-m]-mean(w$y[-m]))^2)#训练集的NMSE:MSE0[i]=mean((w$y[-m]-y0)^2)#测试集的MSE:NMSE[i]=mean((w$y[m]-y1)^2)/mean((w$y[m]-mean(w$y[m]))^2)#测试集的NMSE:MSE[i]=