连续型因变量

1、鞠咪、刘素、白昊龙中国人民大学统计学院教材：孟生旺，回归模型，中国人民大学出版社，2015主要内容主要内容广义线性模型假设因变量服从指数分布族。在指数分布族中，既有连续分布，也有离散分布。本章讨论因变量服从连续分布的广义线性模型。指数分布族中的连续分布主要包括正态分布、伽马分布和逆高斯分布，相应地,与它们对应的广义线性模型称作正态回归模型、伽马回归模型和逆高斯回归模型.3.1 正态回归模型正态回归模型在正态分布假设下,广义线性模型等价于普通的线性回归模型。在广义线性模型的框架下讨论线性回归模型。模型设定模型设定 正态分布也称作高斯分布，是定义在 上的连续分布。如果因变量服从正态分布，则其密度函

2、数可以表示为： （3.1）式中， 是均值参数；是标准差参数， 也称为尺度参数。指数分布族的密度函数可以表示为下述一般形式： （3.2）,22222)(exp21),;(yyf);()(exp),;(ycwbyyfu对式（3.1）变形，正态分布的密度函数也可以表示为指数分布族的形式，即 )2ln(2122/exp)2ln(212)(exp2)(exp21),;(222222222222yuyuuyyyf对比式（3.2）和（3.3）可知 （3.4） 在广义线性模型中，正则连接函数是使得 成立的函数。在正态分布假设下， ，所以相应的正则连接函数就是恒等连接函数，即 。2/)(22b)(gug)(由式

3、（3.4）可得正态分布的均值为： （3.5）正态分布的方差函数为： （3.6）广义线性模型的残差偏差（简称为偏差）定义为： （3.7）式中， 表示当前模型的对数似然函数； 表示饱和模型的对数似然函数。在饱和模型的对数似然函数中，用观察值 代替了当前模型的拟合值 。bbb)(1)()()( bbv);,();,(222yyyD);,(2y);,(2yyy在正态回归模型中，当前模型的对数似然函数为： （3.8）饱和模型的对数似然函数为： （3.9）niiiiiyuuyy1222222)2ln(2122/);,(niniiiiyyyyy1212222222)2ln(21)2ln(2122/);,(因

4、此，正态回归模型的偏差为： （3.10）21122221222221222)(222)2ln(2122/)2ln(212);,();,(2iininiiiiiniiiiiniuyuuyyyuuyyyyD可见，正态回归模型的残差偏差就是线性回归中的残差平方和。换言之，最小化残差平方和等价于最小化偏差。广义线性模型的统计量定义为： （3.11）由于正态分布的方差函数等于1，所以，正态回归模型的统计量与偏差相同，即niiiiuvuy122)()(Duyniii1221)(3.1.2 迭代加权最小二乘估计迭代加权最小二乘估计在广义线性模型中，参数的迭代加权最小二乘估计公式为： （3.12）式中 （3.

5、13） （3.14） )1()1()()1(mmTmmTzWXXWX)(2)()(nniiigvwdiagW1)()(niiiigyz正则连接函数是使得成立的连接函数。在正态分布假设下，式（3.4）表明， ，所以正态回归模型的正则连接函数就是恒等连接函数，即 。在正态分布假设下，如果使用正则连接函数，则有 ，故式（3.13）和（3.14）可以简化为：相应地，式（3.12）可以简化为： )(g1)(gyzdiagW21yXXXTT1)(3.1.3 其他连接函数其他连接函数对数连接函数下的正态回归模型适用于对大于零的连续因变量建模。在广义线性模型提出之前，对于大于零的因变量，通常进行对数变换，然后

6、以因变量的对数为新的因变量建立普通的线性回归模型（即本节所谓的正态回归模型）。对数连接函数下的正态回归模型是对模型的均值进行对数变换，并把变换后的均值表示为线性预测项，即 。对数连接函数下的正态回归模型可以用极大似然法进行估计，但需要用对数连接函数下的均值 代替恒等连接函数下的均值。Tiix)ln()exp(Tix由式（3.3）可知，在对数连接函数下正态回归模型的对数似然函数可以写成：在对数连接函数下，正态回归模型的参数也可以使用迭代加权最小二乘法进行估计，此时，均值参数可以表示为 在正态分布假设下，还可以使用其他连接函数，如倒数连接函数，即令 倒数连接函数可以用于比例型因变量的建模,根据实际

7、需要，还可以使用更加一般的幂连接函数，即令 。niTiTiTiiixxxyy1222222)2ln(212)exp(2/)exp()exp();,()exp(TiixTiiixg /1)( 。Tipiixg)(3.1.4 模拟数据分析模拟数据分析为了说明正态回归模型可能存在的问题，本节首先用伽马分布模拟右偏的因变量，然后用正态回归模型进行拟合。假设因变量y服从伽马分布，受4个解释变量的影响，其中， 是连续型解释变量， 是分类解释变量。模拟数据时使用了6个参数，最后一个参数是和的交互效应，即模拟因变量的均值为：模拟的因变量服从均值为 ，标准差为 的伽马分布。21,xx43,xx)25. 02 .

8、 05 . 08 . 08 . 07exp(324321xxxxxx5.0 本例的因变量实际上服从右偏的伽马分布，模型模型mod1mod1的因变量是y，使用正态分布假设和恒等连接函数使用正态分布假设和恒等连接函数，等价于普通线性回归模型。该模型与模拟数据的机理相去甚远，拟合值与观察值之间的差异很大，出现了负的拟合值，如图所示。 模型模型mod2mod2的因变量是log（y），使用了正态分布使用了正态分布假设和恒等连接函数假设和恒等连接函数，相当对因变量行对数变换以后建立普通线性回归模型。该模型的参数估计值与模拟数据的真实值比较接近，模型的拟合值与观察值也比较接近，如图所示。模型mod2的结果表

9、明，虽然因变量实际上服从伽马分布，但对其做对数变换以后将比较接近正态分布，所以基于log（y）的普通线性回归模型的参数估计值更加接近参数的真实值。 模型模型mod3mod3的因变量是y，使用了正态分布假设和对使用了正态分布假设和对数连接函数数连接函数，属于广义线性模型之一该模型的拟合值与观察值比较接近，如图所示，但参数估计值与真实参数相差较大，这是因为该模型使用了错误的分布假设。绘图比较观察值与模型的拟合值绘图比较观察值与模型的拟合值3.2 伽马回归模型伽马回归模型常见的伽马分布有两个参数：形状参数 和尺度参数 ，其密度函数可以表示为： （3.17）上述伽马分布的均值为 ，方差为 。广义线性模

10、型是对均值参数建立回归模型，所以需要把上述伽马分布的均值设定为一个参数。在式（3.17）中，若令， ，则伽马分布的均值可以表示为 ，方差可以表示为 ，相应地，伽马分布的密度函数可以表示为： （3.18）/1)(1),;(yeyyf2/12)exp()/1 (1),;(/1yyyyf式中， 是伽马分布的均值参数； 称作离散参数。在均值给定的情况下，离散参数越大，伽马分布的离散程度越大，如图3-2所示，三个密度函数的均值都是100，离散参数分别为0.3,0.5和0.8.将式（3.18）的密度函数变形，可以将其表示为一般形式的指数分布族的密度函数： （3.19）上式表明，指数分布族的参数与伽马分布的

11、参数具有如下关系： （3.20）指数分布族的均值与方差是的函数。由式（3.20）可知 （3.21）所以，伽马分布的均值为 ，方差为 。)1(lnlnln1ln/exp),;(yyyf/1)ln()(b22)()()(11)( bbbb)(b2容易求得伽马分布的对数似然函数，从而可以求得伽马回归模型的残差偏差为： （3.22）niiiiiiniiiiiyyyyyyyD11)ln(2)ln(/ln12);();(23.2.2 迭代加权最小二乘估计迭代加权最小二乘估计对大于零的因变量建立伽马回归模型时，较常使用对数连接函数，即 ，故有 。由式（3.20）和（3.21）可知，伽马分布的方差函数为 ，且

12、有 ，将它们代入式（3.13）可得 （3.23） （3.24）离散参数的取值不影响回归参数的估计值，可以令式（3.23）中的 ，此时，对数连接函数下估计伽马回归参数的迭代加权最小二乘算法如表3-1所示。其中均值的初始值设定为 ，可以避免在迭代运算中分母上出现零值。)ln()(g/1)(g2v)(annnniiidiggvwdiagW/1 )()()(211/ )()()(niiiiniiiiygyz12/)(ymeany3.2.3 模拟数据分析模拟数据分析本节使用前面的模拟数据建立伽马回归模型，可以分别应用选代加权最小二乘法和glm函数对伽马回归模型中的参数进行估计。 两者的估计结果完全相同。

13、伽马回归模型与模拟数据使用的两者的估计结果完全相同。伽马回归模型与模拟数据使用的模型完全相同，所以参数估计值与真实值非常接近。模型完全相同，所以参数估计值与真实值非常接近。输出残差图输出残差图 ( 27 )3.3 逆高斯回归模型逆高斯回归模型l适用于对大于零的因变量建模l与伽马分布相比，逆高斯分布具有尖峰厚尾特征( 28 )代码代码l在均值相等（均为100）、方差相等（均为2500）的条件下对伽马分布（GA）和逆高斯分布（IA）的密度函数进行了比较( 29 )模型设定模型设定l常见的逆高斯分布的密度函数如下：l表述为下述指数分布族的形式：2232222111(y; ,)expln 222yfy

14、y 222321()( ; ,)exp2()2yf yyy ( 30 )逆高斯分布逆高斯分布l与指数分布的一般形式对比：l逆高斯分布的均值及方差函数： 均值是 ，方差为221( )12b 3231( )()()( )( )bbbbb23 ( 31 )残差偏差残差偏差l 逆高斯分布假设下的对数似然函数：l逆高斯回归模型的残差偏差：232221/(2) 1/11ln(2)22niiiiiylyy 222222221212212 ( ,; )( ,; )/(2) 1/(2) 1/221()niiiiiiiniiiiiniiiiiDl yylyyyyyyyyy ( 32 )3.3.2 迭代加权最小二乘

15、估计迭代加权最小二乘估计l逆高斯回归模型通常使用对数连接函数：l迭代加权最小二乘算法：( )ln( )g22()(1)(1)1/()()()()()()/iin niin niiiiiiiinnWdiagdiagvgzygy ( 33 )迭代加权最小二乘估计算法迭代加权最小二乘估计算法( 34 )模拟数据分析模拟数据分析l 生成模拟数据：( 35 )参数估计参数估计lGlm函数：l迭代加权最小二乘估计：( 36 )参数估计输出参数估计输出lGlm函数：l迭代加权最小二乘估计：( 37 )残差分析图残差分析图( 38 )3.4 基于基于R的应用的应用在R中有两个函数可以建立广义线性模型：l基础包

16、中的glm函数：专门用于建立广义线性模型lgamlss程序包中的gamlss函数：可以建立更加一般意义上的回归模型( 39 )补充：补充：Generalized Additive Models for Location, Scale and ShapelRigby和Stasinopoulos（2005）提出了基于位置、尺度、形状参数的广义可加模型（gamlss），将分布的假设形式由指数分布族推广到了更加一般的形式， 包括一大类常见的分布类型， 而且可以在各种分布假设下同时对一个分布的位置参数、尺度参数和形状参数建立参数或非参数的回归模型， 具有很大的灵活法( 40 )补充：补充：General

17、ized Additive Models for Location, Scale and Shape( 41 )补充：补充：Generalized Additive Models for Location, Scale and Shape( 42 )补充：补充：GAMLSS模型的优点l一是极大地扩展了分布的种类，提供多种方式产生不同的分布。数据中的不规则特点，如高峰性、厚尾、右偏性等，得到了更多重视。l二是系统性成分有了更为丰富的内容，可以引入更为复杂的半参数或非参数部分及随机效应部分。不但可以建立均值与解释变量间的模型，而且可以建立分布的其他参数与解释变量间的模型。如果保费原理不仅与分布的位

18、置参数有关，还与分布的尺度参数及形状参数有关，那GAMLSS 给出的费率表将是高度风险敏感的，这对风险识别有重要意义。( 43 ) 补充：补充：GAMLSS模型模型lGLM模型：把分布假设从正态分布推广到了指数分布族，并引入了连接函数。lGAMLSS模型：不仅可以对分布的均值参数建立回归模型，还可以对分布的离散参数建立回归模型,1,2,()iiTiiyingx均值为 的指数分布族( ;,)(),1,2,( )iiiTiiTiiyf ygxinhx ( 44 )补充：补充：Gamlss包及包及gamlss函数函数lgamlss程序包中的gamlss函数：( 45 )补充：补充：Gamlss包包:

19、familieslgamlss程序包中的gamlss函数：( 46 )补充：补充：gamlss函数函数l一般线性回归M1=gamlss(yx,data=dt)l存在异方差M2=gamlss(yx,sigma.formula=x,data=dt)l广义线性回归M3=gamlss(yx,family=GA,mu.link=log,data=dt)lGamlss回归M4=gamlss(yx,family=GA,sigma.formula=x,mu.link=log,sigma.link=log) ( 47 )产生模拟数据产生模拟数据( 48 )建立伽马、逆高斯回归模型建立伽马、逆高斯回归模型分别用g

20、lm、gamlss函数对模拟数据集建立伽马回归模型和逆高斯回归模型，全都使用对数连接函数( 49 )模型比较：模型比较：AIC、SBCSBCAICmGA2700.1618 688.6896mIG2 727.725 716.2529( 50 )模型比较：残差图模型比较：残差图Gamma Inverse Gaussian( 51 )模型比较：分位残差图模型比较：分位残差图Gamma Inverse Gaussian( 52 )模型比较：蠕虫图模型比较：蠕虫图( 53 )模型比较：拟合值模型比较：拟合值( 54 )3.5模型推广模型推广l gamlss包补充l删失数据的回归模型l有限混合回归模型Wh

21、at is gamlss lGAMLSS模型: 可以把每个分布参数都表示成解释变量的线性函数与随机效应的线性函数之和，因变量不再局限于 传统模型中的指数分布族 ，只要求密度函数对分布参数 一阶和二阶可导 。它涵盖了传统的线性模型 、 广义线性模型 、广义可加模型和广义线性混合模型等 ， 还可以衍生出半参数可加模型 、非线性半参数模型 、非线性参数模型等 一系列 形式 。( 55 )Gamlss packagelgamlss(拟合gamlss模型)lgamboostLSS(boosting algorithms)lgamlss.cens(拟合删失数据censored)lgamlss.dist（构

22、造分布）lgamlss.mx（拟合有限混合分布）lgamlss.nl（拟合参数化的非线性模型）lgamlss.tr（拟合截断数据truncated） ( 56 )gamlss函数函数lgamlss(formula = formula(data), sigma.formula = 1, nu.formula = 1, tau.formula = 1, family = NO(), data = sys.parent(), weights = NULL, contrasts = NULL, method = RS(), start.from = NULL, mu.start = NULL, sig

23、ma.start = NULL, nu.start = NULL, tau.start = NULL, mu.fix = FALSE, sigma.fix = FALSE, nu.fix = FALSE, tau.fix = FALSE, control = gamlss.control(.), i.control = glim.control(.), .)( 57 )ldata(rent) lmod=c,c,y0)status=ifelse(y=c,0,1)把y表示成为右删失数据的格式y=survival:Surv(time=y, event=status)( 63 )拟合步骤拟合步骤lli

24、brary(gamlss.cens) 生成右删失的伽马分布，默认名称为GArc gen.cens(family=GA,type=right) 建立右删失的伽马回归模型 mGA-gamlss(yx1+x2+x3+x4, data = dt, family = GArc, mu.link=log) plot(mGA) wp(mGA)( 64 )( 65 )( 66 )term.plot(mGA)Error in eval(predvars, data, env) : invalid envir argument of type closure( 67 )截断截断分布分布（truncated dis

25、tribution） library(gamlss.tr) 生成0截断t分布（向左截断）(left,right,both) gen.trun(0,family=“TF”)默认名称：分布+tr dTFtr pTFtr qTFtr rTFtr TFtr 生成0截断观测值 sam-rTFtr(1000,mu=10,sigma=5, nu=5 ) mod1-gamlss(sam1, family=TFtr) mod1 ( 68 )( 69 )有限混合回归模型有限混合回归模型l( 70 ) 1( )niiif xf x模拟数据模拟数据library(gamlss)install.packages(gam

26、lss.mx)library(gamlss.mx)set.seed(111)n=5000#第一个伽马分布的参数mu1=100sigma1=0.2#第二个伽马分布的参数mu2=300sigma2=0.3#模拟混合伽马分布的观察值，第一个伽马分布的占比为40%index=runif(n)x=ifelse(index=0.4,rGA(n,mu=mu1,sigma=sigma1),rGA(n,mu=mu2,sigma=sigma2)#混合伽马分布模拟值得直方图hist(x,breaks=100,main=)( 71 )( 72 )估计混合伽马分布的参数lmod=gamlssMX(x1,family =

27、 GA,k=2)输出参数估计值lmod( 73 )参数mu1sigma1mu2sigma212实际值1000.23000.30.40.6估计值100.580.2049300.0650.30240.4013 0.5987模拟数据模拟数据set.seed(111)n=5000#第一个伽马分布的参数mu1=100sigma1=0.2#第二个伽马分布的参数mu2=300sigma2=0.3#第三个伽马分布的参数mu3=800sigma3=0.1#模拟混合伽马分布的观察值，第一个伽马分布的占比为40%index=runif(n)x=ifelse(index=0.8,rGA(n,mu=mu2,sigma=

28、sigma2),rGA(n,mu=mu3,sigma = sigma3)#混合伽马分布模拟值得直方图hist(x,breaks=100,main=)( 74 )多多峰数据峰数据( 75 )lmod=gamlssMX(x1,family = GA,K=3) mod( 76 )参数mu1sigma1 mu2sigma2 mu3sigma3 123实际值1000.23000.38000.10.40.20.4估计值100.4810.204296.19 0.299800.311 0.10.401 0.190 0.409模拟数据模拟数据set.seed(111)n=5000#模拟次数x1=rgamma(n

29、,2,1)#解释变量x1，x2=rbinom(n,1,0.4)#分类解释变量x2#参数的真实值b0=2，b1=-0.15，b2=-0.2，c0=7，c1=0.2，c2=0.25#两个伽马分布的均值mu1=exp(b0+b1*x1+b2*x2)，mu2=exp(c0+c1*x1+c2*x2)，mu3=exp(d0+d1*x1+d2*x2+d3*x3)，dt=data.frame(x1,x2,x3,mu1,mu2,mu3,index)#模拟因变量的观察值index=runif(n)dt$y=ifelse(index=0.4,rGA(n,mu1), rGA(n,mu)mod=gamlssMX(yx1

30、+x2,data = dt,family = GA,K=2)mod( 77 )应用过程应用过程library（gamlss）library(gamlss.mx)mod=gamlssMX(yx1+x2,data=dt,family=“GA”,K=2)( 78 )参数实际值估计值b021.9684b1-0.15-0.151b2-0.2-0.2272c076.9816c10.20.2144c20.250.2515多多个伽马的混合分布个伽马的混合分布set.seed(111)n=5000#模拟次数x1=rgamma(n,2,1)#解释变量x1；x2=rbinom(n,1,0.4)#分类解释变量x2；x

31、3=rnorm(n,1,0.8)#分类解释变量x3#参数的真实值b0=2；b1=-0.15；b2=-0.2；b3=0.1；c0=7；c1=0.2；c2=0.25；c3=0.3；d0=3；d1=0.1；d2=-0.1；d3=0.8#两个伽马分布的均值mu1=exp(b0+b1*x1+b2*x2)mu2=exp(c0+c1*x1+c2*x2)mu3=exp(d0+d1*x1+d2*x2+d3*x3)index=runif(n)dt=data.frame(x1,x2,x3,mu1,mu2,mu3,index)#模拟因变量的观察值dt$y=ifelse(index=0.8,rGA(n,mu2),rGA

32、(n,mu3)hist(dt$y,breaks=100)mod=gamlssMX(yx1+x2+x3,data = dt,family = GA,K=3)mod( 79 )参数实际值估计值b022.094b1-0.15-0.165b2-0.2-0.179b30.1-0.0831c076.967c10.20.258c20.250.249c30.3-0.0123d032.899d10.10.112d2-0.1-0.076d30.80.82910.40.41820.20.19530.40.387( 80 )从从因变量的直方分布图上并看不出是多个伽马混合分布因变量的直方分布图上并看不出是多个伽马混合分布( 81 )拟合数据时拟合数据时AIC比较比较lfam=c(GA,IG,LOGNO,BCT,NO) m1=list() for (i in 1:5) m1fami=GAIC(gamlssMX(x1,K=2,family=fami,n.cyc=30,trace=F) sort(unlist(m1) AIC比较 GA BCT LOGNO NO IG 59218 59226 59235 59330 61079 ( 82 )stepGAIC,addterm,droptermlstepGAIC(object, sco