R语言GAMLSS模型对艾滋病病例、降雪量数据拟合、预测、置信区间实例可视化

GAMLSS模型是一种半参数回归模型,参数性体现在需要对响应变量作参数化分布的假设,非参数性体现在模型中解释变量的函数可以涉及非参数平滑函数,非参数平滑函数不预先设定函数关系,各个解释变量的非线性影响结果完全取决于样本数据。

由Kaizong Ye,Colin Ge撰写

它克服了GAM模型和广义线性模型(Generalized Linear Models, GLM)的一些局限性。


对连续分布数据拟合的实例–降雪量数据

降雪:63年的年降雪量,每年降雪量数据

目的:帮助客户证明连续分布对单个变量的拟合。结论:正态假设是适当的。

模型的拟合和显示

数据集是降雪数据,数据显示,63年的年降雪量。


> names(parzen)

在这里,我们将数据拟合为正态分布(NO)、γ(GA)、幂指数(PE)分布。

正态与伽马的比较探讨了数据中是否存在正偏性。正态与幂指数的比较表明了峰度的可能性,而BCPE则显示出数据中是否同时显示了偏度和峰度。GAIC将帮助我们在不同的分布之间进行选择。


视频

非线性模型原理与R语言多项式回归、局部平滑样条、 广义相加模型GAM分析

探索见解

去bilibili观看

探索更多视频


> > mBCPE <- histDistsnowfall, "BCPE", density = TRUE, main = "(d)", +
QQ截图20230327124351.png

请注意,选项密度=true请求将非参数核密度估计包含在图中

> GAIC
image.png

GAIC()函数的默认惩罚是k=2,Akaike信息准则(注意,我们可以使用等价函数AIC())。AIC准则表明,正态分布与数据完全吻合。图显示了四个不同的分布。

检验模型

使用R函数ks.test()提供的Kolmogorov-Smirnovness拟合测试测试正态模型(或任何其他模型)的充分性在这里是不可取的,因为我们必须估计分布参数u和o,所以测试无效。 (归一化分位数)残差的检验将提供一种研究适配适足性的方法。归一化分位数残差是独立的标准正态变量。我们期望拟合的(归一化分位数)残差I;近似地表现为正态分布的变量(即使最初的观测值Y不一定是正常的),因此残差的归一化Q-Q图在这里是合适的。r软件提供了用于绘制QQ-绘图的函数。

image.png
image.png

图片

【视频】广义相加模型(GAM)在电力负荷预测中的应用

阅读文章


检验分布拟合参数可靠性的方法有两种:1)汇总函数和Vcov函数。一般来说,这两个值应该是相同的,因为在默认情况下,汇总是vcov获得的标准误差。Vcov()得到的标准误差是通过反演全观测信息矩阵得到的,它们考虑了分布参数估计之间的关系。


随时关注您喜欢的主题


注意,vcov()函数再一次修改最后的模型,以获得Hessian矩阵。

我们修改了所选择的最终模型

> moNO <- gamls
image.png
> summary(moNO)
image.png
> vcov(modNO, type = "se")
image.png

拟合模型由Y~NO(i,a)给出,其中ji=bo=80.3,log()=PO2=3.158,因此6=23.52。注意,j和o是u和o的极大似然估计。

使用vcov()结果,log(O)=Bo2的95%置信区间(CI)为[3.158-(1.960.08922),3.158+(1.960.08922)]=(2.983,3.333),由此[exp(2.983),exp(3.333)]=(19.75,28.02)给出了o的95%CI置信区间。可以与图中的剖面偏差区间(19.96,28.32)进行比较,得到了用下列R脚本得到的[exp(3.021),exp(3.33)]=(20.51,27.93)所给出的自举CI。

image.png
> library(boot)  
>
> funB <- function(data, i) {  
+ d <- dtaframe(swfall = data[i, ])  
+ coef(updae(modNO, dat = d), "sigma")  
+ }  
> boot(paren, funB R  199))
image.png
> plot(modOboot)
> boot.ci
image.png

艾滋病病例数据

在这里,我们使用季度报告的艾滋病病例组成的数据,这些数据来自传染病监测中心公共卫生实验室服务部门。我们首先帮助客户使用泊松族来建模报告病例的数量(响应变量),针对时间(一个连续的解释变量),我们用一个三次样条平滑器,使用5有效自由度,针对Qrt,一个代表季度季节性效应的因子。

然后(I)将该族转换为负二项式(I型)(II),用df=8(Ji)更新平滑参数,去掉季度季节效应(Iv),最后拟合一个响应log(Y)的正态族模型。


GAMLSS-RS iteration 1: Global Deviance = 448.1315 GAMLSS-RS iteration 2: Global Deviance = 448.1315 > # 用负二项分布 > h.nb <-update(h.po, family=NBI) GAMLSS-RS iteration 1: Global Deviance = 388.8086 GAMLSS-RS iteration 2: Global Deviance = 390.7803 GAMLSS-RS iteration 3: Global Deviance = 391.396 GAMLSS-RS iteration 4: Global Deviance = 391.3996 GAMLSS-RS iteration 5: Global Deviance = 391.3965 GAMLSS-RS iteration 6: Global Deviance = 391.3962 > # 用平滑项 > h.nb1 <-update(h.nb,~cs(x,8)+qrt) GAMLSS-RS iteration 1: Global Deviance = 362.943 GAMLSS-RS iteration 2: Global Deviance = 359.1257 GAMLSS-RS iteration 3: Global Deviance = 359.229 GAMLSS-RS iteration 4: Global Deviance = 359.2342 GAMLSS-RS iteration 5: Global Deviance = 359.2348 GAMLSS-RS iteration 2: Global Deviance = -42.3446

预测

使用函数也可以提取模型中特定分布参数在解释变量当前数据值处的线性预测、拟合值和特定项。

为模型中的任何参数提供概要偏差,而不仅仅是对于分布参数。


> predict(aids.1) 1 2 3 4 5 6 7 8 1.322524 1.490931 1.996051 2.140244 2.540856 2.64334

让我们首先假设我们有兴趣拟合一个线性的时间项(X)加上季度季节效应的一个因子,QRT,使用负二项式模型(I型)家族。这个模型被拟合为


GAMLSS-RS iteration 1: Global Deviance = 492.7247 GAMLSS-RS iteration 2: Global Deviance = 492.6375 GAMLSS-RS iteration 3: Global Deviance = 492.6373 > summary(aids1)
image.png

线性时间项(Z)的系数为0.08744,其t值表明它是高度显着的。使用0.08744+-1.96x 0.005382可以得到该参数的95%置信区间,从而得到近似区间(0.0769,0.0980)。现在,我们将使用函数Pror项来为线性项参数找到一个更精确的95%置信区间。请注意,模型公式中的此值指示要配置文件的参数。

prof.term(mo
image.png
image.png
image.png

如我们所见,95%的剖面偏差置信区间为(0.0746,0.1009),与近似的(0.0769,0.0980)相差不远。


可下载资源

关于作者

Kaizong Ye拓端研究室(TRL)的研究员。

本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。

​非常感谢您阅读本文,如需帮助请联系我们!

 
QQ在线咨询
售前咨询热线
15121130882
售后咨询热线
0571-63341498

关注有关新文章的微信公众号


永远不要错过任何见解。当新文章发表时,我们会通过微信公众号向您推送。

技术干货

最新洞察

This will close in 0 seconds