用R语言进行分位数回归Word文档格式.docx

用R语言进行分位数回归Word文档格式.docx

《用R语言进行分位数回归Word文档格式.docx》由会员分享，可在线阅读，更多相关《用R语言进行分位数回归Word文档格式.docx（26页珍藏版）》请在冰豆网上搜索。

98/100））

fit2=summary（rq（foodexp~income,tau=c,,,,））

windows（5,5） 

#新建一个图形窗口，可以去掉这句

plot（fit1）

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="

） 

#画中位数回归的拟合直线，颜色蓝

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%分位点的收入的消费估计值

#90%分位点的收入的消费估计值

windows（10,5）

par（mfrow=c（1,2）） 

#把绘图区域划分为一行两列

plot（c（ps,ps）,c,,type="

#type=”n”表示初始化图形区域，但不画图

xlab=expression（tau）,ylab="

quantile"

plot（stepfun（ps,c[1],）,=F,

add=T）

plot（stepfun（ps,c[1],）,=F,

add=T,="

="

=（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="

xlab="

Density"

lines,ar$dens,col="

lines,ap$dens,col="

black"

legend（"

topright"

c（"

poor"

"

rich"

）,lty=c（1,1）,

col=c（"

））

#比较不同分位点下，收入对食品支出的影响机制是否相同

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="

se="

ker"

#位置尺度漂移模型的检验

T3=KhmaladzeTest（y~X,taus=-1,nullH="

location-scale"

T4=KhmaladzeTest（y~X,taus=10:

##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="

cex=.25）#散点图

Dat<

-（x=x,y=y） 

#基本数据集

us<

-c（.25,.5,.75）

length（us））{

-vFrank（x,df,delta,u=us[i]）

lines（x,qt（v,df）） 

#v为概率，计算每个概率对应的T分布统计量

cfMat<

-matrix（0,3,length（us）+1）#初始矩阵，用于保存结果的系数

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="

#绘制预测曲线

cfMat[i,1]<

-tau 

#保存分位点的值

cfMat[i,2:

4]<

-coef（fit） 

#保存系数到cfMat矩阵的第i行

\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<

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）

#非参数的结果一般是通过画图查看的

plot（times,accel,xlab="

milliseconds"

ylab="

acceleration"

hs<

-c（1,2,3,4） 

#选择不同窗宽进行估计

for（iinhs）{

h=hs[i]

-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）"

Speedinlog（Km/hour）"

points（x[hoppers],y[hoppers],pch="

h"

col="

points（x[specials],y[specials],pch="

s"

others<

-（!

hoppers&

!

specials）

points（x[others],y[others],col="

cex=

-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）!

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="

main="

两组分位数回归结果比较"

points（engel1,cex=,col=2）

points（engel2,cex=,col=3）

for（iin1:

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）

）#将中的函数放入内存中

dat=（"

.dta"

dat=（dat） 

#去掉有缺失值的行

dim（dat） 

#这个数据包含10000个样本，3个指标

names（dat）

#[1]"

urtype"

"

q413"

q410a"

#q413每月家庭总收入

#q410a家庭食品支出

#urtype样本类型，1农村或2城镇

#这里考察上个月的收入与食品支出的关系，并分性别进行讨论

#把变量的名称改一下

names（dat）=c（"

income"

foodexp"

#进行log变换

dat2=dat[dat$income>

0&

dat$foodexp>

0,]

dat2[,"

lnincome"

]=log（dat2[,"

]）

lnfoodexp"

dat2=（dat2） 

#去掉log变换以后的缺失值行

dim（dat2）

#散点图，观察是否需要log变换

attach（dat）

#不变换

plot（income,foodexp,cex=,main="

散点图"

detach（dat）

#变换

attach（dat2）

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="

fit3=rq（lnfoodexp~lnincome,data=dat2[dat2$urtype==2,],

#所有数据放在一起是参数结果

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

Value

Std.Error

tvalue

Pr（>

|t|）

#参数结果的检验

T1=KhmaladzeTest（lnfoodexp~lnincome,taus=seq（.1,.9,by=.1）,

T2=KhmaladzeTest（lnfoo