DLNM时间序列数据分析（2017年1月版）“Example 2

翻译自：《Distributed lag linear and non-linear models for time series data》（Antonio Gasparrini， 2017-01-16）

（译者：分布滞后非线形模型是近几年兴起的时间序列数据研究模型，目前国内对该模型的研究还不多。我在使用该模型研究健康问题时遇到了很多问题，对于模型的理解和参数的设置有很多不清楚的地方。R包DLNM的作者Antonio Gasparrini撰写了大量的文章，对该包的使用进行了说明，我将其2017年的一篇文章中的部分内容进行了翻译，希望与学习和使用该包的同学共同探讨）

---

第二个例子用于展示如何分析特定季节数据。这一问题的特点在于数据是由不同年份的多个等距的有序季节数据组成的，并非一个单一的了连续序列。在这个例子里，我将使用和section3中类似的步骤，分别评估臭氧和温度对于死亡5至10日的滞后影响。

首先，我创建了一个限制在夏季（6月～9月）的季节性的时间序列数据，并把它以data frame 形式保存在 chicagoNMMAPS 里：

chicagoNMMAPSseas <- subset(chicagoNMMAPS, month %in% 6:9)

第一次创建交叉积矩阵

cb2.o3 <- crossbasis(chicagoNMMAPSseas$o3, lag=5,
       argvar=list(fun="thr",thr=40.3), arglag=list(fun="integer"),
       group=chicagoNMMAPSseas$year)
cb2.temp <- crossbasis(chicagoNMMAPSseas$temp, lag=10，
       argvar=list(fun="thr",thr=c(15,25)),
       arglag=list(fun="strata",breaks=c(2,6)), 
       group=chicagoNMMAPSseas$year)

“group”指明分组变量：计算机函数会在每一组的末端打断序列，并将第一行（first row）替换为交叉积矩阵的最大滞后值，其后的季节则填充NA。每个系列应为连续，完整和有序的。我假设臭氧（O3）的效应在低于40.3 μgr/m3 时无影响，之后则为线性，代码中设置为高阈值参数化（high threshold parameterization）（fun="thr"）。温度方面，我选择了双阈值，即假设其效应在低于15◦C和高于25◦C时呈线性，在15◦C至25◦C间时无效应。阈值使用语句thr.value设置（缩写为thr），参数side在第一个交叉积中默认值为“h”，第二个交叉积中默认值为“d”（在制定双阈值的时候使用）。关于滞后维度，我为臭氧定义了一个非结构方程（unconstrained function），对每个滞后使用一个参数设置(fun="integer") ，最大滞后天数为5天（默认的最小滞后值为0）。我为温度变量定义了一个3层区间，分别滞后0-1天，2-5天和6-10天。所有交叉积的设置可以通过语句 summary()查看。

为研究年份中的日和时间构建包含自然样条回归模型，分别描述每年的季节效应和长期趋势。与前者相比，后者的自由度要少很多，因为其只需要捕捉年均趋势。除此之外，其他部分与section3中步骤类似。代码如下：

model2 <- glm(death ~  cb2.o3 + cb2.temp + ns(doy, 4) + ns(time,3) + dow,
          family=quasipoisson(), chicagoNMMAPSseas)
pred2.o3 <- crosspred(cb2.o3, model2, at=c(0:65,40.3,50.3))

必须计算的预测值在at语句中说明，本例中，我定义了从0至65μgr/m3 的整数（这是臭氧分布的大致范围），另外增加了阈值和阈值增加10个单位的值。向量是自动排序。参考通过语句thr()自动选择暴露--反应模型曲线，参数cen也为默认值。

我绘制了10个臭氧单位的给定值的滞后效应预测关系曲线，与section 3中的类似，但本例添加了80%置信区间，以及累计暴露--反应关系，代码如下：

plot(pred2.o3, "slices", var=50.3, ci="bars", type="p", col=2, pch=19,
    ci.level=0.80, main="Lag-response a 10-unit increase above threshold (80CI)")
plot(pred2.o3,"overall",xlab="Ozone", ci="l", col=3, ylim=c(0.9,1.3), lwd=2, 
    ci.arg=list(col=1,lty=3), main="Overall cumulative association for 5 lags")

代码第一行中，参数 ci="bars" 指定了置信区间以bars的方式绘制，这点与默认不同，默认值为“area”。参数 ci.level=0.80 指定了80%置信区间需绘制出来。我选择了“点”而不是“线”绘制，“点”的设置参数为 type 和 pch。
代码的第二行，参数 type="overall" 表示需绘制带有置信区间的累计估计曲线，ylim定义y轴数值范围，lwd定义线条宽度。本例中置信区间以线条的形式显示，这是通过给参数 ci 设置 "l" 来实现的。与前例类似，置信区间曲线通过参数ci.arg精细化处理。

Figure 2
与前例类似，我们可以从 pred2.o3 中提取臭氧阈值之上10个单元的带95%置信区间的累积效应估计值：

pred2.o3$allRRfit["50.3"]

显示为：

    50.3
1.047313

置信区间：

cbind(pred2.o3$allRRlow, pred2.o3$allRRhigh)["50.3",]

显示为：

[1] 1.004775 1.091652

同样的图形和计算可应用于冷及热效应。举个例子，我们可以描述温度高低阈值之外每升高1◦C的风险，读者可以参照上面的步骤实现。