R语言gmm模型

转载

墨色天香 2024-11-26 18:40:08

文章标签 R语言gmm模型深度学习样条数据建模 文章分类 R语言后端开发

前言

最近一段时间在做流行病数据的分析，期间学习DLNM模型的过程中碰到了挺多问题百度搜不到，笔者通过慢慢读原文献、看代码解决了一部分，当然还是有一些没太搞懂hhh。DLNM模型是个比较新的模型，中文版本的介绍也比较少，所以想写一篇推文，做一下知识输出，这也是笔者第一篇关于专业知识的推文，欢迎大家点赞、评论、多多支持，文中言语不当之处还请多多指教，谢谢！

1. 分布滞后模型与自回归模型

分布滞后线性模型最初由 Almon 于1965年提出，并应用于经济学研究。2000年，shcwartz 和 Braga 等人将该模型引入环境健康效应的定量化评估。其原理是用解释变量第t至t-s期的值来预测被解释变量的表现，其表达式如下：

与分布滞后模型不同的是，自回归模型的解释变量和被解释变量是同一个变量（不要弄混啦），其原理是用变量第t-1至t-q期的值来预测本期的表现，并假设它们为一线性关系，其表达式如下：

2. 分布滞后非线性模型

2.1 DLNM的历史及应用场景

过去，分布滞后线性模型（DLM）普遍应用于研究空气污染的健康效应，但由于该模型假设暴露-反应呈线性关系，而现实情况中很多暴露-反应关系（比如：温度-死亡）呈现U、V等非线性关系，导致分布滞后模型并不适用。针对上述问题，2006年 Armstrong 最早将分布滞后非线性模型（DLNM）提出并应用于流行病学当中。2010年 Gasparrin、 Armstong 等人进一步说明了如何利用交叉基（cross—basis）对解释变量进行特征变换，进而以广义线性模型和广义相加模型等传统模型的思想为基础，阐述了分布滞后非线性模型的理论。2012年 Gasparrin、 Armstong 等人介绍了如何将 DLNM 模型框架如何应用于生物医药领域的研究设计和队列研究，研究长期暴露于药物、治癌剂等环境下对健康的影响。2013年，Gasparrin 和 Armstong 继续拓展 DLNM 模型，提出了一种改进的方法来研究 DLNM 复杂的非线性和延迟关联，其中涉及将二维拟合简化到一维。2016年，Gasparrin 通过模拟比较了移动平均模型和分布滞后模型。2017年，Gasparrin、Armstong 等人将惩罚样条的理论引入了 DLNM 模型中。

2.2 DLNM模型理论

本文暂时只对DLNM的理论进行简单解释（笔者学识尚浅，在学习模型时感觉理论的学习非常依赖 R 语言中 DLNM 包的介绍，一边写程序一边读原文献会对学习模型有很大帮助），分布滞后非线性模型公式如下：

其中，代表时间，，是连接函数族，为多种可选概率分布（gammma、Poission等），因变量，，，通常是整形变量（如每日死亡人数，每日患病人数等），自变量，，，通常是同期的空气污染物浓度、温度、相对湿度、药物剂量等环境因素，表示其他混杂因素的线性效应，、为相应的系数。表示各种解释变量的交叉基函数，即对自变量与因变量的关系、滞后效应的分布分别选择合适的基函数，常用的基函数有正交函数、线性阈值函数和样条函数等，具体理论请参看文献。（笔者这部分看原文献没完全搞清楚，后续会继续研究）

虽然作者在原文献中说 DLNM 的代数理论非常复杂，但感觉复杂的部分实际上就在交叉基那里，整体的建模流程简单来说就两步：

利用交叉基进行特征转换，得到新的数据。
用广义线性模型对变换后的数据进行建模。

3. ChicagoNMMAPS数据集介绍及可视化

本文使用ChicagoNMMAPS数据集进行建模分析演示，该数据集包含1987-2000年芝加哥每日的死亡率、天气（温度、露点温度、相对湿度）和污染（PM10和臭氧）等情况。

date: Date in the period 1987-2000.
time: The sequence of observations
year: Year
month: Month (numeric)
doy: Day of the year
dow: Day of the week (factor)
death: Counts of all cause mortality excluding accident
cvd: Cardiovascular Deaths
resp: Respiratory Deaths
temp: Mean temperature (in Celsius degrees)
dptp: Dew point temperature
rhum: Mean relative humidity
pm10: PM10
o3: Ozone

4. DLNM模型的R语言实现

建模之前首先做一下数据清洗，这里笔者将离群值赋值为NA，然后使用mice包对缺失值进行了插补，这里可以参考知乎中的：“R中数据缺失值的处理--基于mice包”。

# 调用包
library(dlnm);library(ggplot2);library(mice);library(reshape2)

# 调用芝加哥数据
data = chicagoNMMAPS
data = as.data.frame(data)

# 数据预处理
data[data$death>200,'death'] <- NA  # 离群值赋值为NA
data[data$cvd>100,'cvd'] <- NA
data[data$resp>200,'resp'] <- NA
data$temp_dif[2:nrow(data)] <- diff(data$temp)  # 生成两日间温差序列
miceMod <- mice(data)  # 进行mice插值
data <- complete(miceMod)  # 生成完整数据
data[data$pm10>150,'pm10'] <- NA  # 这里比较奇怪，pm10中有缺失值所以不能赋值为NA，插补后就可以了，没太弄懂哈哈哈
miceMod <- mice(data)  # 进行mice插值
data <- complete(miceMod)  # 生成完整数据

清洗完数据后，我们先看一下ChicagoNMMAPS数据集大体长什么样子,图片及程序代码如下所示。从Fig.2可以看出，死亡人数、温度和污染物都以12个月为周期，且冬季温度较低、污染较严重、死亡人数较多，夏季呈相反趋势。

#### 数据可视化 ####
plot_varible <- c("death", "cvd","resp","temp","temp_dif","pm10","pm10")
plot_data <- melt(data,id=c('date'), measure.vars = plot_varible)
plot_data$variable <- factor(plot_data$variable, levels = levels(plot_data$variable), labels = c('死亡人数', '心血管疾病致死', '呼吸道疾病致死', '平均温度', '两天间温度差', expression(paste(PM[10], " (", mu, g, "/", m^3, ")", sep = "")), expression(paste(CO, " (", m, g, "/", m^3, ")", sep = "")))
)

work_dir <- "D:/Latex/Latex_project/Markdown/DLNM推文"
setwd(work_dir)
png(file = "chicagoNMMAPS数据集可视化.png", width = 4000, height = 6000, res = 300)

ggplot(plot_data,aes(x = as.factor(date), y = value, fill=variable, color=variable)) +
  geom_point(alpha=0.5,shape=21) +
  scale_y_continuous(name = "Values of varibles")+
  scale_x_discrete(name = "Time", breaks=c('1988-01-01','1989-01-01','1990-01-01','1991-01-01','1992-01-01','1993-01-01','1994-01-01','1995-01-01','1996-01-01','1997-01-01','1998-01-01','1999-01-01','2000-01-01','2001-01-01')) +
  theme_bw() +
  theme(legend.position = "none")+
  facet_wrap(~variable,ncol=1,scales='free', labeller = label_parsed)+
  theme(axis.title = element_text(face="bold"),
        axis.title.x = element_text(colour="black",size=20),
        axis.title.y = element_text(colour="black",size=20),
        strip.background = element_blank(),
        panel.border = element_rect(colour = "black"),
        legend.position = "none",
        axis.text.y = element_text(colour="black",size=12),  
        strip.text = element_text(size=16),
        axis.text.x=element_text(color = "black", size=10, vjust=.8, hjust=0.8)) 
dev.off()

R语言gmm模型_样条

Fig.2 chicagoNMMAPS数据集可视化

通过观察可以得知温度与死亡人数有很大的相关性，根据文献[1]的建模流程，先画一下温度和发病数的散点图看一下不同温度下的发病情况。

png(file = "温度散点图.png", width = 1500, height = 1000, res = 300)
ggplot(data,aes(x = temp, y = death, color='#99FFCC')) +
  geom_point(alpha=0.4)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_y_continuous(name = "Death")+
  scale_x_continuous(name = "Temperature")
dev.off()

R语言gmm模型_深度学习_02

Fig.3 温度-死亡人数散点图

整体来看，高温和低温的情况下死亡人数较多，数据整体呈现一个V形趋势。因为数据点太多了并不直观，所以笔者将数据集中温度的取值范围划分为200个等间隔的小区间，并在每个区间内求均值，绘制了下面的气泡图。气泡的大小和颜色的深浅反映了所在区间内数据点的数量。

scatter_plot <- function(input1, piece) {
  input = data[,paste(input1,"",sep="")]
  m = max(input, na.rm = TRUE)
  n = min(input, na.rm = TRUE)
  output=c()
  for (i in ceiling(piece*n/m):(piece-1)) {
    y_rate =  mean(data[m*(i)/piece < input & input < m*(i+1)/piece, 'death'] , na.rm = TRUE)
    number = length(data[m*(i)/piece < input & input < m*(i+1)/piece, 'death'])
    a <- cbind((i/piece)*m+m/(2*piece), y_rate, number)
    output = rbind(output,a)
    
  }
  output <- as.data.frame(output)
  colnames(output) = c(input1,'death','number')
  # print(output)
  ggplot(output, aes(x=output[,input1], y=death, size=number, color=number)) +
    scale_color_gradient2(low = "#deebf7", mid= '#6baed6',high = "#08306b") +
    geom_point(alpha=0.7) +
    theme(legend.position = "bottom")+
    xlab(input1) + 
    theme_bw()+
    scale_size_continuous(range = c(0.5, 5))
}

png(file = "温度均值散点图.png", width = 1500, height = 800, res = 300)
scatter_plot('temp',200)
dev.off()

R语言gmm模型_数据_03

Fig.4 温度-死亡人数气泡图

观察完温度-死亡人数的关系后，我们首先使用 DLNM 包中的 crossbasis 函数对原始温度数据进行转换。

x 通常是一列有序的时间序列向量
lag 表示最大滞后时长或滞后范围，通常根据 QAIC、QBIC 准则进行选择
argvar、arglag 分别用于生成用于横向特征维度和纵向滞后维度的两个基矩阵
group 在数据中分为多组的情况下会用到，比如用一个DLNM模型同时研究北上广深四个城市长期暴露在高污染情况下对人体健康情况的影响

其中 argvar、arglag 中基矩阵的参数选择，首先要根据数据情况选择不同的 fun 参数，进而选择不同 fun 对应的参数，我们来看几个例子。(为了方便理解，我们先假设滞后维度上都选取线性，最大滞后时长设置为10，只变换特征维度的基函数)

首先选取bs（B-样条），不了解 B-样条的同学可以参考知乎中的“简单粗暴：B-样条曲线入门”，这里 degree 代表每一段样条的最大阶数，knots 代表样条节点，我们首先将 degree 设置为 2 ，knots 设置为 20（在温度为 20°C 的时候断一下），因为只有芝加哥一个城市所以不用设置 group。

cb.temp = crossbasis(data$temp, lag=10, argvar=list(fun="bs", degree=2, knots= c(20)),
                    arglag=list(fun="lin"))

可以看到 crossbasis 将原本的一列温度数据转换成了6列，其中v1，v2，v3分别表示截距项、一阶项、二阶项，l1、l2分别表示样条断开的两部分，因为通过观察在节点两侧温度和死亡人数基本都是线性关系，所以v3.l1和v3.l2两个二阶部分基本都变成了0

R语言gmm模型_数据_04

Fig.5 温度-交叉基转换

类似的我们也可以根据数据特征选取自然样条 ns、惩罚样条ps、多项式函数 poly、线性函数 lin、阈值函数 thr、地层函数 strata 等进行基矩阵变换。不同函数的详细介绍参考 DLNM包的说明，以及 Generalized Additive Models 书籍中的介绍。提取码：1234

# 自然样条 ns，knots与bs样条类似，df决定了自然样条的自由度
cb.temp = crossbasis(data$temp, lag=10, argvar=list(fun="ns", degree=2, knots= c(20)),
                    arglag=list(fun="lin"))

# 多项式函数 poly，只需要选取最高的阶数 degree
cb.temp = crossbasis(data$temp, lag=10, argvar=list(fun="ns", degree=2, knots= c(20)),
                    arglag=list(fun="lin"))

# 阈值函数 thr， 需要选取双侧或单侧阈值
cb.temp = crossbasis(data$temp, lag=10, argvar=list(fun="thr", thr.value=c(19,23), side= 'd'),
                    arglag=list(fun="lin"))

# 地层函数 strata，选取不同的断点，断点间的解释变量取值是常熟
cb.temp = crossbasis(data$temp, lag=10, argvar=list(fun="strata", breaks=c(19,23)),
                    arglag=list(fun="lin"))

本次建模使用了特征维度使用了ns，滞后维度使用了poly。根据数据特征选择完交叉基后需要使用转换得到的新特征拟合模型，本次建模我们选取类 possion 模型。原文献说最大滞后时长的确定根据选取最小的QAIC、QBIC准则。然而QAIC、QBIC原本是用于类 poisson 模型，类 poisson 模型通常不涉及滞后，很多情况下滞后时长越长，QAIC、QBIC越小，但是太长的滞后时长可能并不符合现实情况。所以笔者根据局部最小 QAIC、QBIC 选取最大滞后时长为30天。

library(splines) # 使用样条要用到这个package
cb.temp = crossbasis(data$temp, lag=30, argvar=list(fun="ns", knots= c(20)),
                     arglag=list(fun="poly",degree=4)) # 本次建模使用了特征维度使用了ns，滞后维度使用了poly
model1 = glm(death ~ cb.temp + ns(time,10*1)，family=quasipoisson(), data) # 拟合类 poission 模型
summary(model1)
pred1.temp = crosspred(cb.temp, model1, cen=round(median(data$temp)), bylag=0.2) # 拟合模型计算 RR 以及 RR 的置信区间， cen 选取对照点作为 RR 分母，bylag 表示 lag 的切片单位。

# 根据 QBIC、QAIC 选取最大滞后时长
qaicbic <- function(model) {
  phi <- summary(model)$dispersion
  logll <- sum(dpois(ceiling(model$y), lambda=exp(predict(model)), log=TRUE))
  cbind((-2*logll + 2*summary(model)$df[3]*phi),
        (-2*logll + log(length(resid(model)))*phi*summary(model)$df[3]))}

tqba <- data.frame()
for (ii in 7:100) {
  cb.temp = crossbasis(data$temp, lag=ii, argvar=list(fun="ns", knots= c(20)),
                       arglag=list(fun="poly",degree=4))
  model1 = glm(death ~ cb.temp + ns(time,10*1), family=quasipoisson(), data)
  summary(model1)
  q <- cbind(ii,qaicbic(model1))
  tqba <- rbind(tqba,q)
}
print(tqba)

接下来我们通过2D和3D两种展示形式来看一下不同温度、滞后下的相对危险度，3D图的话更立体，通过3D图看某一2D截面的话趋势会比较清晰，但是3D图只能展示某一个角度，会出现某些部分看不清楚的情况，2D等高线图的话可以看清整个温度-滞后范围下的RR值。

png(file = "single_temp.png", width = 3000, height = 2500, res = 300)
plot(pred1.temp, "contour", xlab="MeanTemp", key.title=title("RR"),cex.axis=2,
     plot.axes={axis(1,cex.axis=2)
       axis(2,cex.axis=2)},
     key.axes = axis(4,cex.axis=2),
     plot.title=title(xlab="MeanTemp (°C)",ylab="Lag (days)",cex.main=2,cex.lab=1.5))
dev.off()

R语言gmm模型_R语言gmm模型_05

可以看出来，高温和低温情况下的 RR 都是比较高的，且高温情况下 RR 会随着滞后时长变长快速下降，影响较短，低温情况下 RR 下降的更加缓慢，对人的影响更加持久。

png(file = "single_temp_3D.png", width = 3000, height = 2500, res = 300)
plot(pred1.temp,ticktype='detailed',border='#3366FF',xlab="MeanTemp (°C)",ylab="Lag (days)",zlab="RR",col='#99FFCC',shade = 0.1,cex.lab=1.3,cex.axis=1.3,lwd=1,theta = 20, phi = 25,ltheta = -35)
dev.off()

R语言gmm模型_样条_06

可以看出过去3天的累积 RR 在高低温情况下都比较大，因为低温比高温对死亡率的影响更加持久，随着时长的增加，高温环境下的危险度逐渐降低。

png(file = "累积效应.png", width = 4500, height = 1600, res = 300)
par(mfrow=c(1,3))
crall <- crossreduce(cb.temp,model1,cen=20,type="overall",lag=c(1,3))
plot(crall,xlab="Meantemp",ylab="RR",col=2,lwd=2,cex.lab=1.2,cex.axis=1.2,mar=c(1,2,0,1))
mtext(text="Overall cumulative association 1-3 days",cex=0.89)

crall <- crossreduce(cb.temp,model1,cen=20,type="overall",lag=c(1,7))
plot(crall,xlab="Meantemp",ylab="RR",col=2,lwd=2,cex.lab=1.2,cex.axis=1.2,mar=c(1,2,0,1))
mtext(text="Overall cumulative association 1-7 days",cex=0.89)

crall <- crossreduce(cb.temp,model1,cen=20,type="overall",lag=c(1,30))
plot(crall,xlab="Meantemp",ylab="RR",col=2,lwd=2,cex.lab=1.2,cex.axis=1.2,mar=c(1,2,0,1))
mtext(text="Overall cumulative association 1-30 days",cex=0.89)
dev.off()

R语言gmm模型_R语言gmm模型_07

Fig.8 累积效应

5. 待讨论的问题

原文献中代数表达式的理解问题
DLNM包中源代码的理解问题
滞后单位不同时的最大滞后时长选择问题
QAIC、QBIC准则的拓展研究
惩罚样条的理解及使用方法
分布滞后非线性模型灵敏度研究
多元情况下的组变量选择问题
分布滞后非线性模型与多元时间序列分析之间的联系

总结

很早之前就想写一篇 DLNM 的推文了，但是由于事情比较多一直也没来得及写。这是第一次在知乎发推文，虽然花了不少的时间和精力，但是在知识输出的过程中还是学习了不少新的知识和技能，发现了很多之前没理解的问题，也有了很多新的 idea。写文章的过程还是非常开心的，以后会不断做知识输出~

参考文献

[1] 杨军, 欧春泉, 丁研,等. 分布滞后非线性模型[J]. 中国卫生统计, 2012, 29(005):772-773,777.

[2] Armstrong B. Models for the relationship between ambient temperature and daily mortality[J]. Epidemiology, 2006: 624-631.

[3] A, Gasparrini, B, et al. Distributed lag non-linear models[J]. Statistics in Medicine, 2010.

[4] Gasparrini, Antonio. Reducing and meta-analysing estimates from distributed lag non-linear models.[J]. Bmc Medical Research Methodology, 2013.

[5] Gasparrini A , Armstrong B . Reducing and meta-analysing estimates from distributed lag non-linear models[J]. Bmc Medical Research Methodology, 2013, 13(1):1-1.

[6] Gasparrini A . Modelling Lagged Associations in Environmental Time Series Data: A Simulation Study[J]. Epidemiology, 2016, 27.

[7] Antonio, Gasparrini, Fabian, et al. A penalized framework for distributed lag non-linear models.[J]. Biometrics, 2017.

本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。