《风险模型》课件 第三部分 GLM

GLM：引言孟生旺线性回归模型（LM）广义线性模型（GLM）线性混合模型（LMM）广义线性混合模型（GLMM）广义可加模型（GAM）位置、尺度和形状参数的GAMLSS回顾：线性回归模型线性回归模型（LinearRegressionModel,LRM）是统计建模的基石。几乎所有复杂模型——广义线性模型（GLM）、广义加性模型（GAM）、生存模型、混合模型、甚至神经网络——都在某种意义上是线性模型的推广。理解线性回归=理解GLM的“原型”。

假设含义线性性（Linearity）对参数的线性，而非对自变量的线性独立性（Independence）同方差性（Homoscedasticity）正态性（Normality）线性回归模型的基本假设

拟合优度

诊断内容方法典型表现线性关系残差图vs.拟合值曲线形态说明非线性异方差性Breusch–Pagan检验残差方差随拟合值变化正态性Q-Q图/Shapiro-Wilk检验重尾或偏斜异常值Cook’sDistance局部高影响观测点模型诊断问题原因GLM解决方案响应变量不服从正态分布（如二项、泊松）正态假设被破坏指数离散族预测值可能超出取值范围（如概率<0或>1）恒等链接函数不适用使用对数/逻辑连接函数方差不恒定均值-方差关系变化方差函数随分布变化线性回归的局限性LM在精算应用中的局限性正态假设不能满足保险精算领域的实际需要。例：保险损失可能是计数型变量、大于零的右偏型变量。因变量的方差常常随均值而变化。例：损失金额数据。自变量之间可能通过乘法关系对因变量产生影响。例：保险费率通常表示为各种费率因子的乘积

线性回归广义线性模型响应变量连续型，正态指数离散族链接函数恒等可选（logit,log,probit等）方差结构常数与均值相关拟合方法最小二乘极大似然应用领域普通连续数据分类数据、计数数据、比例数据等小结：从线性回归到广义线性回归GLM的三个构成要素随机成分：指数离散族（ExponentialDispersionFamily,EDF）。指数离散族的方差随其均值的变化而变化。系统成分：连接函数：

g单调可导。

分布假设（EDF）：连接函数：注：g表示任意连接h表示自然连接指数离散族（EDF）孟生旺单参数指数族

（ExponentialFamily，EF）密度函数：矩母函数：均值和方差：指数离散族（ExponentialDispersionFamily，EDF）X

~单参数指数族~指数离散族

EDF中的常用分布：正态分布泊松分布二项分布伽马分布逆高斯分布例：正态分布属于EDF

例：泊松分布属于EDF

例：伽马分布属于EDF

例：逆高斯分布属于EDF分布应用场景正态分布连续型响应二项分布比例/分类泊松分布计数型响应伽马分布正偏连续型响应逆高斯分布时间/持续期建模EDF小结EDF的均值和方差孟生旺EDF的均值和方差

证明（1）：证明（2）：两边关于

求导数：注：Fisher信息与信息矩阵

称作方差函数，记为例：泊松分布的方差函数为

故

注：积分中的常数项与μ无关，可以吸收进a(y,

)伽马分布的密度函数方差函数决定EDF的具体形式练习：伽马分布的方差函数为

2，即方差可以表示为

2逆高斯分布的方差函数为

3，即方差可以表示为

3注：在伽马分布中，

在逆高斯分布中，分布方差函数正态二项泊松伽马逆高斯TweedieEDF的矩母函数孟生旺矩母函数：EDF的矩母函数密度函数：

证明（可略）：EDF的累积量母函数（CumulantGeneratingFunction）

例：计算EDF的均值和方差。例（EDF的卷积）：若相互独立，则其中证：例：Tweedie分布=泊松-伽马复合分布EDF的矩母函数：Y~Tweedie分布：的矩母函数：

的矩母函数（Y~Tweedie）：S的矩母函数：EDF的KL散度和Bregman散度孟生旺在评价GLM时，均方误差（MSE）不再适用，因为GLM的方差结构依赖于均值。取而代之的，是一种由似然原理（LikelihoodPrinciple）

导出的度量——偏差（Deviance）偏差是以KL散度为基础的似然距离，而KL散度又是Bregman散度的一种特例。本节主要内容：KL散度Bregman散度KL散度（Kullback–LeiblerDivergence）

KL散度的统计学意义：EDF的KL散度：EDF的KL散度应用正则连接，令h(

)

=将

0看做饱和模型，

0=h(y);将

1看做当前模型，

1=h(

)，则注：上式乘以2，即为EDF的单位偏差（下节内容）。Bregman散度

Bregman散度的统计意义：统一框架：为各种损失函数和偏差度量提供统一视角最优性保证：确保在正确模型设定下估计的最优性模型诊断：提供基于散度的模型诊断工具不同的统计模型本质上是在不同的几何结构下寻找数据的最佳近似。所有EDF的负对数似然都对应一个Bregman散度！若令：则：EDF的Bregman散度：注：上式乘以2，即为EDF的单位偏差（下节内容）。

Bregman散度与KL散度的关系

分布KL散度Bregman散度二项分布

(伯努利)泊松分布单位偏差

（UnitDeviance）孟生旺单位偏差（UnitDeviance）上式在μ=y处达到最大（饱和模型）：单位偏差：当前模型与饱和模型之间的距离

单位偏差的三种解释（正比于）：饱和模型与当前模型的对数似然之差饱和模型与当前模型的KL散度饱和模型与当前模型的Bregman散度单位偏差可表示为对数似然之差：应用单位偏差，EDF的密度函数可以表示为：最大化似然=最小化单位偏差饱和模型与当前模型的对数似然之差单位偏差可表示为KL散度：应用自然连接，令h(

)

=

将

0看做饱和模型，

0=h(y);将

1看做当前模型，

1=h(

)，则注：上式乘以2，即为单位偏差。单位偏差可表示为Bregman散度：Bregman散度：注：是严格凸函数令：则：注：上式乘以2，即为单位偏差。例：求泊松分布的单位偏差例：求正态分布的单位偏差单位偏差的比较：均值=2练习：求伽马分布的单位偏差。应用均值参数和离散参数1/练习：求逆高斯分布的单位偏差,应用均值参数和离散参数1/

EDF的鞍点近似孟生旺EDF的鞍点近似将方差函数修正为，可以提高精度，且适用于y=0的情形。（鞍点近似）饱和模型的密度

例：鞍点近似对正态分布是精确的。在正态分布中，V(μ)=1,故V(y)=1练习：鞍点近似对于逆高斯是精确的。定理：如果鞍点近似成立，尺度化单位偏差近似服从自由度为1的卡方分布，即需要证明：尺度化单位偏差的矩母函数为只需证明：单位偏差的矩母函数为单位偏差的矩母函数（如何消除积分？）例：鞍点近似对正态分布是精确的，故其尺度化单位偏差服从自由度为1的卡方分布。讨论题：鞍点近似对逆高斯分布是精确的，故尺度化单位偏差服从自由度为1的卡方分布。鞍点近似成立的经验准则（Smyth&Verbyle，1999）：

例：对于泊松分布，

=1，V(y)=y，例：对于伽马分布，V(y)=y2形状参数

=1/

3例：假设泊松分布的参数为1.2，伽马分布的参数为（μ=100，

=0.5），分别用鞍点近似计算它们的概率，并绘图比较鞍点近似的精度。

lam=

c(0.0001,0.001,0.01,seq(0.1,10,by=0.1))

fED=

function(mu){

y=

seq(0,100,1)

#对于给定的lambda，计算单位偏差的均值

sum(dpois(y,lambda=mu)*

poisson()$dev.resids(y,mu,wt=1))

}

ED=

Vectorize(fED)(lam)#对于不同的lambda，计算单位偏差的均值

plot(ED~

lam,type="n",main=

'泊松分布的单位偏差',xlab=

expression(lambda),ylab=

expression(E(d(y,mu))))

polygon(x=c(-1,-1,12,12),y=c(0.95,1.05,1.05,0.95),col="gray",border=NA)

lines(ED~

lam,col=2,type="l",lty=2,lwd=2)

abline(h=1)phi=

seq(0.01,0.5,by=0.01)

fED=

function(phi){

y=

seq(0,2000,1)

mu=

100

shape=

1/phi;scale=

phi*mu

#对于给定的phi，计算单位偏差的均值

sum(dgamma(y,shape=shape,scale=scale)*

Gamma()$dev.resids(y,mu=mu,wt=1))

}

ED=

Vectorize(fED)(phi)#对于不同的phi，计算单位偏差的均值

plot(ED~

phi,type="n",main=

'伽马分布的单位偏差',xlab=

expression(phi),ylab=

expression(E(d(y,mu))))

polygon(x=c(0,0,0.5,0.5),y=c(0,0,1.05*0.5,0.95*0.5),col="gray",border=NA)

lines(ED~

phi,col=2,type="l",lty=2,lwd=2)

lines(phi,phi,lty=1,lwd=2)R中计算单位偏差：poisson()$dev.resids(y,mu,wt=1)Gamma()$dev.resids(y,mu=mu,wt=1)从单位偏差到偏差（

μ

已知）偏差：尺度化偏差：（假设鞍点近似成立）从偏差到剩余偏差（

μ

未知）（residualdeviance）剩余偏差：尺度化剩余偏差：k+1是模型的参数个数注：饱和模型的剩余偏差等于0（假设鞍点近似成立）连接函数孟生旺连接函数:名称连接函数等值对数平方根幂函数logitprobitlog-logComplementarylog-log

泊松分布：伽马分布：分布自然连接函数正态分布二项分布泊松分布伽马分布逆高斯分布常见分布的自然连接函数

自然连接函数的作用：简化回归系数的估计方程（参见回归系数估计部分）EDF的参数估计孟生旺EDF的参数估计自然参数的极大似然估计：均值参数的极大似然估计：自然参数与均值参数的关系：MLE的平衡性：注：确保总体预测水平的合理性。是一致最小方差无偏估计（UMVU）观察值：的无偏估计：Cramer-Rao下界:Fisher信息：的方差达到Cramer-Rao下界：得分：Fisher信息：下界密度函数：模型性能评价孟生旺模型性能评价泛化损失如何计算泛化损失？交叉验证Bootstrap泛化损失（GeneralizationLoss）历史观察数据：对新数据Y

的预测值：如何评价模型的预测性能？计算泛化损失（generalizationloss）历史观察数据：对新数据Y

的预测：最常见的期望泛化损失：均方误差（meansquarederror）。BiasEstimationvarianceProcess

variance定理：均方误差可以表示为（证明见下页）期望泛化损失与泛化损失均方误差的分解：注：对Y的最优预测是E(Y)，此时，偏差和方差都等于零。但实际上E(Y)未知,需用A(Yn)预测。在给定Yn的条件下，条件均方误差（conditionalMSE）称作泛化损失：注：均方误差（MSE）是期望泛化损失例：指数离散族的MSE假设：自然参数的极大似然估计值为：均值参数的极大似然估计值：即，则MSE为：（证明见下页）用均值参数的极大似然估计值预测Y：注：当时，估计方差趋于零。估计方差过程方差MSE源于正态，不适用于EDF：

例：如果索赔次数的真实值为1

，预测值为1.2，则MSE=0.04，如果索赔次数的真实值为0，预测值为0.2，则MSE=0.04，但从业务逻辑看，第二种情况下的误差更大。

MSE不适用于EDF，如何解决？基于偏差定义期望泛化损失和泛化损失例：若用预测Y，则期望偏差为其中：称作估计风险定理：若用

预测Y，则期望偏差为类似于过程方差不能消除估计风险无法分解为估计的偏差和方差在给定Yn的条件下，偏差为：例：正态分布的估计风险（estimationrisk）biasEstimationvariance注：正态分布下最小化估计风险=最小化均方误差例：泊松分布的估计风险（estimationrisk）假设

A(Yn)是

的无偏估计，则若在无偏估计中增加偏差c，则（*）变为泊松分布的估计风险:（*）注：无偏变有偏，估计风险变小！建议：应用期望偏差对模型进行评价，应基于无偏预测。常用分布的单位偏差对于Tweedie单位偏差的统一表达式：定理：求解均值函数，Bregman散度是唯一的一致性损失函数，故最小化Bregman散度，即可求得均值函数。推论：在EDF中，单位偏差是一致性损失函数，故最小化期望偏差即可求得均值函数。求解上式，需要已知Y的分布，通常是未知的，故考虑经验估计：交叉验证（CrossValidation）问题：当自然参数

未知时，如何估计下面的期望偏差泛化损失？经验估计：以均值参数的MLE为例对于给定的，如何估计偏差泛化损失？（上式积分号内部）偏差泛化损失的经验估计：样本外偏差样本内偏差现实问题：样本量不大，不宜划分为训练集和测试集？解决之道：留一交叉验证（Leave-one-outCross-Validation）K

折交叉验证，通常K=10分层K

折交叉验证留一交叉验证（Leave-one-outCross-Validation）K折交叉验证，随机划分K组，通常K=10分层K

折交叉验证：前K

个随机分到K组，随后K

个再随机分到K

组,……练习：在MTPL数据集中，ClaimNb是索赔次数（Yi），Exposue是风险暴露（wi），索赔频率定义为假设索赔次数服从泊松分布，通过最小化泊松偏差估计泊松参数。基于所有样本数据，计算样本内泊松偏差计算十折交叉验证泊松偏差基于十折交叉验证，计算泊松偏差的标准误library(CASdatasets)data(freMTPL2freq)dat<-freMTPL2freq[,-2]dat$VehGas<-factor(dat$VehGas)data(freMTPL2sev)sev<-freMTPL2sevsev$ClaimNb<-1dat0<-aggregate(sev,by=list(IDpol=sev$IDpol),FUN=sum)[c(1,3:4)]names(dat0)[2]<-"ClaimTotal"dat<-merge(x=dat,y=dat0,by="IDpol",all.x=TRUE)dat[is.na(dat)]<-0dat<-dat[which(dat$ClaimNb<=5),]dat$Exposure<-pmin(dat$Exposure,1)sev<-sev[which(sev$IDpol%in%dat$IDpol),c(1,2)]dat$VehBrand<-factor(dat$VehBrand,levels=c("B1","B2","B3","B4","B5","B6","B10","B11","B12","B13","B14"))

library(data.table)fwrite(dat,"D:\\dat.csv")基于相同的样本观察值，如何评价模型？AIC可以证明：在一定条件下，最小化AIC

=最小化留一交叉验证均方误差。Bootstrap非参数Boostrap参数Bootstrap非参数Boostrap参数BoostrapGLM回归系数估计孟生旺GLM回归系数的估计

Y的分布：对数似然：假设有n个相互独立的观测值，则联合对数似然：求各个参数的一阶偏导，令其等于0，即得参数的估计方程。【定理】GLM回归系数的极大似然估计是下列估计方程的解：其中：注：离散参数

不影响回归系数的估计值！下面证明：【证明】回归系数的极大似然估计是下述方程组的解：应用链式法则：上式表明，自然连接满足平衡性：每个类别的预测值之和=观察值之和特别地，如果使用自然连接，g(.)=h(.)将其带入可得：求解估计方程的方法之一：Newton迭代法估计方程：泰勒展开：变形：迭代公式：一阶导数向量称作得分向量，第j个元素：二阶导数矩阵称作海森矩阵（Hessianmatrix），(

j,

h)位置的元素：Newton迭代过程（不易使用）：（1）设置初始值

（2）根据初始值计算和（3）应用迭代公式：

（4）计算，判断其是否小于阈值（如）（5）如果大于阈值，则将作为新的初始值带入迭代公式。重复步骤（3）~（5），直到两个估计值之差的绝对值小于阈值。

求解估计方程的方法之二：迭代加权最小二乘法

其中：得分向量和信息矩阵：

z称为工作因变量（workingresponse），R

称为工作残差（working

residual）【定理】估计GLM回归系数的迭代加权最小二乘法：其中Newton迭代法：证明：将海森矩阵用期望海森矩阵（负的信息矩阵）代替，即：回归系数估计的标准误参数估计值的方差：协方差矩阵对角线元素参数估计值的标准误：参数估计值的协方差矩阵：例：Logistic回归的迭代加权最小二乘法流程

极大似然=最小偏差极大似然：极大似然=最小偏差：例：yi133567910xi-1-1000011建立如下GLM，估计模型的回归系数：yi133567910xi-1-1000011dt=

data.frame(y=

c(1,3,3,5,6,7,9,10),

x=

c(-1,-1,0,0,0,0,1,1))

b0=

c(2,1)

X=

model.matrix(~

x,data=dt)#设计矩阵，大写X

m=

0

repeat{

m=

m+

1

b=

b0

W=

diag(c(1/(X%*%b)))

z=

dt$y

b0=

solve((t(X)%*%W%*%X))%*%t(X)%*%W%*%z

if(max(abs(b-

b0))<

1e-8)break

}

b#迭代加权最小二乘法的参数估计值##[,1]

##(Intercept)5.500000

##x3.586957m#迭代次数##[1]7glm(y~

x,data=dt,family=

poisson(link=

'identity'))$coef##(Intercept)x

##5.5000003.586958讨论题：yi133567910xi-1-1000011wi0.40.80.91.21.10.70.90.8建立如下GLM，估计模型的回归系数：

#数据输入

y<-

c(1,3,3,5,6,7,9,10)

x<-

c(-1,-1,0,0,0,0,1,1)

w<-

c(0.4,0.8,0.9,1.2,1.1,0.7,0.9,0.8)

#建立数据框

data<-

data.frame(y=y,x=x,w=w)

#方法1：使用glm函数，设置offset参数

model<-

glm(y~x,family=

poisson(link=

"log"),

offset=

log(w),data=data)

#显示结果

summary(model)

#提取系数

beta_hat<-

coef(model)

cat("\n===回归系数估计结果===\n")

cat("β₀=",round(beta_hat[1],4),"\n")

cat("β₁=",round(beta_hat[2],4),"\n")

#模型诊断和预测

cat("\n===模型拟合效果===\n")

fitted_values<-

fitted(model)

residuals<-

residuals(model,type=

"pearson")

results<-

data.frame(

Observed=y,

Weight=w,

x=x,

Fitted=

round(fitted_values,3),

Residuals=

round(residuals,3)

)

print(results)

#计算模型的偏差

cat("\n模型偏差:",round(deviance(model),4),"\n")

cat("AIC:",round(AIC(model),4),"\n")#方法2：手动验证-使用迭代加权最小二乘法

#初始化参数

beta<-

c(0,0)#初始值

max_iter<-

100

tolerance<-

1e-8

for(iterin

1:max_iter){

#计算当前预测值

eta<-beta[1]+beta[2]*x#线性预测器

mu<-w*

exp(eta)#均值

#计算权重矩阵W(对角矩阵)

W<-

diag(mu)

#计算工作响应变量z

z<-eta+(y-mu)/mu

#构建带偏移量的设计矩阵

X<-

cbind(1,x)

#迭代加权最小二乘更新

beta_new<-

solve(t(X)%*%W%*%X)%*%

t(X)%*%W%*%z

#检查收敛

if(max(abs(beta_new-beta))<tolerance){

break

}

beta<-beta_new

}

cat("\n===手动IWLS验证结果===\n")

cat("β₀=",round(beta[1],4),"\n")

cat("β₁=",round(beta[2],4),"\n")

cat("收敛于第",iter,"次迭代\n")

应用glm函数求得的估计值：β₀=1.7354β₁=0.4815编写IWLS的程序代码求得的结果：β₀=1.7354β₁=0.4815收敛于第6次迭代GLM离散参数的估计孟生旺离散参数的估计离散参数的取值不会影响回归系数的极大似然估计值。理论上，可以应用极大似然法估计离散参数，即求解方程：但应用极大似然法，没有统一形式的估计方程。实际应用中，通常使用矩估计法：有两种方法若因变量服从EDF，其方差可以表示为：离散参数的一个无偏估计为：（证明参见下页）离散参数的矩估计近似无偏：参数个数=k+1此式的期望值=

注：估计离散参数时，最好使用未经汇总的数据。

数据汇总后，会减少观测值个数n，影响离散参数估计的稳定性。数据汇总不会对回归系数的估计值产生影响。索赔次数性别1男2女2男0女未汇总平均索赔次数性别1.5男1女汇总GLM的尺度化剩余偏差：

是饱和模型的对数似然

是当前模型的对数似然注：饱和模型的偏差=0估计离散参数的第二种方法：基于偏差鞍点近似成立的情况下，尺度化剩余偏差近似服从n–k–1的卡方分布：基于剩余偏差估计离散参数（无偏估计：下式的期望值=

）GLM的检验和比较孟生旺GLM的检验和比较嵌套模型：卡方检验和F检验非嵌套模型：Vuong检验嵌套模型的检验

离散参数

未知：F

检验嵌套模型（例）：尺度化偏差相减，近似服从卡方分布（

r参数个数之差）嵌套模型的卡方检验（

已知，如

=1）嵌套模型的F检验（

未知）非嵌套模型：两个模型无法通过参数约束从其中一个得到另一个。例：一个泊松回归vs.一个负二项回归一个Logit模型vs.一个Probit模型一个Gamma回归vs.一个对数正态回归非嵌套模型非嵌套模型的比较：Vuong检验是一种用于比较两个非嵌套模型的统计检验方法，由Vuong(1989)提出。如果两个模型对数据的拟合能力“没有差异”，那么它们在每个观测点上的对数似然值之差的均值应该为零。它解决了传统似然比检验只能用于嵌套模型比较的局限性。

Vuong检验的统计量

Vuong检验统计量的应用基于AIC或BIC修正的Vuong检验：

Vuong检验方法的选择：如果旨在选择最佳预测模型，使用基于AIC的Vuong修正。如果旨在找到“真实模型”（当它存在于候选模型中时），使用基于BIC的Vuong修正。如果只关心纯粹的样本内拟合优度，而不考虑复杂度，可以使用标准Vuong检验，但这在现实中较少使用。Vuong检验示例：#加载必要包

library(MASS)

library(pscl)

#模拟数据

set.seed(123)

n<-

200

age<-

rnorm(n,50,10)

bmi<-

rnorm(n,25,3)

smoking<-

rbinom(n,1,0.3)

#生成索赔次数（存在过度离散）

mu<-

exp(0.1*age+

0.2*bmi+

0.5*smoking)

claims<-

rnbinom(n,size=1,mu=mu)#负二项分布

insurance_data<-

data.frame(claims,age,bmi,smoking)

#拟合两个竞争模型

model_pois<-

glm(claims~age+bmi+smoking,family=poisson,data=insurance_data)

model_nb<-

glm.nb(claims~age+bmi+smoking,data=insurance_data)

#执行Vuong检验

vuong_result<-

vuongtest(model_pois,model_nb)

print(vuong_result)检验结果：VuongNon-NestedHypothesisTest-Statistic:

(test-statisticisasymptoticallydistributedN(0,1)underthe

nullthatthemodelsareindistinguishible)

-------------------------------------------------------------

Vuongz-statisticH_Ap-value

Raw-4.572model2>model12.45e-06

AIC-corrected-4.572model2>model12.45e-06

BIC-corrected-4.572model2>model12.45e-06解释:检验统计量V=-4.572p值<0.001，强烈拒绝两个模型无差异的原假设由于V为负且显著，结论是模型2（负二项回归）优于模型1（泊松回归）不拒绝H₀并不意味着两个模型一样好，也可能意味着它们一样差，或者样本量太小，无法检测出它们之间的差异。并非“模型正确性”检验：Vuong检验只告诉我们哪个模型相对更好，有可能两个模型都是错误的。Vuong检验依赖于模型似然函数的正确设定。如果模型本身的似然函数计算有误，检验结果也将无效。在R中，vuong()

函数（通常在

pscl

包中）可以方便地执行此检验。应用Vuong检验注意事项：非嵌套模型的比较：直接使用AIC或BIC

应用AIC的经验规则：（1）AIC之差小于2.5时，表明两个模型没有明显差异。（2）AIC之差大于10时，AIC较小的模型明显较优（3）样本量大于256且AIC之差大于2.5小于6，则AIC较小的模型较优。（4）样本量大于64且AIC之差大于6小于9，则AIC较小的模型较优应用BIC的经验规则：（1）BIC之差在0~2之间，表明两个模型存在微弱差异。（2）BIC之差在2~6之间，表明两个模型存在一定差异。（3）BIC之差在6~10之间，表明两个模型存在显著差异（4）BIC之差大于10，表明两个模型存在非常显著的差异。GLM残差与模型诊断孟生旺GLM的基本假设数据中没有异常值.使用了正确的连接函数.线性预测项中包含了所有的重要解释变量，每个解释变量使用了正确的尺度.使用了正确的方差函数.离散参数是常数.因变量的观察值相互独立.因变量来自特定的指数离散族.残差残差01原始残差02皮尔逊残差03偏差残差04分位残差

y1异常，y2正常

原始残差（因变量残差）

局限性：在分布偏斜严重时（如二项分布中概率接近0或1，泊松分布中均值很小），即使模型正确，皮尔逊残差的分布也可能不对称，不太像正态分布。Pearson残差Deviance残差

偏差衡量当前模型与饱和模型之间的差距。偏差残差是考虑了似然函数本身，更符合GLM的似然框架。与皮尔逊残差相比，偏差残差通常更接近正态分布，尤其是在小到中等规模的样本中。在模型诊断中，应优先使用偏差残差。近似服从标准正态分布。当需要严格检验分布假设（特别是通过Q-Q图）时，分位残差是最强大、最可靠的工具。

par(mfrow=c(1,2))x=seq(0,400,1)y=pgamma(x,shape=2,scale=50)plot(x,y,type='l',col=2,xlab='y',ylab='pgamma(y)',lwd=2)x1=seq(-3,3,0.1)y1=pnorm(x1)plot(x1,y1,type='l',col=3,xlab='r=qnorm(pgamma(y))',ylab='pgamma(y)',lwd=2)分位残差

par(mfrow=c(1,2))x=seq(0,400,1)y=pgamma(x,shape=2,scale=50)plot(x,y,type='l',col=2,xlab='y',ylab='pgamma(y)',lwd=2)x1=seq(-3,3,0.1)y1=pnorm(x1)plot(x1,y1,type='l',col=3,xlab='r=qnorm(pgamma(y))',ylab='pgamma(y)',lwd=2)Gamma(shape=2,scale=50)Norm(0,1)

dpois(lambda=1)Norm(0,1)标准化残差将杠杆效应考虑在内后，对普通残差（如偏差残差）的进行修正，能够：公平地比较不同杠杆位置上的观测点的拟合优劣。精确地识别出真正的统计异常值。作为计算Cook距离的基础，从而系统地发现那些对模型参数估计有巨大影响力的数据点。

基于残差的模型诊断流程dt=data.frame(y=c(1,3,3,5,6,7,9,10),x=c(-1,-1,0,0,0,0,1,1))mod=glm(y~x,data=dt,family=poisson(link=log))#残差resid(mod,type='pearson')#Pearson残差resid(mod,type='deviance')#Deviance残差rstandard(mod)#标准化Deviance残差mod$residuals#相对残差,等价于下式(dt$y-fitted(mod))/fitted(mod)#分位残差library(statmod)qresid(mod)R中的残差dt=data.frame(y=c(1,3,3,5,6,7,9,10),x=c(-1,-1,0,0,0,0,1,1))

mod=glm(y~x,data=dt,family=poisson(link=

'log'))

#残差

resid(mod,type=

'pearson')#Pearson残差##12345678

##-0.89880.3952-0.84400.06310.51670.9703-0.28280.0352resid(mod,type=

'deviance')#Deviance残差##12345678

##-1.01810.3799-0.90890.06280.49830.9097-0.28720.0352rstandard(mod)#标准化Deviance残差

##12345678

##-1.19590.4462-0.97980.06770.53720.9807-0.38570.0472mod$residuals#相对残差,等价于下式##12345678

##-0.58150.2556-0.38280.02860.23440.4401-0.08990.0112(dt$y-fitted(mod))/fitted(mod)##12345678

##-0.58150.2556-0.38280.02860.23440.4401-0.08990.0112#分位残差

library(statmod)qresid(mod)##[1]-1.15860.3137-0.8726-0.07770.59471.0112-0.09950.0469练习根据索赔次数数据，在负二项分布假设下，建立仅含截距项的GLM，估计模型参数，绘制分位残差的QQ图。n=0:5#索赔次数policy=c(1235,521,98,32,4,1)#保单数#负二项的极大似然估计与随机分位残差

n=

0:5

policy=

c(1235,521,98,32,4,1)

library(fitdistrplus)data=

rep(n,policy)

fit=

fitdist(data,‘nbinom’,

method=

‘mle’)

r=

fit$estimate[1]

mu=

fit$estimate[2]

set.seed(2203)

Fn=

c(0,pnbinom(n,size=r,mu=mu))#负二项分布函数

u=

NULL

for(i

in

1:6){

u=

c(u,runif(policy[i],Fn[i],Fn[i+1]))

}

Qn=

qnorm(u,0,1)#随机分位残差

qqnorm(Qn);qqline(Qn)方法2：应用statmod程序包n=

0:5

#索赔次数

policy=

c(1235,521,98,32,4,1)#保单数

num=

rep(n,policy)

dt=

data.frame(num=num)

library(statmod)

library(MASS)

mNB=

glm.nb(num~

1,data=dt)

qqnorm(qresid(mNB));qqline(qresid(mNB))n=

0:5

#索赔次数

policy=

c(1235,521,98,32,4,1)#保单数

num=

rep(n,policy)

dat=

data.frame(num=num)

library(gamlss)mod=

gamlss(num~

1,family=NBI,data=

dat)#QQ图

plot(mod)方法3：应用gamlss程序包GLM拟似然方法孟生旺GLM需要完全指定概率分布，但在实际应用中，常常遇到两个标准GLM无法妥善处理的问题：过度离散：如在拟合泊松回归时，假设方差等于均值。但现实数据中，方差往往大于均值。忽略过度离散会导致标准误被低估，导致错误的推断。分布形式未知：我们可能对均值和方差之间的关系有一个合理的猜测（例如，方差与均值的平方成正比），但无法确定具体的概率分布形式。拟似然（Quasi-Likelihood）方法应运而生，它提供了一种强大而灵活的“分布无关”的替代方案。背景拟似然的核心思想核心思想：不必完全指定整个概率分布，而只需指定其前两阶矩（均值和方差），即可进行可靠的统计推断。放弃构建一个基于真实概率分布的似然函数，而是构造一个“拟似然函数”。这个构造的函数具备标准似然函数的某些关键性质，从而使得基于它的推断（如参数估计、假设检验）是渐近有效的。拟似然的构造

基于拟似然的参数估计

应用场景局限性无法使用基于似然的工具进行模型比较，如：不能使用似然比检验：拟似然不是真实似然，其差值不再服从卡方分布。不能使用AIC/BIC：这些信息准则依赖于真实似然值。不能进行概率性预测：无法计算预测区间，因为缺乏完整的分布。如果方差函数设定错误，估计结果可能有偏或不一致。其优良性质（如估计量的渐近正态性）依赖于大样本理论。当过度离散有明确的机制时（如存在过多零值、存在随机效应），使用一个完全指定的混合模型（如负二项回归、零膨胀模型、广义线性混合模型）可能在理论上更严谨。双广义线性模型

（DGLM）孟生旺双广义线性模型（DoubleGLM）模型假设：问题：如何估计这两个模型的回归系数？鞍点近似：鞍点近似的对数似然：鞍点近似的对数似然：上式中，令则鞍点近似的对数似然为：注：服从伽马分布，自然参数为，离散参数为2上述伽马中：以

为因变量，建立伽马回归（GLM2），其均值即为

。GLM2：因变量服从伽马分布，均值参数，离散参数为2例：正态分布和逆高斯分布，鞍点近似是精确的故精确服从伽马分布，均值为

，离散参数为2双GLM的参数估计：迭代求解给定，计算，估计GLM1的回归系数：均值的预测值：工作权重工作残差计算GLM1的偏差（作为GLM2的因变量）：服从伽马分布，均值参数为，离散参数为2。以为因变量构建伽马回归（GLM2）。估计GLM2的回归系数：Z

是GLM2的设计矩阵工作残差：工作权重：GLM2使用对数连接贝叶斯GLM与正则化

（版本1）孟生旺

贝叶斯GLM参数估计过程：最大后验估计与正则化

正则化的作用：当模型中解释变量过多时，模型易过拟合，导致预测方差偏大；正则化通过引入适当的偏差，显著降低预测方差，最终提升模型的整体泛化性能；注意：截距项通常不进行正则化，以免破坏GLM自然连接函数下的平衡性；若协变量的尺度差异较大，需先对协变量进行标准化或尺度变换，避免正则化强度受变量单位影响。先验分布选择与常用正则化

岭回归（RidgeRegression，L2正则化）

LASSO（L1正则化）

弹性网（ElasticNet）

弹性网的主要特性：分组选择：能实现分组选择（对线性相关的变量组，同时保留组内变量），解决LASSO的分组选择缺陷；变量选择与系数压缩：兼具LASSO的变量选择功能和岭回归的系数压缩特性，泛化性能更优；适用场景：高维数据、存在多重共线性的数据集。三类正则化方法对比方法先验分布正则化项核心功能局限性岭回归高斯分布L2范数压缩系数、缓解共线性无法变量选择LASSO拉普拉斯分布L1范数变量选择、简化模型无法分组选择、一致性不足弹性网混合分布L1+L2范数变量选择、分组选择、压缩系数正则化GLM的实现（以R语言为例）glmnet：核心包，支持岭回归、LASSO、弹性网的快速计算，适用于高维数据；glmnetUtils：辅助包，提供公式接口、数据预处理等功能，简化glmnet的使用；caret：建模框架包，支持交叉验证、模型评估等统一流程；跨语言支持：glmnet另有MATLAB和Python版本，功能一致。正则化GLM小结贝叶斯GLM中，后验分布是推断核心，复杂模型需通过MCMC方法近似；正则化的本质是通过先验分布（或惩罚项）实现偏差-方差权衡，解决过拟合问题；岭回归适合缓解共线性，LASSO适合变量选择，弹性网兼具两者优势；实践中注意：协变量标准化、截距项不惩罚、通过交叉验证选择惩罚参数；R语言的glmnet包提供了高效的正则化GLM实现。贝叶斯GLM与正则化

（版本2）孟生旺贝叶斯GLM参数估计方法

的先验分布：

与

的联合分布：

的后验分布：

的贝叶斯估计=后验均值：在给定

的条件下，计算随机变量Yn+1的期望值：给定

的条件下，Y与Yn+1相互独立问题：如何求后验分布

与后验均值

？MCMC近似计算：GibbssamplingMetropolis–Hastings(MH)algorithmsequentialMonteCarlo(SMC)samplingnon-linearparticlefiltersHamiltonMonteCarlo(HMC)algorithm参数的最大后验估计贝叶斯GLM的后验对数似然：正则化项，防止过拟合最大后验估计（MAP，MaximalaPosteriorEstimator）：先验分布的选择：注：在实际应用中，不对截距项进行正则化，以免丧失自然连接函数下GLM的平衡性。注：正则化方法中，如果协变量的尺度不同，需对协变量进行标准化或尺度变换常用正则化方法对应的先验分布：岭回归：LASSO：弹性网：先验：均值为零的高斯分布先验：Laplace分布岭回归LM的参数估计岭回归：通过正则化将部分回归系数的估计值压缩至接近零值。目标函数（需对解释变量进行标准化处理）：注：对截距项

0不惩罚。LM截距项的估计值为因变量的均值。注：LM中,如果对y

进行中心化处理，截距项为零。LM岭回归的参数估计（假设对y进行了中心化处理，截距项为零）：上式关于参数向量

求导，可得：岭回归的有效自由度=H

矩阵的迹，即H

对角线上元素之和：注：若基于数据进行选择（如交叉验证），则岭回归的有效自由度需要修正：/article/10.1007%2Fs11135-013-9949-7注：修正前后差别不是很大注：

越大，有效自由度越小。若

=0，则自由度等于模型的参数个数。岭回归的预测值：岭回归参数估计值的偏差与方差：惩罚参数

=0时，偏差为零，惩罚参数越大，偏差越大惩罚参数越大，方差越小问题：理想的平衡点在哪里？如何选择

？交叉验证/fundamentals-of-statistics/ridge-regression岭回归GLM的参数估计得分方程：负的期望Hessian矩阵：求解

的迭代公式：求解

的迭代公式：（证明见下页）帽子矩阵：预测值：LASSO：LM特例：假设Y服从高斯分布，且仅有一个协变量（忽略截距项）（困难：正则项在零点不可微）其中：得分方程：得分方程：S表示软阈值算子

=4软阈值算子推广到k个协变量的情形：迭代计算,每次求解一个beta：LASSO：GLM近端梯度下降算法（Proximalgradientdescentalgorithm）：（1）无约束梯度下降算法：（2）应用软阈值算子：注：LASSOGLM通常不满足下述的一致性：（1）当n

时，模型的参数估计值趋于其真实值（2）当n

时，可以筛选出真实的特征（变量）建议：首先用LASSO选择变量，然后建立无约束的GLM练习数据集：程序包insuranceData中的数据集dataCar。以泊松回归为例，建立索赔频率的正则化GLM。以伽马回归为例，建立案均赔款的正则化GLM。混合分布的GLM孟生旺

形状参数=(1,20,40)，尺度参数=100，每个成分密度函数的占比p

=(0.7,0.1,0.2)指数离散族的混合：示例（每个成分的密度函数是乘以概率p以后的结果）指数离散族的混合：对数似然：MLE求解困难

潜在类别视角：完全数据的对数似然（一个样本观察值的情况）：完全数据的对数似然：

EM算法在混合GLM中的应用

EM算法的迭代步骤

对于不包含协变量的空模型：

右删失数据的GLM孟生旺右删失数据：示例右删失数据不完全（右删失）数据的对数似然：完全数据的对数似然：基于EM算法的右删失数据GLM参数估计对于右删失数据的GLM，其核心思路的是：E步：基于当前的参数估计值，计算右删失数据（缺失数据）的条件期望，用该期望代替缺失的真实响应值，构造“完整数据”；M步：基于E步构造的“完整数据”，使用传统的极大似然估计方法估计GLM的自然参数；重复E步和M步，直到参数估计值收敛，得到最终的参数估计结果。右删失数据的EM算法E步：估计右删失数据的均值，用其代替右删失数据M步：定义完全数据，基于完全数据用MLE估计自然参数例：右删失Gamma左截断数据的GLM孟生旺左截断数据：示例左截断数据的对数似然函数

左截断数据个数的分布假设

对数似然函数的化简

左截断数据的EM算法估计

例：左截断Gamma

例：零截断泊松

GLM在精算中的应用孟生旺GLM在精算中的应用案均赔款索赔频率出险概率纯保费准备金评估、死亡率预测、……伽马回归、逆高斯回归案均赔款GLM孟生旺304伽马与逆高斯的比较逆高斯的优点：灵活,从对称到尖峰厚尾与伽玛分布的比较：

305Gamma：shape=3;scale=4;

IG：mu=12;phi=1/36;IG具有尖峰伽马与逆高斯的比较均值＝12,方差＝48IG具有厚尾案均赔款的含义：观测到次例：假设每次赔款服从正态分布，均值为，方差为，则案均赔款仍然服从正态分布，且有：损失金额的GLM分布广义线性模型正态分布线性回归模型伽马分布伽马回归模型逆高斯分布逆高斯回归模型假设案均赔款服从高斯分布，GLM等价于普通的线性回归：

高斯回归

迭代加权最小二乘法：

在高斯分布假设下，使用恒等连接函数，则

伽马回归伽马分布（表示为EDF形式）：均值和方差：

伽马分布(

，

)的密度函数案均赔款：假设每次赔款服从伽马分布，均值为，方差为，则案均赔款仍然服从伽马分布，且有：案均赔款的伽马回归：使用对数连接函数：

离散参数的估计：注：p为模型的参数个数案均赔款的伽马回归注：

无影响逆高斯回归与伽马分布相比，逆高斯分布具有尖峰厚尾的特征。密度函数（EDF形式）：均值和方差

案均赔款：假设每次赔款服从逆高斯分布,均值为，方差为，则案均赔款仍然服从逆高斯分布，且有：案均赔款的逆高斯回归：使用对数连接函数：

离散参数的估计：注：p为模型的参数个数案均赔款的逆高斯回归索赔频率GLM孟生旺索赔频率GLM：泊松回归负二项回归零截断泊松\负二项回归零膨胀泊松\负二项回归零调整泊松\负二项回归混合回归模型（略）glm(N~x1+x2+offset(log(t)),family=poisson(link=log))泊松回归泊松回归的迭代加权最小二乘法：使用对数连接：注：风险暴露期越长，权重越大

set.seed(111)n=500

#模拟次数x1=rgamma(n,2,1)

#解释变量x1x2=rgamma(n,2,3)

#解释变量x2x3=rbinom(n,1,0.4)

#解释变量x3x4=rbinom(n,1,0.7)

#解释变量x4#参数的真实值b0=-2;b1=0.45;b2=-0.8;b3=0.3;b4=-0.2

#线性预测项eta=b0+b1*x1+b2*x2+b3*x3+b4*x4#因变量的均值使用对数链接mu=exp(eta)#用负二项分布模拟因变量library(gamlss)y=rNBI(n,mu=mu,sigma

=

2)#把x3和x4转化为因子x3=as.factor(x3)x4=as.factor(x4)#模拟的数据集dt=data.frame(y,x1,x2,x3,x4)#输出部分模拟数据head(dt)模拟数据#初始值D0=0mu=(y+mean(y))/2#设计矩阵X=model.matrix(~x1+x2+x3+x4)#迭代运算repeat{

W=diag(mu)

z=log(mu)+(y-mu)/mu

beta

=

solve((t(X)%*%W%*%X))%*%t(X)%*%W%*%z

eta=c(X%*%beta)

mu=exp(eta)

D1=D0

di=ifelse(y==0,2*mu,2*(y*log(y/mu)-(y-mu)))

D0=sum(di)

diff=D1–D0

if(abs(diff)<1e-8)

break

}泊松回归的迭代加权最小二乘法#输出参数估计值beta#输出参数估计值:

(Intercept)-1.9085x10.4642x2-0.6805x310.3721x41-0.3267#应用glm估计模型参数mod1=glm(y~x1+x2+x3+x4,data=dt,family

=

poisson(link

=

log))summary(mod1)#泊松回归##

##Call:

##glm(formula=y~x1+x2+x3+x4,family=poisson(link=log),

##data=dt)

##

##DevianceResiduals:

##Min1QMedian3QMax

##-2.4119-0.7377-0.5795-0.41754.2155

##

##Coefficients:

##