基于时空变系数模型的台湾呼吸系统疾病研究
陈天雷，中南大学数学与统计学院统计学2020级本科生，中国人民大学统计学院流行病与卫生统计学2024级研究生。
今天跟大家分享的是一篇发表于2022年Biostatistics的文章，文章提出了一个时空疾病制图模型，以研究台湾地区与呼吸系统疾病医院就诊率之间的关联。
Wang, F., Duan, C., Li, Y., Huang, H., Shia, B.C. (2022). Spatiotemporal Varying Coefficient Model for Respiratory Disease Mapping in Taiwan. Biostatistics. DOI: 10.1093/biostatistics/kxac046
一、引言
在公共卫生领域，呼吸系统疾病一直是一个全球范围内的重要问题。最近，由于(直径小于2.5 m的颗粒物）的普遍存在，有关其与呼吸系统疾病关系的研究引起了广泛关注。有文献研究了长期和短期暴露于对健康的影响。例如，Tumer和其他人(2011年)研究了美国非吸烟者长期平均环境浓度与肺癌死亡率之间的关系。在台湾，一些新出现的研究也验证了对呼吸系统疾病的影响。例如，长期暴露于环境与慢性阻塞性肺病(COPD)的风险增加有关(Guo等人，2018年)；非吸烟者也有类似发现(Huang等人，2019年)。
在我们的工作中，我们研究了台湾地区与一般呼吸系统疾病医院就诊率之间的空间变化关系，为的危害提供了更全面的证据。探究呼吸系统疾病流行情况的一个常用方法是疾病制图模型，该模型在公共卫生和流行病学领域有着悠久的历史(Banerjee等,2014年；Wang等，2019年)。通过一些探索分析，我们发现对医院就诊率的影响存在时空差异。因此，要研究对台湾呼吸系统疾病发病率的影响，考虑空间和时间变化的疾病制图模型非常有意义。
在文献中，已有模型可以描述风险因素影响的空间异质性。然而，文献中并未考虑“change of support problem”(Gelfand等人，2001年)，即空间观测数据可以在不同的空间分辨率下收集。具体来说，疾病数据(如呼吸系统疾病的医院就诊人数)是以区域为单位进行观测的，称为区域数据。协变量数据(如浓度)通常是在一组监测站记录的数值，称为点数据。为了研究区域层面的疾病结果与点层面协变量之间的关系，我们需要首先将点层面的协变量重新调整到区域层面。以往文献中常见的做法是采用插值法(如克里金法)预测每个区域中心点的协变量值。然而，这种预测会给随后的疾病制图模型带来新的不确定性。
为了解决上述问题，我们在本研究中提出了一个时空变系数模型，用于处理空间错位的疾病结果(即台湾328个地区的呼吸系统疾病医院就诊率)和协变量数据(即台湾69个监测站记录的浓度)。第一步，我们将点级浓度视为具有Matérn空间方差函数的连续高斯过程。然后，应用随机偏微分方程 (Stochastic Partial Differential Equation, SPDE) 方法 (Lindgren 等，2019年)，通过有限基函数将连续高斯过程近似为离散高斯马尔可夫随机场。因此可以预测未观测地点的浓度。接着，就可以利用蒙特卡罗积分来计算的区域浓度。第二步，我们提出了一种新的时空变系数模型，用于分析区域级浓度与呼吸系统疾病发病率之间的关联。在所提出的模型中，空间变系数采用固有条件自回归 (Conditional Autoregression, CAR)模型，而时间变化则采用参数趋势模型。最后，为了考虑预测的不确定性，进行了多次蒙特卡罗积分，然后对具有变化系数的时空疾病制图模型进行了迭代估计。
由于提出的时空模型是在贝叶斯框架下建立的，因此模型拟合通常采用马尔科夫链蒙特卡罗(MCMC)方法。然而，MCMC 通常计算成本高、耗时长，不适用于大规模数据集。特别是，为了考虑估计点级协变量的不确定性，需要进行多次蒙特卡罗积分，这又增加了计算负担。在这项工作中我们应用了集成嵌套拉普拉斯近似 (Integrated Nested Laplace Approximation, INLA)方法(Rue等，2009年)来提高计算效率，该方法常用于疾病制图模型。
总之，我们对这项工作做出了以下贡献。首先，我们提出了一个新的时空变系数模型来描述浓度对台湾呼吸系统疾病发病率的时空变化作用。具体而言，空间变化系数采用CAR模型，而时间变化则采用参数趋势模型。其次，为解决“change of support problem”，我们将点水平协变量视为连续高斯过程，然后应用SPDE方法求得其区域平均值。此外，为了考虑预测区域值的不确定性，我们进行了多次蒙特卡罗积分，然后对疾病制图模型进行迭代估计。最后，为了减轻计算负担，采用了INLA方法来提高计算效率。
二、数据
1.就诊率与PM2.5浓度
本文研究范围包括台湾的349个三级行政区。对于每个地区，我们从国民健康保险研究数据库(NHIRD，https://nhird.nhri.org.tw/)中收集了2012年所有医院不同疾病的每周就诊人数。在这项工作中，我们将所有医院的就诊次数视为风险人数。有28个地区的所有疾病的每周就诊人数为零，因此我们将这些地区合并到属于同一上级行政单位的邻近地区。这样，共有328个地区。
浓度(单位：)的原始数据来自台湾环境保护署(http://taqmepa.gov.tw)台湾的空气质量监测网络由69个监测站组成。浓度的原始值每小时记录一次。图1(a)显示了2012年所有地区各周的平滑平均趋势。平均趋势呈粗略的“T”形，表明呼吸系统疾病在寒冷季节(1月、2月和12月)比温暖季节(6月、7月和8月)更加严重，这一时间规律与图1(b)所示的的时间规律相吻合。我们在实证分析中将季节作为控制变量。
图2(a)显示了呼吸系统疾病医院原始就诊率的空间模式，其中显示了所有地区的年平均就诊率。具体而言，大多数地区的呼吸系统疾病原始就诊率介于0.1和0.4之间。台湾中部地区的就诊率通常较高。最高的就诊率为0.45，出现在位于台湾中西部的水上乡；而最低的就诊率正好为0，出现在面积37平方公里、人口55886人的大埔乡。由于大埔是一个人口较少的小地区，该地区在1或2周内偶尔出现没用病人的情况。然而，其真实比率可能会更高，因为该地区大多数社区的医院原始就诊率都非常高。
图1(b)显示了所有地区每周内插的，以及平滑平均曲线的时间趋势，台湾地区全年的浓度一般在10到60之间变化。此外，春季和冬季的浓度值高于夏季和秋季。图2(b)中展示了四季浓度的区域分布。从图中可以看出，台湾西南部地区的通常相对高于其他地区。
2.人口构成、气象数据、吸烟率以及医院密度
328个地区按5岁年龄组划分的人口构成数据来自中华民国内政部户政司https:/wwwrisgov.tw/app/portal/346)。为了描述异质性人口构成的特征，我们计算了每个行政区域的儿童 (0-14岁)比例。此外，我们还从台湾中央气象局(https://www.cwb.gov.tw/)收集了包括温度和相对湿度在内的气象数据。台湾共有537个气象监测站，其中193个提供温度数据，113个提供相对湿度数据。温度和相对湿度的原始值也是每小时记录一次。最后，我们还收集了吸烟率和医院数量这两个控制变量。
在我们的研究中，我们只掌握了台湾18个二级行政区的总体吸烟数据。这些数据来自健康促进署每年进行的“成人吸烟行为调查”。在这项调查中，吸烟率被定义为在过去30天内吸过100支以上香烟并使用过烟草制品的18岁以上人群的百分比。台湾的平均吸烟率为0.18。医院数量是指各地区诊所、地区医院、区域医院和医疗中心的总数。328个地区的医院数量数据来自台湾国家卫生与健康研究所。平均而言，328个地区的医院数量为126.25家。在本文中，我们使用医院密度(即医院数量除以地区总人口)作为控制变量。
三、方法
1.疾病制图模型综述
假设我们在一组区域上观测到疾病计数，为了建立的时空模型，我们假定了一个表征时空效应的简单加法模型：
其中，是区域在时刻协变量，是疾病计数，是风险人口总数，是真实疾病发病率，是相应的系数，是特定区域随机效应，是时间随机效应。另外，常被假设服从CAR分布，即：
其中，是控制空间相关性随距离衰减的比例参数，如果区域,相邻，则1，否则为0。在高斯假设下，通过一些简单的计算，我们很容易得出每个空间随机效应的条件分布，即：其中，。如(3.3)所示，的条件分布是正态分布，其均值是区域相邻所以的加权平均值，方差也受其邻居数量的影响。这样，CAR分布就能反映各地区的空间结构。2.系数变化的时空疾病制图模型
为了反映风险因素对疾病发病率总体和特定地区的影响，我们采用了变系数的时空疾病制图模型。如前所述，空间观测数据通常是在不同的空间分辨率下采集的。
以台湾呼吸系统疾病数据为例。每个地区收集到的呼吸系统疾病医院就诊人数是原始数据。然而，浓度是在台湾69个监测站收集的点数据。这就导致了“change of support problem”，即以不同空间分辨率收集的数据的融合问题(Gelfand等人，2001年)。为解决这一问题，我们采用了SPDE方法(Lindgren等，2019年)，该方法常用有限元法求解，得到协变量(如PM2.5浓度)在任意区域的预测值。但是考虑到预测值的不确定性，需要进行多次蒙特卡洛积分，即不断利用新的协变量值和网格点参数值进行迭代更新。
按照混合效应模型的标准设置，我们将具有固定效应的点级协变量记为，将具有固定效应的区域协变量记为。进一步，将具有随机效应的点级协变量记为，将具有固定效应的区域协变量记为。我们提出的模型如下：
其中，和是固定效应，反映了和在所有区域和时间的总体影响；和反映了和在第和区域的特定效应；是不能被和解释的空间随机效应，我们进一步假设它服从参数为的CAR分布(详见(3.2))。另外，我们用参数时间趋势来描述时空交互作用：时间效应由未知参数的线性项表示，交互作用由表示。我们用表示与中第个协变量对应的变系数，类似定义。假设它们服从参数为的CAR模型，与前述类似。接下来，我们考虑变系数时空疾病制图模型的估计问题，所有未知参数为。我们采用MCMC方法估计模型（计算细节可在生物统计学网站的补充材料附录A中找到），并进行多次蒙特卡洛积分解决预测点级协变量的不确定性。另外，数据具有高维的空间结构且SPDE方法也增加了计算负担，因此，我们采用INLA方法代替MCMC。MCMC方法试图估计所有潜变量和参数的联合后验分布，而INLA方法侧重于每个单独变量或参数的边际后验分布。INLA方法由R中的INLA软件包实现，该方法大大降低计算成本，在我们的工作中，一次蒙特卡洛积分预测区域值，INLA方法仅需3.5分钟，而MCMC方法需要约53.7分钟。
四、模拟研究
1.模拟设计
我们利用台湾的地理结构生成数据。，，每个时间段的两个协变量和。截距和相应的系数为，和。时空交互效应为，其中从参数为，的CAR先验分布中生成。最后，我们考虑协变量依赖性和不同类型的系数。
协变量依赖性。在实践中，协变量之间可能存在相关性。为了验证模型在依赖性协变量情况下的性能，我们考虑了三种协变量依赖情况，即不相关、弱相关和高度相关。
不同类型的系数。在整个模拟研究中，所有区域的截距都是固定的。至于和的系数类型，我们考虑了两种情况。第一种情况，假设固定为0.5，而假设在不同地区各不相同，它是由规模参数的CAR先验分布生成的。在第二种情况下，假设和在不同地区都是不同的。具体来说，我们从CAR分布中生成和，参数分别为和。
基于上述数值，我们根据(3.5)计算出真实疾病率。生成的在0.02和0.13间变化。然后，疾病计数从均值为的泊松分布中独立同分布生成。这里，独立同分布，服从均值为，方差为的正态分布。我们考虑了两种均值设置：，并在每个设置中重复生成过程100次。除上述实验设置外，我们还研究了点水平协变量情景下的估计性能，详情参见文章的补充材料附录B.1。另外，还对预测性能进行了研究，详情参见文章的补充材料附录B.2。
2.模拟结果
为了评估我们提出的模型的性能，我们将我们的模型与三个对象进行了比较，分别是：(i)具有固定系数的时空疾病制图模型；(ii)Bakar等人(2016年)提出的空间变化系数模型，其中变化系数由高斯过程建模；(iii)经典疾病制图模型，但不包含空间随机效应。对于不同模型，我们采用INLA进行估计，然后我们从三个方面对不同模型进行评估。
首先，对于参数，我们使用均方误差(MSE)和平均绝对相对误差(MARE)。另外，为了评估不同方法在量化不确定性方面的性能，我们计算了系数和疾病发生率的蒙特卡洛方差以及参数的标准误差(SE)。最后，我们比较不同模型的计算效率，记录平均计算时间。表1和表2列出了相关情况。表中“VC”、“FC”、“BK”和“NAS”分别表示我们的空间变系数模型、固定系数模型、Bakar的变系数模型和无空间随机效应的疾病制图模型。我们得出以下结论。
表1 固定、变化情况下的模拟结果以及平均计算时间第一，对于和的估计，我们的时空变化系数模型几乎在所有情况下都能获得最小的MSE和MARE值。第二，随着相关系数的增大，我们的模型获得的MSE和MARE值也相差无几。这一结果表明，我们提出的模型在处理高度相关的协变量时具有稳健性。此外，随着的均值(即)的增加，我们发现模型的MSE和MARE值呈下降趋势。第三，对于真正的固定系数参数(如），不同模型得到的估计值的方差都相对较小，这表明这些模型的不确定性很小。对于变化系数参数(如），我们提出的模型和Bakar的模型都比固定系数模型有更大的方差，因为前者可以对不同地区产生不同的系数估计值。尽管如此，我们提出的模型得到的估计方差还是小于Bakar的模型。第四，比较经验方差和估计方差，我们发现在我们提出的模型下，它们是相似的。然而，在其他模型下，由于模型规格不正确，它们的表现截然不同。
最后，我们关注不同方法所消耗的计算时间。一般来说，我们提出的模型和Bakar的模型都会消耗更多的时间。这是因为变化系数模型的结构更为复杂，需要估计的参数也更多。与Bakar的模型相比，我们的时空变化系数模型表现得更快，这意味着我们模型的计算效率更高。
五、台湾呼吸系统疾病
在将疾病制图模型应用于台湾的呼吸系统疾病数据之前，我们首先利用卡方检验法进行了初步分析，详见文章的补充材料附录C.1。通过卡方检验，我们发现呼吸系统疾病的医院就诊率与浓度之间存在显著的相关性。然后，我们应用所提出模型来描述在台湾的时空变化效应。请注意，浓度、温度和相对湿度都是点水平协变量。我们采用前述SPDE方法来预测它们的区域值，详情请参见文章的补充材料附录C.2。
接下来，我们讨论疾病制图模型的公式。我们共考虑了七个协变量。它们分别是浓度、温度、相对湿度、儿童比例、吸烟率和医院密度。为了使、温度和相对湿度的分布近似于正态分布，符合SPDE方法的假设，我们对这三个协变量进行了对数变换。在季节方面，我们以春季为基线，因此定义了三个虚拟变量，即夏季、秋季和冬季。
为了研究是否对每个协变量设置变系数，我们进行了一些初步分析，结果表明，只有浓度的相关系数具有明显的空间模式。因此，我们只假设浓度的系数是变化的，而其他协变量的系数都是固定的。此外，浓度是我们关注的主要协变量，我们假设它既有总体效应，也有偏离总体效应的效应。因此，我们假设具有固定效应和空间随机效应。我们的模型如下：
如前所述，考虑到参数集很大、预测区域值不确定等问题，我们使用INLA进行模型估计。在配备2.6GHz\ CPU和128Gb内存的台式机上，用INLA对模型进行一次蒙特卡洛积分，拟合浓度、温度和相对湿度的区域预测值大约需要3.5分钟。因此，要复制模型估计100次，INLA总耗时约为5.8小时。我们还使用R软件包CARBayesST(Lee等，2018年)应用了基于MCMC的方法。单次蒙特卡洛积分耗时约53.7分钟，100次蒙特卡洛积分的总耗时约为89.5小时。表3按INLA列出了各参数的估计后验均值、标准差和95%置信区间。
接下来，我们讨论所提出模型的系数估计情况。如表3所示。
表3 台湾呼吸系统疾病数据的模型估计结果(的固定效应)的后验均值为0.017，95%置信区间为(0.012,0.023)。这说明PM2.5浓度总体上与呼吸系统疾病医院就诊率呈正相关。(空间变化系数)表征了PM2.5对呼吸系统疾病医院就诊率的区域性影响。因此的值衡量了在区域内的总影响。在所有328个地区中，有213个地区的值为正。结果表明，在大多数地区，较高的浓度会导致较高的呼吸系统疾病发病率。此外，不同地区的实际发病率也不尽相同。
(温度的固定效应）的后验均值为-0.539，95%可信区间为(-0.552，-0.525)，(相对湿度的固定效应)的后验均值为-0.236，95%置信区间为(-0.251,-0.222)。结果表明，气象协变量温度和相对湿度都对呼吸系统疾病医院就诊率有显著负向作用。(儿童比例的固定效应)的后验均值为-1.737，相应的95%可信区间为(-2.239，-1.235)。这一结果表明，一个地区的成年人越多，其呼吸系统疾病的医院就诊率就越高。(吸烟的固定效应)的后验均值为0.748，相应的95%可信区间为(0.242，1.254)。这有力地证明了吸烟与呼吸系统疾病呈正相关。(医院密度)的后验均值为-110.804，相应的95%可信区间为(-121.582,-100.025)，这表明它对台湾呼吸系统疾病的医院就诊率有负向作用。(夏季固定效应)的后验均值为-0.02，(秋季固定效应)的后验均值为0.019，(冬季固定效应)的后验均值为0.050。因此，与春季相比，呼吸系统疾病夏季较轻，而冬季则严重得多，这与图1(a)中时间序列的"U"型吻合。
具有负的后验均值，表明通过控制、温度、相对湿度、儿童比例、吸烟率和季节，呼吸系统疾病的医院就诊率会随着时间的推移而略有下降。(空间效应)范围从-0.02到0.02，空间分布详见附录中的图C.2(b)。负空间效应表明，这些地区的粗略医院就诊率被高估了，负空间效应可以降低它。相反，正空间效应表明，这些地区的粗略医院就诊率被低估了，正空间效应可以提高它。
六、结论和讨论
本文提出了一个时空疾病制图模型，以研究台湾地区与呼吸系统疾病医院就诊率之间的关联。初步分析表明，对医院就诊率的影响在不同地区存在空间差异。因此，我们引入了空间变化系数来描述空间异质性。为了解决“change of support problem”，我们对点水平的、温度和相对湿度采用了SPDE方法。此外，为了考虑预测区域值的不确定性，我们使用SPDE进行了多次蒙特卡罗积分，然后对提出的疾病制图模型进行了迭代估计。鉴于台湾的实际应用是一个大规模的估计问题，即328个三级行政区加上54个每周时间点，我们采用INLA方法进行模型估计，以提高计算效率。通过各种模拟验证了有限样本的性能。在实际数据分析中，我们发现在台湾328个三级行政区域中，大部分地区的呼吸系统疾病就诊率与之间存在正相关关系，但其关系在不同地区和季节存在差异。
不过，我们的研究也有一些局限性，这可以启发我们今后进一步改进。我们只调查了对某一年呼吸系统疾病的影响。如果我们有更多的数据，就可以揭示多年的长期影响。此外，如果提供更多的环境和社会经济协变量，还可以进行更广泛的研究，还可以采用非线性和非参数模型来更好地理解PM2.5与呼吸系统疾病之间的关系。
参考文献
Bakar, K. S., Kokic, P. and Jin, H. (2016). Hierarchical spatially varying coefficient and temporal dynamic process models using spTDyn. Journal of Statistical Computation and Simulation 86, 820–840.
Banerjee, S., Carlin, B. P. and Gelfand, A. E. (2014). Hierarchical Modeling and Analysis for Spatial Data, 2nd edition. Boca Raton, FL: Chapman & Hall/CRC Press.
Blangiardo, M. and Cameletti, M. (2015). Spatial and Spatio-temporal Bayesian Models with R-INLA.Hoboken (HQ), NJ, USA: Wiley.
Cameletti, M., Gomez-Rubio, V. and Blangiardo, M. (2019). Bayesian modeling for spatially misaligned health and air pollution data through the INLA-SPDE approach. Spatial Statistics 31, 100353.
Carroll, R., Lawson, A. B., Faes, C., Kirby, R. S., Aregay, M. andWatjou, K. (2015). Comparing INLA and openbugs for hierarchical poisson modeling in disease mapping. Spatial and Spatio-Temporal Epidemiology 14, 45–54.
Carroll, R., Lawson, A. B., Kirby, R. S., Faes, C., Aregay,M. andWatjou, K. (2017). Space-time variation of respiratory cancers in south carolina: a flexible multivariate mixture modeling approach to risk estimation. Annals of Epidemiology 27, 42–51.
Chen, C. H., Wu, C. D., Chiang, H. C., Chu, D., Lee, K. Y., Lin, W. Y., Yeh, J. I., Tsai, K. W. and Guo, Y. L. (2019). The effects of fine and coarse particulate matter on lung function among the elderly. International Journal of Environmental Research and Public Health 9, 1–8.
Franco-Villoria, M., Ventrucci, M. and Rue, H. (2016). Hierarchical nearest-neighbor gaussian process models for large geostatistical datasets. Journal of the American Statistical Association 111, 800–812.
Gelfand, A. E., Diggle, P. J., Fuentes,M. and Guttorp, P. (2010). Handbook of Spatial Statistics. BocaRaton, FL: Chapman & Hall/CRC Press. Gelfand, A. E., Li, Z. and Carlin, B. P. (2001). On the change of support problem for spatio-temporal data. Biostatistics 2, 31–45.
Guo, C., Zhang, Z. L., Lau, A. K. H., Lin, C. Q. and Lao, X. Q. (2018). Effect of long-term exposure to fine particulate matter on lung function decline and risk of chronic obstructive pulmonary disease in Taiwan: a longitudinal, cohort study. Lancet Planet Health 2, 114–125.
Huang, H. C., Lin, F. C. F.,Wu,M. F., Nfor, O. N, Hsu, S. Y., Lung, C. C., Ho, C. C., Chen, C. Y. and Liaw, Y. P. (2019). Association between chronic obstructive pulmonary disease and PM2.5 in Taiwanese nonsmokers. International Journal of Hygiene and Environmental Health 222, 884–888.
Hwang, S. L., Guo, S. E., Chi,M. C., Chou, C. T., Lin, Y. C., Lin, C.M. and Chou, Y. L. (2016). Association between atmospheric fine particulate matter and hospital admissions for chronic obstructive pulmonary disease in southwestern Taiwan: a population-based study. International Journal of Environmental Research & Public Health 13, 366–374.
Jayes, L., Haslam, P. L., Gratziou, C. G., Powell, P., Britton, J., Vardavas, C., Jimenez-Ruiz, C. and Leonardi-Bee, J. (2016). SmokeHaz: systematic reviews and meta-analyses of the effects of smoking on respiratory health. Chest 150, 164–179.
Katanoda, K., Sobue, T., Satoh, H., Tajima, K., Suzuki, T.,Nakatsuka, H., Takezaki, T., Nakayama, T., Nitta, H., Tanabe, K. and others. (2011). An association between long-term exposure to ambient air pollution and mortality from lung cancer and respiratory diseases in Japan. Journal of Epidemiology, 1102090211.
Kuo, C. Y., Pan, R. H., Chan, C. K., Wu, C. Y., Van Phan, D. and Chan, C. L. (2019). Application of a time-stratified case-crossover design to explore the effects of air pollution and season on childhood asthma hospitalization in cities of differing urban patterns: big data analytics of government open data. International Journal of Environmental Research and Public Health 15, 647–659.
Kushta, J., Pozzer, A. and Lelieveld, J. (2018). Uncertainties in estimates of mortality attributable to ambient PM2.5 in Europe. Environmental Research Letters 13, 064029.
Lee, D., Rushworth, A. and Napier, G. (2018). Spatio-temporal areal unit modeling in R with conditional autoregressive priors using the CARBayesST package. Journal of Statistical Software 84, 1–39.
Lindgren, F., Rue, H. and Lindström, J. (2019). An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73, 423–498.
Oliveira, B. F. A., Ignotti, Eliane, Artaxo, Paulo, do Nascimento Saldiva, Paulo Hilário, Junger,Washington Leite and Hacon, Sandra. (2012). Risk assessment of PM2.5 to child residents in Brazilian amazon region with biofuel production. Environmental Health 11, 64.
Rue, H., Martino, S. and Chopin, N. (2009). Approximate Bayesian inference for latent gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society: Series B (statistical methodology) 71, 319–392.
Schrödle, B. and Held, L. (2011). A primer on disease mapping and ecological regression using INLA. Computational Statistics 26, 241–258.
Turner, M. C., Krewski, D., Iii, Cap, Chen, Y., Gapstur, S. M. and Thun, M. J. (2011). Long-term ambient fine particulate matter air pollution and lung cancer in a large cohort of never-smokers. American Journal of Respiratory & Critical Care Medicine 184, 1374–1381.
Wang, F., Wang, J., Gelfand, A. E. and Li, F. (2019). Disease mapping with generative models. The American Statistician 73, 213–223.
Wang, K. Y., Tang-Tat, C. and Gracie, S.M. (2013). An association between air pollution and daily outpatient visits for respiratory disease in a heavy industry area. PLoS One 8, e75220.
Wang, Y. C. and Lin, Y. K. (2015). Temperature effects on outpatient visits of respiratory diseases, asthma, and chronic airway obstruction in taiwan. International Journal of Biometeorology 59, 815–825.
Xing, Y. F., Xu, Y. H., Shi, M. H. and Lian, Y. X. (2016). The impact of PM2.5 on the human respiratory system. Journal of Thoracic Disease 8, E69.