Time series analysis: Case-crossover analysis
2024-07-22
Case-crossover analysis
In the case-crossover design, the study population consists of subjects who have experienced an episode of the health outcome of interest. Similar to a crossover study, each subject serves as his or her own control. As in a matched case-control study, the inference is based on a comparison of exposure distribution rather than the risk of disease. The case-crossover study is most suitable for studying relations with the following characteristics: 1) the individual exposure varies within short time intervals; 2) the disease has abrupt onset and short latency for detection; and 3) the induction period is short.
Case-crossover design allows use of routinely monitored air pollution information and at the same time makes it possible to study individuals rather than days as the unit of observation. .
##read data
data <- read.csv("/Users/kennyorielolana/Documents/DATA ANALYSIS and Mapping/Training in Time series analysis /data_urolithiasis.csv")
##Reading R packages
library(dlnm)## This is dlnm 2.4.7. For details: help(dlnm) and vignette('dlnmOverview').library(splines)
library(survival)
head(data)## date studyid sex age event time dow l0_wbgt l1_wbgt l2_wbgt
## 1 2010-01-01 382405_1 2 66 0 2 6 13.292202 12.711058 12.414207
## 2 2010-01-01 86767 2 53 0 2 6 7.000611 5.941786 8.097641
## 3 2010-01-01 307636_1 1 41 0 2 6 13.292202 12.711058 12.414207
## 4 2010-01-01 6127201 1 28 0 2 6 7.124393 20.376735 18.107819
## 5 2010-01-01 5859 1 75 0 2 6 2.210671 18.645376 18.426286
## 6 2010-01-01 6088901 2 18 0 2 6 7.124393 20.376735 18.107819
## l3_wbgt l4_wbgt l5_wbgt l6_wbgt l7_wbgt holiday
## 1 11.202209 8.232980 9.561680 7.606777 9.417363 1
## 2 8.731275 8.817321 8.068007 6.320272 4.392536 1
## 3 11.202209 8.232980 9.561680 7.606777 9.417363 1
## 4 17.110460 17.114432 15.423693 12.200597 14.839795 1
## 5 17.955405 14.915938 13.821673 11.252827 12.311683 1
## 6 17.110460 17.114432 15.423693 12.200597 14.839795 1## date~date, sex~1 (male) and 2 (female), event~1 refers to the case period and 0 refers to the control period; wbgt~ Wet Bulb Globe Temperature
summary(data)## date studyid sex age
## Length:185777 Length:185777 Min. :1.000 Min. :-2.00
## Class :character Class :character 1st Qu.:1.000 1st Qu.:36.00
## Mode :character Mode :character Median :1.000 Median :46.00
## Mean :1.332 Mean :47.02
## 3rd Qu.:2.000 3rd Qu.:58.00
## Max. :2.000 Max. :98.00
##
## event time dow l0_wbgt
## Min. :0.0000 Min. :1.000 Min. :1.000 Min. :-15.49
## 1st Qu.:0.0000 1st Qu.:2.000 1st Qu.:2.000 1st Qu.: 13.69
## Median :0.0000 Median :2.000 Median :4.000 Median : 20.40
## Mean :0.2272 Mean :1.773 Mean :3.931 Mean : 18.68
## 3rd Qu.:0.0000 3rd Qu.:2.000 3rd Qu.:6.000 3rd Qu.: 24.99
## Max. :1.0000 Max. :2.000 Max. :7.000 Max. : 30.20
## NA's :25
## l1_wbgt l2_wbgt l3_wbgt l4_wbgt
## Min. :-15.49 Min. :-15.49 Min. :-15.49 Min. :-15.49
## 1st Qu.: 13.73 1st Qu.: 13.73 1st Qu.: 13.72 1st Qu.: 13.75
## Median : 20.42 Median : 20.42 Median : 20.43 Median : 20.42
## Mean : 18.70 Mean : 18.71 Mean : 18.71 Mean : 18.72
## 3rd Qu.: 25.01 3rd Qu.: 24.99 3rd Qu.: 25.01 3rd Qu.: 25.04
## Max. : 30.20 Max. : 30.20 Max. : 30.20 Max. : 30.20
## NA's :15 NA's :19 NA's :28 NA's :15
## l5_wbgt l6_wbgt l7_wbgt holiday
## Min. :-15.49 Min. :-14.38 Min. :-15.49 Min. :0.00000
## 1st Qu.: 13.75 1st Qu.: 13.75 1st Qu.: 13.76 1st Qu.:0.00000
## Median : 20.39 Median : 20.40 Median : 20.43 Median :0.00000
## Mean : 18.72 Mean : 18.72 Mean : 18.72 Mean :0.07141
## 3rd Qu.: 25.06 3rd Qu.: 25.02 3rd Qu.: 25.01 3rd Qu.:0.00000
## Max. : 30.20 Max. : 30.20 Max. : 30.20 Max. :1.00000
## NA's :28 NA's :29 NA's :23###DLNM function
## 7-day lag with 2 degrees of freedom
lag <- 7
lagnk <- 2
fmla<-as.formula(event ~ cb+strata(studyid)+as.factor(holiday))
## Removing Temperature Extremes
wbgt <- data[,8:(8+7)]
p99<-quantile(wbgt, probs = 0.999,na.rm = T)
p1<-quantile(wbgt, probs = 0.001,na.rm = T)
wbgt[wbgt > p99] <- NA
wbgt[wbgt < p1] <- NA
bound<-c(round(p1,1),round(p99,1))
##
t1<-round(quantile(wbgt,0.01,na.rm = T),1)
t2<-round(quantile(wbgt,0.025,na.rm = T),1)
t3<-round(quantile(wbgt,0.975,na.rm = T),1)
t4<-round(quantile(wbgt,0.99,na.rm = T),1)
## building up cross-bases
argvar <- list(fun="ns", knots = equalknots(wbgt,fun = "ns",df=4),
Bound=bound)
cb <- crossbasis(wbgt,lag=lag,argvar=argvar,
arglag=list(knots=logknots(lag,lagnk)))
model<-clogit(fmla,data=data,na.action=na.exclude)
pred <- crosspred(cb,model,by=0.1,cumul = TRUE)## centering value unspecified. Automatically set to 10cen<-pred$predvar[which.min(pred$allfit)]
print(cen) ## MMT## [1] 15pred_S <- crosspred(cb,model,by=0.1,cumul = TRUE,cen=cen,
at=c(seq(bound[1],bound[2],by=0.1),t3,t4))
##Calculating OR
x1<-t(as.matrix(round(with(pred,cbind(allRRfit,allRRlow,allRRhigh)[as.character(t1),]),3)))
x2<-t(as.matrix(round(with(pred,cbind(allRRfit,allRRlow,allRRhigh)[as.character(t2),]),3)))
x3<-t(as.matrix(round(with(pred,cbind(allRRfit,allRRlow,allRRhigh)[as.character(t3),]),3)))
x4<-t(as.matrix(round(with(pred,cbind(allRRfit,allRRlow,allRRhigh)[as.character(t4),]),3)))
result <- data.frame(rbind(x1,x2,x3,x4))
rownames(result) <- c(t1,t2,t3,t4)
##Exposure-response curve
# Define the x-axis label with the degree symbol
xlab <- expression(paste("WBGT (", degree, "C)"))
# Set up plotting parameters
par(mar = c(4, 4, 3, 1), cex.axis = 0.9, mgp = c(2.5, 1, 0), cex.lab = 1)
# Ensure that pred_S is properly defined and that the plot function supports your arguments
if (exists("pred_S")) {
plot(pred_S, "overall", col = 2, cex.axis = 1, xlab = xlab, ylab = "Odds ratio", ylim = c(0.8, 2.0),
main = paste0("Exposure-response curve"), cex.main = 1, ci = "area", ci.arg = list(col = grey(0.85)))
lines(pred_S, col = 2, lwd = 2)
} else {
warning("The object 'pred_S' does not exist.")
}
# Save the plot as a TIFF file
tiff("figure_ER.tiff", res = 300, height = 1200, width = 1200)
# Reproduce the plot for saving
par(mar = c(4, 4, 3, 1), cex.axis = 0.9, mgp = c(2.5, 1, 0), cex.lab = 1)
plot(pred_S, "overall", col = 2, cex.axis = 1, xlab = xlab, ylab = "Odds ratio", ylim = c(0.8, 2.0),
main = paste0("Exposure-response curve"), cex.main = 1, ci = "area", ci.arg = list(col = grey(0.85)))
lines(pred_S, col = 2, lwd = 2)
# Close the TIFF device
dev.off()## quartz_off_screen
## 2## lag pattern curve
# Define the x-axis label
xlab <- "Lag (day)"
ylab <- "Odds ratio"
# Set up plotting parameters
par(mar = c(4, 4, 3, 1), cex.axis = 0.9, mgp = c(2.5, 1, 0), cex.lab = 1)
# Ensure that pred_S and t4 are properly defined
if (exists("pred_S") && exists("t4")) {
plot(pred_S, "slices", var = t4,
xlim = c(0, 7), cex.main = 1,
col = 2, lwd = 2,
xlab = xlab,
ylab = ylab, cex.lab = 1,
ci = "area", ci.arg = list(col = grey(0.85)),
main = paste0("Lag pattern"))
lines(pred_S, col = 1, lwd = 2)
} else {
warning("The objects 'pred_S' or 't4' do not exist.")
}
# Save the plot as a TIFF file
tiff("figure_lag.tiff", res = 300, height = 1200, width = 1200)
# Set up plotting parameters
par(mar = c(4, 4, 3, 1), cex.axis = 0.9, mgp = c(2.5, 1, 0), cex.lab = 1)
# Ensure that pred_S and t4 are properly defined
if (exists("pred_S") && exists("t4")) {
plot(pred_S, "slices", var = t4,
xlim = c(0, 7), cex.main = 1,
col = 2, lwd = 2,
xlab = xlab,
ylab = ylab, cex.lab = 1,
ci = "area", ci.arg = list(col = grey(0.85)),
main = paste0("Lag pattern"))
lines(pred_S, col = 1, lwd = 2)
} else {
warning("The objects 'pred_S' or 't4' do not exist.")
}
# Close the TIFF device
dev.off()## quartz_off_screen
## 2