时间病例交叉设计-急性病与短期暴露

写在前面：案例数据与SAS/R代码见文末。

一、背景

病例交叉研究（Case-Crossover Study），或叫病例交叉设计（Case-Crossover design），是由Maclure在1991年首次提出的，在环境流行病学领域大放异彩，当然临床流行病也偶见文献报道。它是一种用于研究短暂暴露对急性病的瞬间影响的流行病学方法。病例交叉研究是选择发生某种急性结局事件，分别调查事件发生时及发生前的暴露情况及程度，以判断暴露危险因素与结局事件有无关联及关联程度的一种研究方法。

病例交叉研究可以看成是自身前后对照的病例对照设计，也可以看成是回顾性队列研究。它以病例自身作为对照不仅避免了选择对照所引起的偏倚，而且可以避免各病例之间一些不可控制的因素( 如年龄、智力、遗传等) 所引起的偏倚，既省时省力，又节约样本，便于实施。但是病例交叉研究要求暴露与事件的发生时间间隔不能太久，且暴露不能有遗留效应。

1.1基本思想

比较相同研究对象在急性事件(如剧烈活动与心梗的发生；存在暴露累积效应滞留不适合; 慢性事件不适合)发生前一段时间的暴露情况与未发生事件的某时间段内的暴露情况是否一致。如果暴露与事件有关，那么在事件发生前的一段时间里的暴露量应该大一些。交叉在病例交叉研究里的体现就是，研究对象在不同的时间段里分别是病例和对照两种不同的状态。

1.2 分类

病例交叉研究可分为单向病例交叉研究和双向病例交叉研究。近年又拓展出时间分层病例交叉设计、时间空间分层病例交叉设计。

二、时间分层病例交叉研究

单向和双向策略引入了来自各种来源的偏差，例如暴露或健康事件的时间趋势(长期趋势)和季节性模式(季节趋势)以及对照日期选择并非独立。时间分层病例交叉设计是控制上述偏倚的最佳方法。

时间分层病例交叉研究的基本原理:将时间进行分层，病例期和对照期处于同一年、同一个月和同一个星期几，而且在同一时间层内，几个对照期是随机分布的，病例期并非固定在某一位置。例如，假设病例期发生在2014 年12 月12 日( 星期五) ，则2014 年12 月其他的星期五均被选为对照期。

三、理论

病例交叉资料分析方法采用配对检验(配对X2、配对t等)来进行。病例交叉研究中对照数据可以是个数(1：n配对条件logistic)，也可以是暴露人时(poisson回归)，故可选用条件logistics回归和条件泊松回归等方法进行分析，二者结果近似。

时间分层病例交叉设计相当于按时间层匹配的病例对照研究，时间分层病例交叉研究可以同时控制季节性与星期几效应等混杂因素、消除时间趋势偏倚，并能得到参数的无偏估计( 条件logistic 回归) 。而Poisson 回归和条件Poisson 回归通过设置哑变量stratum(年、月、星期几) 的方法，同样能达到按时间层匹配的目的。

四、范文

JACC论文PMID: 33478650。内容为短期暴露于环境空气污染与心梗（MI）死亡风险增加显著相关。研究调查了2013-18年间湖北省15w例MI死亡病例，通过结合居住地址及其周边空气质量监测站点数据，采用反距离权重插值法IDW计算各研究对象在病例日和对照日大气污染物暴露水平，拟合条件logistic回归模型定量评估环境空气污染诱发心梗死亡的暴露反应关系。

调整了气象因素（日均温度和露点温度）的暴露反应条件logistic回归模型中，可以观察到颗粒物短期暴露与心梗死亡呈非线性关联，在折点前，PM2.5（33.3 μg/m3）和PM10（57.3 μg/m3）浓度每增加10 μg/m3，心梗死亡风险分别增加4.14%（95% CI：1.25%–7.12%）和2.67%（0.80%–4.57%），而在较高浓度的折点后关联性逐渐减弱。双污染物模型结果提示，经调整气态污染物的影响后PM2.5和PM10与心梗死亡之间的关联仍保持稳定。不难发现，当嵌套模型中包含NO2时，颗粒物的健康损害效应更显著（异质性检验P<0.05）。

五、实战

5.1 sas宏 %Time_Stratified

%Time_Stratified( IN_DSN =, OUT_DSN =, DATE_VAR =, LAG_EFFECT =, LAG_VAR = ，COUNTS_VAR= )，中国医院统计杂志发表在2019年的%Time_Stratified宏程序，代码有部分问题，小编已做修改。%Time_Stratified适合时间分层病例交叉设计。只需指定输入数据集、输出数据集、日期变量、滞后效应、滞后效应相关变量、频数变量，就可自动完成对时间分层病例交叉设计资料的整理，以便进一步使用条件logistic 回归等方法分析。

另外%case_cross 宏程序，适合向前、向后、双向法，本文不赘述。

5.2 R代码实战

原始数据时间序列格式“londondataset2002_2006.dta”，是伦敦的一个空气污染与人群死亡数据集为例演示的时间分层病例交叉设计。数据集包含从2002 年1 月1 日到2006 年12 月31 日期间，伦敦市每天的臭氧浓度( ozone) 、温度( temperature) 、相对湿度( relative _ humidity ) 以及总死亡数( numdeaths)。所研究的暴露因素是臭氧浓度，健康结局是总死亡数，两个潜在的混杂因素是温度和相对湿度。

具体数据框与函数解析如下：

R代码如下：

setwd("D:/聂个人文件/我的公众号/时间分层病例交叉设计")
library(foreign)
data <- read.dta("londondataset2002_2006.dta")
#剔除缺失
# (MISSING EXCLUDED IN ESTIMATION BUT RE-INSERTED IN PREDICTION/RESIDUALS)
options(na.action="na.exclude")

#-------- 初步分析-------------#
# FIGURE 1
oldpar <- par(no.readonly=TRUE)
par(mex=0.8,mfrow=c(2,1))

#每日死亡子图SUB-PLOT FOR DAILY DEATHS, WITH VERTICAL LINES DEFINING YEARS
plot(data$date,data$numdeaths,pch=".",main="Daily deaths over time",
  ylab="Daily number of deaths",xlab="Date")
abline(v=data$date[grep("-01-01",data$date)],col=grey(0.6),lty=2)

# 每日臭氧子图
plot(data$date,data$ozone,pch=".",main="Ozone levels over time",
     ylab="Daily mean ozone level (ug/m3)",xlab="Date")
abline(v=data$date[grep("-01-01",data$date)],col=grey(0.6),lty=2)
par(oldpar)
layout(1)

summary(data)
#pacman::p_load("epicalc") 
des(data)
summ(data)
#相关系数
cor(data[,2:4])
#设置量纲
data$ozone10 <- data$ozone/10

#------------------------ 时间分层交叉设计，生成年月周 -------------------- 
data$temp <- data$temperature
data$rh <- data$relative_humidity
# 原始数据必须有月份年星期需要因子化
data$dow <- as.factor(weekdays(data$date))
data$month <- as.factor(months(data$date,abbr=TRUE))
data$year <- as.factor(substr(data$date,1,4))
#时间分层
data$stratum <- as.factor(data$year:data$month:data$dow)

#---------------方法1  条件logistic；3种方法结果一致---------
#条件logistic 回归得到的结果为OR,Poisson 回归和条件Poisson回归得到的是RR，结果一致！

#------  条件logistic需要个体数据形态，不可计算自相关、过离散 -------
library(season) 
model1<- casecross ( numdeaths~ozone + temp +rh,matchdow =T,stratamonth = T,data = data)
#同一年同一月同一个星期几
summary( model1)  

#---------------方法2 glm函数实现 poisson回归 主要必须 factor(stratum)加权---------
model2<- glm(numdeaths ~ozone + temp + rh +factor(stratum),data = data,family = poisson)
#输出stratum几百行
summary( model2)

#---------------方法3 gn函数实现 条件poisson回归---------
library(gnm) 
model3<-gnm(numdeaths ~ozone + temp + rh,data = data,family = poisson,
              eliminate = factor(stratum))
summary(model3)

#---------------方法4 quasipoisson代替poisson分布族控制过离散---------
model4 <- glm(numdeaths ~ozone + temp + rh +factor(stratum),data = data,family = quasipoisson)
summary( model4)
#---------------方法5 控制自相关---------
library(sModel)  
#提取模型残差项,加入残差项控制1阶自相关
reslag1 <- Lag (residuals(model4,type ="deviance"),1) 
model5 <- glm (numdeaths~ozone + temp + rh +reslag1,data = data,family = quasipoisson) 

#--------------- 方法6 单滞后模型------------
library(Epi)
library(tsModel) 
tablag1 <- matrix(NA,7+1,3,dimnames=list
                (paste("Lag",0:7),c("RR","ci.low","ci.hi")))
#建立7行3列的表格
for(i in 0:7){
  numdeathslag <- Lag(data$numdeaths,-i)
  model6 <- glm(numdeathslag~ozone+temp+rh+factor(stratum),data=data,family=quasipoisson)
  tablag1[i+1,] <- ci.lin(model6,subset="ozone",Exp=T)[5:7]
}
tablag1
##建立滞后变量，并依次估计不同滞后的参数

#--------------- 方法7 分布滞后模型 ---------------
library(dlnm)
cbozone<-crossbasis(data$ozone,lag=c(0,7),
                    argvar=list(type="lin",cen=0),arglag=list(type="integer"))
#建立臭氧与滞后的交叉基
model7<-glm(numdeaths~cbozone+temp+rh+factor(stratum),data=data,family=quasipoisson)
pred<-crosspred(cbozone,model7,at=10)
#cen=0指臭氧浓度参照水平为0ug/m3,at=10指相对于0μg/m3,臭氧浓度为10ug/m3的效应
tablag2<-with(pred,t(rbind(matRRfit,matRRlow,
                           matRRhigh)))
#提取RR值及可信区间
colnames(tablag2)<-c("RR","ci.low","ci.hi")
tablag2
pred$allRRfit
pred$allRRlow 
pred$allRRhigh
#不同滞后的累积效应

上述R code利用clogit条件logistic、glm带weight的一般poisson、gnm条件poisson 3种方法分析结果一致！时间分层病例交叉设计的条件logistic回归与泊松回归分析是近似等价的！其实就是罕见疾病的OR与RR类型的道理。发病率<< 5%的时候，OR近似RR。

六、发散思维

6.1 病例交叉优缺点：

6.2过离散、自相关、滞后等相关问题：

条件logistic回归是病例交叉设计最常用的统计分析方法，但其并不能控制时间序列数据的自相关、过离散等问题，对此我们可以利用Poisson 回归来控制这些问题，但Poisson回归也存在估计参数较多的缺点。

条件Poisson回归是Poisson 回归的进一步发展，其显著优点是减少了需要估计的参数，缩短了软件运算的时间，同时又不影响参数估计的准确性。

小结

时间分层病例交叉设计用于短暂暴露对急性发病关联研究。能很好的控制bias，包括长期趋势、季节趋势和人口学行为遗传因素。方法推荐条件logistic(OR值)、条件Poisson(RR值)，后者更可以适配过离散自相关问题。