########################################################################################################################### ########################################################################################################################### ##### Workshop "Modelling desert dust exposure events for epidemiological short-term health effects studies". ##### ##### Massimo Stafoggia (m.stafoggia@deplazio.it) and Aurelio Tobías (aurelio.tobias@idaea.csic.es) ##### ##### ##### ##### 31st Annual Conference of the International Society of Environmental Epidemiology. ##### ##### Utrecht, The Netherlands, 24th August 2019. ##### ##### ##### ##### LAST UPDATE: 17/08/2019 ##### ########################################################################################################################### ########################################################################################################################### ### Remove objects rm(list=ls()) ### Call the relevant libraries library(splines) ##### Define the folder path #setwd("C:\\Max\\DUST\\ISEE2019\\") ##### Import the stata dataset with simulated data # Call the data and check names of the variables and dimension of the dataset data <- readRDS("data_isee2019dust.rds") names(data) dim(data) # Set the dataframe as a data.table setDT(data) # Sort it by date setkey(data, date) # Describe the data str(data) # date: current date (from 1/1/2005 to 31/12/2010) # trend: progressive number from 1 to 2191 # yy: year # mm: month # dd: day of the month # dow: day of the week # allnat: daily mortality counts for non-accidental causes (simulated data) # temp: daily mean air temperature # dust: 0/1 variable for dust advection days # pm10: daily PM10 concentrations # pm10natural: daily PM10 concentrations from natural sources (>0 only when DUST=1) # pm10local: daily PM10 concentrations from non-natural sources (equal to PM10-PM10natural) ### Create lag terms (up to 6 days, and cumulative lags) for cold and hot temperatures adjustment n <- dim(data)[1] data$l1temp <- c(NA,data$temp[1:(n-1)]) data$l2temp <- c(NA,NA,data$temp[1:(n-2)]) data$l3temp <- c(NA,NA,NA,data$temp[1:(n-3)]) data$l4temp <- c(NA,NA,NA,NA,data$temp[1:(n-4)]) data$l5temp <- c(NA,NA,NA,NA,NA,data$temp[1:(n-5)]) data$l6temp <- c(NA,NA,NA,NA,NA,NA,data$temp[1:(n-6)]) data$l16temp <- (data$l1temp + data$l2temp + data$l3temp + data$l4temp + data$l5temp + data$l6temp)/6 data$l16tempcold <- ((data$l16temp - median(data$l16temp,na.rm=T))*(data$l16temp <= median(data$l16temp,na.rm=T)))+median(data$l16temp,na.rm=T) data$l01temp <- (data$temp + data$l1temp)/2 data$l01temphot <- ((data$l01temp - median(data$l01temp,na.rm=T))*(data$l01temp >= median(data$l01temp,na.rm=T)))+median(data$l01temp,na.rm=T) names(data) dim(data) # 2191 x 22 ### Adjustment for temperature using separate terms for cold temperatures and hot temperatures: # 1. Cold temperatures: defined on days with lag 1-6 temp < median, spline with one inner knot at 25th percentile # 2. Hot temperatures: defined on days with lag 0-1 temp > median, spline with two inner knots at 75th and 90th percentile # See MED-PARTICLES papers for reference (Stafoggia et al. 2013; Stafoggia et al. 2016) perc.hot <- c(75/100, 90/100) perc.cold <- c(25/100) ktemphot.lag01 <- as.vector(quantile(data$l01temphot, perc.hot, na.rm=T)) ktempcold.lag16 <- as.vector(quantile(data$l16tempcold,perc.cold,na.rm=T)) ktemphot.lag01 # 21.83359 25.80250 ktempcold.lag16 # 9.988542 ###################################### ### Definition of the adjustment model ###################################### mod <- glm(allnat ~ ns(trend, df=4*6)+as.factor(dow)+ns(l01temphot,k=ktemphot.lag01)+ns(l16tempcold,k=ktempcold.lag16), data=data, family=quasipoisson, na=na.omit) ############################################ ### 1. Effect of DUST as binary (0/1) metric ############################################ ##### RESEARCH QUESTION: Is mortality higher on DUST days compared with NO-DUST days? mod.1 <- update(mod, .~. + dust) summary(mod.1) # Estimate Std. Error t value Pr(>|t|) # dust 0.03411 0.01094 3.118 0.001847 ** ### INTERPRETATION: # On days WITH DUST EVENTS mortality increases by 3.5% (95%CI = 1.3%; 5.7%) ### CALCULATIONS: # (exp(0.03411)-1)*100 (exp(0.03411-1.96*0.01094)-1)*100 (exp(0.03411+1.96*0.01094)-1)*100 ############################################################################################################# ##### RESEARCH QUESTION: Is there an association of DUST events with mortality independent on PM10 increases? mod.1b <- update(mod, .~. + dust + pm10) summary(mod.1b) # Estimate Std. Error t value Pr(>|t|) # dust 0.0311292 0.0111738 2.786 0.005387 ** ### INTERPRETATION: # On days WITH DUST EVENTS mortality increases by 3.1% (95%CI = 0.9%; 5.5%), regardless of PM10 concentrations ### CALCULATIONS: # (exp(0.0311292)-1)*100 (exp(0.0311292-1.96*0.0111738)-1)*100 (exp(0.0311292+1.96*0.0111738)-1)*100 ############################################################################################################# ################################################# ### 2. Effect of PM10 by DUST binary (0/1) metric ################################################# ### As the first step, we want to check whether PM10 is associated to mortality. ##### RESEARCH QUESTION: Is PM10 associated with natural mortality? mod.2a <- update(mod, .~. + pm10) summary(mod.2a) # Estimate Std. Error t value Pr(>|t|) # pm10 0.0005337 0.0002996 1.781 0.075059 . ### INTERPRETATION: # Mortality increases by 0.5% (95%CI = -0.1%; 1.1%) per each 10 ug/m3 increment in total PM10 ### CALCULATIONS: # (exp(0.0005337*10)-1)*100 (exp((0.0005337-1.96*0.0002996)*10)-1)*100 (exp((0.0005337+1.96*0.0002996)*10)-1)*100 ############################################################################################################# ### As the second step, we check whether PM10-mortality association is confounded by DUST advection ##### RESEARCH QUESTION: Is there an association between PM10 and mortality independent on DUST events? mod.2b <- update(mod, .~. + dust + pm10) summary(mod.2b) # Estimate Std. Error t value Pr(>|t|) # pm10 0.0003900 0.0003071 1.270 0.204242 ### INTERPRETATION: # Mortality increases by 0.4% (95%CI = -0.21; 1.0%) per each 10 ug/m3 increment in total PM10, independently on DUST advections ### CALCULATIONS: # (exp(0.0003900*10)-1)*100 (exp((0.0003900-1.96*0.0003071)*10)-1)*100 (exp((0.0003900+1.96*0.0003071)*10)-1)*100 ############################################################################################################# ### Now we consider DUST as an effect modifer of the PM10-mortality association. ##### RESEARCH QUESTION: Is the association between PM10 and mortality different between DUST and NO-DUST days? mod.3 <- update(mod, .~. + dust*pm10) summary(mod.3) # Estimate Std. Error t value Pr(>|t|) # pm10 0.0003179 0.0003329 0.955 0.339749 # dust:pm10 0.0004336 0.0007705 0.563 0.573687 ### INTERPRETATION: # 1. On days with DUST = 0, per each 10 ug/m3 increase in PM10 the following day mortality increases by 0.3% (95%CI = -0.3%; 1.00%) # 2. On days WITH DUST = 1, per each 10 ug/m3 increase in PM10 the following day mortality increases by 0.8% (95%CI = -0.6%; 2.2%) ### CALCULATIONS: # 1. (exp(0.0003179*10)-1)*100 (exp((0.0003179-1.96*0.0003329)*10)-1)*100 (exp((0.0003179+1.96*0.0003329)*10)-1)*100 # 2. : # n <- length(mod.3$coefficients) # b.sum <- mod.3$coefficients[n-1] + mod.3$coefficients[n] # se.sum <- sqrt(sum(summary(mod.3)$cov.scaled[c(n-1,n),c(n-1,n)])) # z <- mod.3$coefficients[n]/(sqrt(summary(mod.3)$cov.scaled[n,n])) # p <- 2*(1-pnorm(abs(z))) # perc.inc.10 <- (exp(b.sum*10)-1)*100 # low.10 <- (exp((b.sum-1.96*se.sum)*10)-1)*100 # up.10 <- (exp((b.sum+1.96*se.sum)*10)-1)*100 # results <- c(b.sum, se.sum, p, perc.inc.10, low.10, up.10) # names(results) <- c("beta","se","p.interaction","%IR","95% LOW", "95% HIGH") # results # beta se p.interaction %IR 95% LOW 95% HIGH # 0.0007514410 0.0007113507 0.5736246221 0.7542713765 -0.6407447751 2.1688737263 ############################################################################################################# ##################################### ### 3. Effect of source-specific PM10 ##################################### ### Now we consider natural and local PM10 as two different exposures ##### RESEARCH QUESTION: Are natural (dust) and local (anthropogenic) sources of PM10 independently associated with mortality? mod.4 <- update(mod,.~. + pm10natural + pm10local) summary(mod.4) # Estimate Std. Error t value Pr(>|t|) # pm10natural 0.0013220 0.0006727 1.965 0.049526 * # pm10local 0.0004338 0.0003180 1.364 0.172709 ### INTERPRETATION: # 1. Mortality increases by 1.3% (95%CI = 0.0%; 2.7%) per each 10 ug/m3 increase in natural PM10 # 2. Mortality increases by 0.4% (95%CI = -0.2%; 1.1%) per each 10 ug/m3 increase in local PM10 ### CALCULATIONS: # 1. (exp(0.0013220*10)-1)*100 (exp((0.0013220-1.96*0.0006727)*10)-1)*100 (exp((0.0013220+1.96*0.0006727)*10)-1)*100 # 2. (exp(0.0004338*10)-1)*100 (exp((0.0004338-1.96*0.0003180)*10)-1)*100 (exp((0.0004338+1.96*0.0003180)*10)-1)*100 ######################################################################################################### ############################################# ### 4. Effect of source-specific PM10 by DUST ############################################# ##### RESEARCH QUESTIONS: ##### 1. is the association between local (non-desert) sources of PM10 with mortality different on DUST versus NO-DUST days? ##### 2. are these associations independent from natural (desert) sources? mod.5 <- update(mod,.~. + dust + pm10natural + pm10local + dust:pm10local) summary(mod.5) # Estimate Std. Error t value Pr(>|t|) # pm10natural -5.459e-05 9.299e-04 -0.059 0.953197 # pm10local 3.189e-04 3.329e-04 0.958 0.338102 # dust:pm10local 1.304e-03 9.958e-04 1.309 0.190682 ### INTERPRETATION: # 1. On days WITHOUT DUST EVENTS, mortality increases by 0.3% (95%CI = -0.3%; 1.0%) per 10 ug/m3 increases in LOCAL PM10 # On days WITH DUST EVENTS, mortality increases by 1.6% (95%CI = -0.24%; 3.55%) per 10 ug/m3 increases in LOCAL PM10 # 1. Mortality decreases by -0.1% (95%CI = -1.9%; 1.8%) per 10 ug/m3 increases in NATURAL PM10 ### CALCULATIONS: # 1. (exp(-5.459e-05*10)-1)*100 (exp((-5.459e-05-1.96*9.299e-04)*10)-1)*100 (exp((-5.459e-05+1.96*9.299e-04)*10)-1)*100 # 2. (exp(3.189e-04*10)-1)*100 (exp((3.189e-04-1.96*3.329e-04)*10)-1)*100 (exp((3.189e-04+1.96*3.329e-04)*10)-1)*100 # 3. : # n <- length(mod.5$coefficients) # b.sum <- mod.5$coefficients[n-1] + mod.5$coefficients[n] # se.sum <- sqrt(sum(summary(mod.5)$cov.scaled[c(n-1,n),c(n-1,n)])) # z <- mod.5$coefficients[n]/(sqrt(summary(mod.5)$cov.scaled[n,n])) # p <- 2*(1-pnorm(abs(z))) # perc.inc.10 <- (exp(b.sum*10)-1)*100 # low.10 <- (exp((b.sum-1.96*se.sum)*10)-1)*100 # up.10 <- (exp((b.sum+1.96*se.sum)*10)-1)*100 # results <- c(b.sum, se.sum, p, perc.inc.10, low.10, up.10) # names(results) <- c("beta","se","p.interaction","%IR","95% LOW", "95% HIGH") # results # beta se p.interaction %IR 95% LOW 95% HIGH # 0.0016224929 0.0009514915 0.1905337124 1.6357267902 -0.2421368714 3.5489397649 ########################################################################################################################### ########################################################################################################################### ##### ##### ##### END OF EXAMPLE ##### ##### ##### ########################################################################################################################### ###########################################################################################################################