用R语言进行分位数回归.docx
- 文档编号:7412084
- 上传时间:2023-01-23
- 格式:DOCX
- 页数:26
- 大小:159.02KB
用R语言进行分位数回归.docx
《用R语言进行分位数回归.docx》由会员分享,可在线阅读,更多相关《用R语言进行分位数回归.docx(26页珍藏版)》请在冰豆网上搜索。
用R语言进行分位数回归
用R语言进行分位数回归:
基础篇
第一节 分位数回归介绍
[2][3]
[4]
[5][6]
[7]
[8]
[9]
[10]
[11]
第二节 用R语言进行分位数回归
(“quantreg”) #保持联网的情况下安装包
library(“quantreg”) #加载包
() #进入R帮助首页
help(rq) #获取rq函数的帮助,也可以写成:
rq
example(rq) #显示分位数回归函数rq()的一个简单示例代码
data(engel) #加载quantreg包自带的数据集,见说明①
fit1=rq(foodexp~income,tau=,data=engel) #进行分位数回归,见说明②
fit1 #直接显示分位数回归的模型和系数,见说明③
summary(fit1) #得到更加详细的显示结果,见说明④
r1=resid(fit1) #得到残差序列,并赋值为变量r1
c1=coef(fit1) #得到模型的系数,并赋值给变量c1,见说明⑤
summary(fit1,se=“nid”) #通过设置参数se,可以得到系数的假设检验,说明⑥
Frisch–Newton内点算法;C.“pfn”,针对特别大数据,使用经过预处理的Frisch–Newton逼近方法;D.“fnc”,针对被拟合系数特殊的线性不等式约束情况;E.“lasso”和“scad”,基于特定惩罚函数的平滑算法进行拟合。
#不同分位点下的系数估计值的比较
fit1=summary(rq(foodexp~income,tau=2:
98/100))
fit2=summary(rq(foodexp~income,tau=c,,,,))
windows(5,5) #新建一个图形窗口,可以去掉这句
plot(fit1)
windows(5,5) #新建一个图形窗口,可以去掉这句
plot(fit2)
#散点图
attach(engel) #打开engel数据集,直接运行其中的列名,就可以调用相应列
plot(income,foodexp,cex=,type="n", #画图,说明①
xlab="HouseholdIncome",ylab="FoodExpenditure")
points(income,foodexp,cex=,col="blue") #添加点,点的大小为
abline(rq(foodexp~income,tau=,col="blue") #画中位数回归的拟合直线,颜色蓝
abline(lm(foodexp~income),lty=2,col="red") #画普通最小二乘法拟合直线,颜色红
taus=c,,,,,
for(iin1:
length(taus)){ #绘制不同分位点下的拟合直线,颜色为灰色
abline(rq(foodexp~income,tau=taus[i]),col="gray")
}
detach(engel)
#比较穷人(收入在10%分位点的那个人)和富人(收入在90%分位点的那个人)的估计结果
#rq函数中,tau不在[0,1]时,表示按最细的分位点划分方式得到分位点序列
z=rq(foodexp~income,tau=-1)
z$sol #这里包含了每个分位点下的系数估计结果
=quantile(income,#10%分位点的收入
=quantile(income,#90%分位点的收入
ps=z$sol[1,] #每个分位点的tau值
=c(c(1,%*%z$sol[4:
5,]) #10%分位点的收入的消费估计值
=c(c(1,%*%z$sol[4:
5,]) #90%分位点的收入的消费估计值
windows(10,5)
par(mfrow=c(1,2)) #把绘图区域划分为一行两列
plot(c(ps,ps),c,,type="n", #type=”n”表示初始化图形区域,但不画图
xlab=expression(tau),ylab="quantile")
plot(stepfun(ps,c[1],),=F,
add=T)
plot(stepfun(ps,c[1],),=F,
add=T,="gray",="gray")
=(c(0,diff(ps))+c(diff(ps),0))/2
ap=akj,z=,p=
ar=akj,z=,p=
plot(c,,c(ap$dens,ar$dens),
type="n",xlab="FoodExpenditure",ylab="Density")
lines,ar$dens,col="gray")
lines,ap$dens,col="black")
legend("topright",c("poor","rich"),lty=c(1,1),
col=c("black","gray"))
#比较不同分位点下,收入对食品支出的影响机制是否相同
fit1=rq(foodexp~income,tau=
fit2=rq(foodexp~income,tau=
fit3=rq(foodexp~income,tau=
anova(fit1,fit2,fit3)
#残差形态的检验
source("C:
/ProgramFiles/R/")
x=gasprice
n=length(x)
p=5
X=cbind(x[(p-1):
(n-1)],x[(p-2):
(n-2)],x[(p-3):
(n-3)],x[(p-4):
(n-4)])
y=x[p:
n]
#位置漂移模型的检验
T1=KhmaladzeTest(y~X,taus=-1,nullH="location")
T2=KhmaladzeTest(y~X,taus=10:
290/300,
nullH="location",se="ker")
#位置尺度漂移模型的检验
T3=KhmaladzeTest(y~X,taus=-1,nullH="location-scale")
T4=KhmaladzeTest(y~X,taus=10:
290/300,
nullH="location-scale",se="ker")
##DemoofnonlinearquantileregressionmodelbasedonFrankcopula
vFrank<-function(x,df,delta,u) #某个非线性过程,得到的是[0,1]的值
-log(1-(1-exp(-delta))/(1+exp(-delta*pt(x,df))*((1/u)-1)))/delta
#非线性模型
FrankModel<-function(x,delta,mu,sigma,df,tau){
z<-qt(vFrank(x,df,delta,u=tau),df)
mu+sigma*z
}
n<-200 #样本量
df<-8 #自由度
delta<-8#初始参数
(1989)
x<-sort(rt(n,df)) #生成基于T分布的随机数
v<-vFrank(x,df,delta,u=runif(n)) #基于x生成理论上的非参数对应值
y<-qt(v,df) #v对应的T分布统计量
windows(5,5)
plot(x,y,pch="o",col="blue",cex=.25)#散点图
Dat<-(x=x,y=y) #基本数据集
us<-c(.25,.5,.75)
for(iin1:
length(us)){
v<-vFrank(x,df,delta,u=us[i])
lines(x,qt(v,df)) #v为概率,计算每个概率对应的T分布统计量
}
cfMat<-matrix(0,3,length(us)+1)#初始矩阵,用于保存结果的系数
for(iin1:
length(us)){
tau<-us[i]
cat("tau=",format(tau),"..")
fit<-nlrq(y~FrankModel(x,delta,mu,sigma,df=8,tau=tau),#非参数模型
data=Dat,tau=tau, #data表明数据集,tau分位数回归的分位点
start=list(delta=5,mu=0,sigma=1),#初始值
trace=T)#每次运行后是否把结果显示出来
lines(x,predict(fit,newdata=x),lty=2,col="red") #绘制预测曲线
cfMat[i,1]<-tau #保存分位点的值
cfMat[i,2:
4]<-coef(fit) #保存系数到cfMat矩阵的第i行
cat("\n") #如果前面把每步的结果显示出来,则每次的结果之间添加换行符
}
colnames(cfMat)<-c("分位点",names(coef(fit))) #给保存系数的矩阵添加列名
cfMat
#2-7-1半参数模型----
fit<-rq(y~bs(x,df=5)+z,tau=.33)
#其中bs()表示按b-spline的非参数拟合
#2-7-2非参数方法
lprq<-function(x,y,h,m=50,tau={
#这是自定义的一个非参数计算函数,在其他数据下同样可以使用
xx<-seq(min(x),max(x),length=m) #m个监测点
fv<-xx
dv<-xx
for(iin1:
length(xx)){
z<-x-xx[i]
wx<-dnorm(z/h) #核函数为正态分布,dnorm计算标准正态分布的密度值
r<-rq(y~z,weights=wx,tau=tau, #上面计算得到的密度值为权重
ci=FALSE)
fv[i]<-r$coef[1]
dv[i]<-r$coef[2]
}
list(xx=xx,fv=fv,dv=dv) #输出结果
}
library(MASS)
data(mcycle)
attach(mcycle)
windows(5,5) #非参数的结果一般是通过画图查看的
plot(times,accel,xlab="milliseconds",ylab="acceleration")
hs<-c(1,2,3,4) #选择不同窗宽进行估计
for(iinhs){
h=hs[i]
fit<-lprq(times,accel,h=h,tau= #关键拟合函数
lines(fit$xx,fit$fv,lty=i)
}
legend(45,-70,c("h=1","h=2","h=3","h=4"),
lty=1:
length(hs))
#2-7-3非参数回归的另一个方法----
#考察最大的跑步速度与体重的关系
data(Mammals)
attach(Mammals)
x<-log(weight) #取得自变量的值
y<-log(speed) #取得因变量的值
plot(x,y,xlab="Weightinlog(Kg)",ylab="Speedinlog(Km/hour)",
type="n")
points(x[hoppers],y[hoppers],pch="h",col="red")
points(x[specials],y[specials],pch="s",col="blue")
others<-(!
hoppers&!
specials)
points(x[others],y[others],col="black",cex=
fit<-rqss(y~qss(x,lambda=1),tau= #关键的拟合步骤
plot(fit,add=T) #add=T表示在上图中添加这里的曲线
#MM2005分位数分解的函数,改变参数可直接使用
MM2005=function(formu,taus,data,group,pic=F){
#furmu为方程,如foodexp~income
#taus为不同的分位数
#data总的数据集
#group分组指标,是一个向量,用于按行区分data
#pic是否画图,如果分位数比较多,建议不画图
engel1=data[group==1,]
engel2=data[group==2,]
#开始进行分解
fita=summary(rq(formu,tau=taus,data=engel1))
fitb=summary(rq(formu,tau=taus,data=engel2))
tab=matrix(0,length(taus),4)
colnames(tab)=c("分位数","总差异","回报影响","变量影响")
rownames(tab)=rep("",dim(tab)[1])
for(iin1:
length(taus)){
ya=cbind(1,engel1[,names(engel1)!
=formu[[2]]])%*%fita[[i]]$coef[,1]
yb=cbind(1,engel2[,names(engel2)!
=formu[[2]]])%*%fitb[[i]]$coef[,1]
#这里以group==1为基准模型,用group==2的数据计算反常规模型拟合值
ystar=cbind(1,engel2[,names(engel2)!
=formu[[2]]])%*%fita[[i]]$coef[,1]
ya=mean(ya)
yb=mean(yb)
ystar=mean(ystar)
tab[i,1]=fita[[i]]$tau
tab[i,2]=yb-ya
tab[i,3]=yb-ystar#回报影响,数据相同,模型不同:
模型机制的不同所产生的差异
tab[i,4]=ystar-ya#变量影响,数据不同,模型相同:
样本点不同产生的差异
}
#画图
if(pic){
attach(engel)
windows(5,5)
plot(income,foodexp,cex=,type="n",main="两组分位数回归结果比较")
points(engel1,cex=,col=2)
points(engel2,cex=,col=3)
for(iin1:
length(taus)){
abline(fita[[i]],col=2)
abline(fitb[[i]],col=3)
}
detach(engel)
}
#输出结果
tab
}
#下面是用一个样本数据做的测试
data(engel)
group=c(rep(1,100),rep(2,135)) #取前100个为第一组,后135个第二组
taus=c,,,, #需要考察的不同分位点
MM2005(foodexp~income,taus,data=engel,group=group,pic=F)#参数说明,见①
第三节 一个案例
#一个案例
setwd("F:
/实践2/2012_分位数回归等两个任务/分位数回归方法夹")#设置工作目录
setwd(””)
rm(list=ls())
gc()
library(quantreg)
library(foreign)
source("")#将中的函数放入内存中
dat=(".dta")
dat=(dat) #去掉有缺失值的行
dim(dat)
#这个数据包含10000个样本,3个指标
names(dat)
#[1]"urtype""q413" "q410a"
#q413每月家庭总收入
#q410a家庭食品支出
#urtype样本类型,1农村或2城镇
#这里考察上个月的收入与食品支出的关系,并分性别进行讨论
#把变量的名称改一下
names(dat)=c("urtype","income","foodexp")
#进行log变换
dat2=dat[dat$income>0&dat$foodexp>0,]
dat2[,"lnincome"]=log(dat2[,"income"])
dat2[,"lnfoodexp"]=log(dat2[,"foodexp"])
dat2=(dat2) #去掉log变换以后的缺失值行
dim(dat2)
#散点图,观察是否需要log变换
attach(dat)
windows(5,5)
#不变换
plot(income,foodexp,cex=,main="散点图")
detach(dat)
#变换
attach(dat2)
windows(5,5)
plot(lnincome,lnfoodexp,cexp=,main="散点图(对数变换以后)")
detach(dat2)
#建立基础的分位数回归模型
taus=c,,,,
fit1=rq(lnfoodexp~lnincome,data=dat2,tau=taus,method="fn")
fit2=rq(lnfoodexp~lnincome,data=dat2[dat2$urtype==1,],
tau=taus,method="fn")
fit3=rq(lnfoodexp~lnincome,data=dat2[dat2$urtype==2,],
tau=taus,method="fn")
#所有数据放在一起是参数结果
s1=summary(fit1)
#男性
s2=summary(fit2)
#女性
s3=summary(fit3)
#参数结果的直接比较
tabs(s1) #这里的tabs函数是笔者自己写的,方便显示结果
tabs(s2)
tabs(s3)
#具体的系数估计值及假设检验
tabcoef(s1) #这里的tabcoef函数是笔者自己写的,方便显示结果
tabcoef(s2)
tabcoef(s3)
Quantiles
(Intercept)
lnincome
Quantiles
(Intercept)
lnincome
Quantiles
(Intercept)
lnincome
Quantiles
Value
Std.Error
tvalue
Pr(>|t|)
Quantiles
Value
Std.Error
tvalue
Pr(>|t|)
Quantiles
Value
Std.Error
tvalue
Pr(>|t|)
#参数结果的检验
attach(dat2)
T1=KhmaladzeTest(lnfoodexp~lnincome,taus=seq(.1,.9,by=.1),
nullH="location",se="ker")
T2=KhmaladzeTest(lnfoo
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 语言 进行 位数 回归
![提示](https://static.bdocx.com/images/bang_tan.gif)