GAMLSS代码示例

4.2 fitting regression models to continuous distributions

data("CD4")#HIV-1感染者母亲生下的小孩中未感染的小孩数量,age=孩子的年龄
plot(cd4~age,data=CD4)#数量足够多，因此我们可以视之为连续的响应变量
# 根据初步观察数据，有三个问题：1 很难直观地看出因变量和自变量之间的关系是否为线性
#2 响应变量方差的异质性是如何？3 给定age的前提下，cd4的分布是怎样的？是否服从正态分布？
# 以上问题通常的解决方法：对响应变量进行变换或对响应变量和解释变量都进行变换
op<-par(mfrow=c(3,4),mar=par("mar")+c(0,1,0,0),
        pch="+",cex=0.45,cex.lab=1.8,cex.axis=1.6)
page<-c("age^-0.5","log(age)","age^.5","age")#横坐标,对自变量age进行四种变换
pcd4<-c("cd4^-0.5","log(cd4+1)","cd4^.5")#纵坐标，对因变量cd4进行三种变换
for (i in 1:3){
  yy<-with(CD4,eval(parse(text=pcd4[i])))#字符串转命令行及执行操作
  for (j in 1:4){
    xx<- with(CD4,eval(parse(text=page[j])))
    plot(yy~xx,xlab=page[j],ylab=pcd4[i])
  }
}
par(op)
#以上画出的图中，也很难看出来线性，方差的均一性，正态误差的分布等
#GAMLSS framework可以一次性解决以上的难题
#下例中，用gamlss拟合一个方差为常数的，默认正态分布的响应变量
con<-gamlss.control(trace=FALSE)#gamlss fitting 的辅助函数，是否每一次迭代都显示
m1<-gamlss(cd4~age,sigma.fo=~1,data=CD4,control=con)#sigma即方差是固定值，均值不固定
m2<-gamlss(cd4~poly(age,2),sigma.fo=~1,data=CD4,control=con)#正交多项式，次数为2
m3<-gamlss(cd4~poly(age,3),sigma.fo=~1,data=CD4,control=con)
m4<-gamlss(cd4~poly(age,4),sigma.fo=~1,data=CD4,control=con)
m5<-gamlss(cd4~poly(age,5),sigma.fo=~1,data=CD4,control=con)
m6<-gamlss(cd4~poly(age,6),sigma.fo=~1,data=CD4,control=con)
m7<-gamlss(cd4~poly(age,7),sigma.fo=~1,data=CD4,control=con)
m8<-gamlss(cd4~poly(age,8),sigma.fo=~1,data=CD4,control=con)
GAIC(m1,m2,m3,m4,m5,m6,m7,m8)
GAIC(m1,m2,m3,m4,m5,m6,m7,m8,k=log(length(CD4$age)))#得到m7拟合效果最好
op<-par(mfrow=c(1,1))
plot(cd4~age,data=CD4)
lines(CD4$age[order(CD4$age)],fitted(m7)[order(CD4$age)],col="red")#画出m7的拟合曲线
#这个曲线画出的图看起来不是那么令人信服，曲线在尾部为了去适应数据摆动太大不够稳定，这是正交多项式拟合的典型特点。
#现在尝试另外一种方法，两参数，运用分维多项式和分段多项式，还有一个非参数的平滑方法：三次样条函数
m1f<-gamlss(cd4~fp(age,1),sigma.fo=~1,data=CD4,control=con)#fp()用来拟合可加的平滑项
m2f<-gamlss(cd4~fp(age,2),sigma.fo=~1,data=CD4,control=con)
m3f<-gamlss(cd4~fp(age,3),sigma.fo=~1,data=CD4,control=con)
GAIC(m1f,m2f,m3f)
GAIC(m1f,m2f,m3f,k=log(length(CD4$age)))
m3f#AIC和SBC选择了具有三项的多项式，mu=557.5
m3f$mu.coefSmo#coefficients power transformations
plot(cd4~age,data=CD4)
lines(CD4$age[order(CD4$age)],fitted(m1f)[order(CD4$age)],col="blue")
lines(CD4$age[order(CD4$age)],fitted(m2f)[order(CD4$age)],col="green")
lines(CD4$age[order(CD4$age)],fitted(m3f)[order(CD4$age)],col="red")
#该结果拟合的还是不够好
#接下来尝试分段多项式，用多种不同的自由度
m2b<-gamlss(cd4~bs(age),data=CD4,trace=FALSE)
m3b<-gamlss(cd4~bs(age,df=3),data=CD4,trace=FALSE)
m4b<-gamlss(cd4~bs(age,df=4),data=CD4,trace=FALSE)
m5b<-gamlss(cd4~bs(age,df=5),data=CD4,trace=FALSE)
m6b<-gamlss(cd4~bs(age,df=6),data=CD4,trace=FALSE)
m7b<-gamlss(cd4~bs(age,df=7),data=CD4,trace=FALSE)
m8b<-gamlss(cd4~bs(age,df=8),data=CD4,trace=FALSE)
GAIC(m2b,m3b,m4b,m5b,m6b,m7b,m8b)
GAIC(m2b,m3b,m4b,m5b,m6b,m7b,m8b,k=log(length(CD4$age)))
#画图，分段多项式，5,6,7个自由度下的图
plot(cd4~age,data=CD4)
lines(CD4$age[order(CD4$age)],fitted(m5b)[order(CD4$age)],col="blue")
lines(CD4$age[order(CD4$age)],fitted(m6b)[order(CD4$age)],col="green")
lines(CD4$age[order(CD4$age)],fitted(m7b)[order(CD4$age)],col="red")
#继续用三次样条来拟合数据，
fn<-function(p) AIC(gamlss(cd4~cs(age,df=p[1]),data=CD4,trace=FALSE),k=2)#惩罚因子k=2
opAIC<-optim(par=c(3),fn,method="L-BFGS-B",lower=c(1),upper=c(15))#用AIC选择最优模型
fn<-function(p) AIC(gamlss(cd4~cs(age,df=p[1]),data=CD4,trace=FALSE),k=log(length(CD4$age)))#惩罚因子k=log(n)
opSBC<-optim(par=c(3),fn,method="L-BFGS-B",lower=c(1),upper=c(15))#用SBC选择最优模型
opAIC$par#根据AIC，最优的模型是平滑自由度为10.85，约等于11
opSBC$par#根据SBC，最优的模型是平滑自由度为1.854，约等于2
maic<-gamlss(cd4~cs(age,df=10.85),data=CD4,trace=FALSE)#用10.85的自由度拟合
msbc<-gamlss(cd4~cs(age,df=1.85),data=CD4,trace=FALSE)#用1.85的自由度拟合
plot(cd4~age,data=CD4)
lines(CD4$age[order(CD4$age)],fitted(maic)[order(CD4$age)],col="blue")
lines(CD4$age[order(CD4$age)],fitted(msbc)[order(CD4$age)],col="green")
#接下来用pb()函数,Pb()函数自动选择的自由度mu，运用了局部最大似然程序，是8.606（包括常数和线性项的2）。这个数值介于AIC和SBC建议的值之间。
mpb1<-gamlss(cd4~pb(age),data=CD4,trace=FALSE)#pb()自动选择自由度
mpb1$mu.df
#现在运用pb()函数来找到既适合mu又适合log(sigma)的模型(因为sigma与正态分布的默认连接函数是log):
mpb2<-gamlss(cd4~pb(age),sigma.fo=~pb(age),data=CD4,gd.tol=10,trace=FALSE)#默认的容差为5，我们增加到10，以便模型能够持续运算直到收敛
mpb2$mu.df
mpb2$sigma.df
#现在mu的自由度为6.2，sigma的自由度为3.8，mu的自由度低于之前的拟合结果。
plot(cd4~age,data=CD4)
lines(CD4$age[order(CD4$age)],fitted(mpb1)[order(CD4$age)],col="blue")
lines(CD4$age[order(CD4$age)],fitted(mpb2)[order(CD4$age)],col="green")
#上图中展示了mu模型的两个拟合结果
#让我们想象一下，如果我们运用cs()函数，并且试图用GAIC来评估自由度，会是怎样的结果。
m1<-gamlss(cd4~cs(age,df=10),sigma.fo=~cs(age,df=2),data=CD4,trace=FALSE)
fn<-function(p) AIC(gamlss(cd4~cs(age,df=p[1]),sigma.fo=~cs(age,p[2]),data=CD4,trace=FALSE,start.from=m1), k=2)
opAIC<-optim(par=c(8,3),fn,method="L-BFGS-B",lower=c(1,1),upper=c(15,15))
opAIC$par
#得到的总有效自由度for Mu 和 sigma 分别是5.72和3.81（包含常数和线性项的2），这与pb()函数得到的6.2和3.8非常接近
#用估计的总有效自由度（4.55,2.93）来重新跑一遍SBC的代码，需要注意的是，optim的lower界限改变为c(0.1,0.1)
m42<-gamlss(cd4~cs(age,df=3.72),sigma.fo=~cs(age,df=1.81),data=CD4,trace=FALSE)
fitted.plot(m42,mpb2,x=CD4$age,line.type=TRUE)
#广义交叉验证偏差方程VGD（）提供了另外一种调试平滑项自由度的方法。该方法很适合大体量的数据集，
#我们运用一部分数据来拟合这个模型（训练），一部分来验证模型。我们接下来证明该方法如何使用。
#首先，数据被分为训练和验证数据集，六四开。其次，我们使用optim()函数来查询mu和σ的最优平滑自由度
set.seed(1234)
rSample6040<-sample(2,length(CD4$cd4),replace=T,prob=c(0.6,0.4))#将数据分成率定和验证两组
fn<-function(p) VGD(cd4~cs(age,df=p[1]),sigma.fo=~cs(age,df=p[2]),data=CD4,rand=rSample6040)
op<-optim(par=c(3,1),fn,method="L-BFGS-B",lower=c(1,1),upper=c(10,10))
op$par
wp(mpb2,xvar=CD4$age,ylim.worm=1.5)#展示了四种去趋势的正态分布QQ图
#残差图表明，正态分布并不适合，因为U型的残差图表示残差的正偏分布