Creating Database only with summer months

villes_s<-list()

for(i in 1:length(villes)) {
  villes_s[[i]]<-villes[[i]][villes[[i]]$saison=="summer",]
  villes_s[[i]]$temps <- ave(villes_s[[i]]$time,villes_s[[i]]$annee,  FUN = seq_along)
  villes_s[[i]]$time<-1:nrow(villes_s[[i]])
}

names(villes_s)<-cities

table of N obs per months to check to have only summer month

for (i in 1:length(villes)){
print(table(villes_s[[i]]$mois))
}
## 
##   1   2   3   4   5   6   7   8   9  10  11  12 
##   0   0   0   0   0 540 558 558 540   0   0   0 
## 
##   1   2   3   4   5   6   7   8   9  10  11  12 
##   0   0   0   0   0 540 558 558 540   0   0   0 
## 
##   1   2   3   4   5   6   7   8   9  10  11  12 
##   0   0   0   0   0 540 558 558 540   0   0   0 
## 
##   1   2   3   4   5   6   7   8   9  10  11  12 
##   0   0   0   0   0 540 558 558 540   0   0   0 
## 
##   1   2   3   4   5   6   7   8   9  10  11  12 
##   0   0   0   0   0 540 558 558 540   0   0   0 
## 
##   1   2   3   4   5   6   7   8   9  10  11  12 
##   0   0   0   0   0 540 558 558 540   0   0   0 
## 
##   1   2   3   4   5   6   7   8   9  10  11  12 
##   0   0   0   0   0 540 558 558 540   0   0   0 
## 
##   1   2   3   4   5   6   7   8   9  10  11  12 
##   0   0   0   0   0 540 558 558 540   0   0   0 
## 
##   1   2   3   4   5   6   7   8   9  10  11  12 
##   0   0   0   0   0 540 558 558 540   0   0   0 
## 
##   1   2   3   4   5   6   7   8   9  10  11  12 
##   0   0   0   0   0 540 558 558 540   0   0   0 
## 
##   1   2   3   4   5   6   7   8   9  10  11  12 
##   0   0   0   0   0 540 558 558 540   0   0   0 
## 
##   1   2   3   4   5   6   7   8   9  10  11  12 
##   0   0   0   0   0 540 558 558 540   0   0   0 
## 
##   1   2   3   4   5   6   7   8   9  10  11  12 
##   0   0   0   0   0 540 558 558 540   0   0   0 
## 
##   1   2   3   4   5   6   7   8   9  10  11  12 
##   0   0   0   0   0 540 558 558 540   0   0   0 
## 
##   1   2   3   4   5   6   7   8   9  10  11  12 
##   0   0   0   0   0 540 558 558 540   0   0   0 
## 
##   1   2   3   4   5   6   7   8   9  10  11  12 
##   0   0   0   0   0 540 558 558 540   0   0   0 
## 
##   1   2   3   4   5   6   7   8   9  10  11  12 
##   0   0   0   0   0 540 558 558 540   0   0   0 
## 
##   1   2   3   4   5   6   7   8   9  10  11  12 
##   0   0   0   0   0 540 558 558 540   0   0   0

Creating Heat Wave variable

trshld95<-c()

for (i in 1:length(villes_s)){
  trshld95[i]<-quantile(villes_s[[i]]$tempmax, probs=c(0.95),na.rm=TRUE)
  villes_s[[i]]$heat_wave<-0
  villes_s[[i]]$heat_wave[which(villes_s[[i]]$tempmax > trshld95[i])]<- 1 
}

 

Number of Heat Wave for each city:

Heat wave=0 Heat wave=1
bm 2087 109
bordeaux 2087 109
clermont 2088 108
dijon 2088 108
grenoble 2086 110
lehavre 2088 108
lille 2087 109
lyon 2087 109
marseille 2093 103
montpellier 2089 107
nancy 2086 110
nantes 2087 109
nice 2088 108
paris 2090 106
rennes 2089 107
rouen 2089 107
strasbourg 2086 110
toulouse 2086 110 
# creation of no2 variable of today and previous day
#function filter {stats}
for (i in 1:length(villes_s)){
villes_s[[i]]<-transform(villes_s[[i]], no2moy = as.vector(filter(no2,sides = 1, filter = rep(1, 2))/2), no2Lag1=Lag(no2,1),no2Lag2=Lag(no2,2))
}

Analysis

BOOTSTRAP METHOD

Code:

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:kableExtra':
## 
##     group_rows
## The following object is masked from 'package:nlme':
## 
##     collapse
## The following objects are masked from 'package:lubridate':
## 
##     intersect, setdiff, union
## The following objects are masked from 'package:Hmisc':
## 
##     src, summarize
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## create dataframe without missing values
for (i in 1:length(villes_s)){
  villes_s[[i]]<-villes_s[[i]][villes_s[[i]]$Dates>="2007-06-01",]
} 
     
  # start bootstrap
     
  nreps=1000 #number of bootstrap reps
  results<-matrix(NA,nrow = 18,ncol=19)
  colnames(results)<-c("city","CDE","CDE_2.5","CDE_97.5","RefInt","RefInt_2.5","RefInt_97.5","MedInt","MedInt_2.5","MedInt_97.5","PIE","PIE_2.5","PIE_97.5","OPM","OPM_2.5","OPM_97.5","OPAI","OPAI_2.5","OPAI_97.5")
  
for (i in (1:length(villes_s))){
  boot.cde<- rep(NA,nreps)
  boot.ref_int <- rep(NA,nreps)
  boot.med_int <- rep(NA,nreps)
  boot.pie<-rep(NA,nreps)
  boot.terr<-rep(NA,nreps)
  boot.total<-rep(NA,nreps)
  boot.prop_cde<-rep(NA,nreps)
  boot.prop_refint<-rep(NA,nreps)
  boot.prop_medint<-rep(NA,nreps)
  boot.prop_pie<-rep(NA,nreps)
  boot.opm<-rep(NA,nreps)
  boot.opai<-rep(NA,nreps)
  
  for (rep in 1:nreps)
      
    {
      dat_resample <- sample(1:nrow(villes_s[[i]]),nrow(villes_s[[i]]),replace=T)
      dat_boot <- villes_s[[i]][dat_resample,]
        
      # outcome model - Poisson Regression: Mortality, HW and O3
      m.outcome<-glm(nocc_tot~heat_wave + o3 + o3*heat_wave +no2+no2moy+ns(time,df=round(8*length(time)/365.25))+Jours+Vacances+hol,data=dat_boot,family=quasipoisson)
      
      # mediator model - multinomial linear regression
      m.mediator<-lm(o3~heat_wave+no2+no2Lag1+no2Lag2+no2moy+ns(time,df=round(8*length(time)/365.25))+Jours+Vacances+hol,data=dat_boot) 
      
      # save coefficients
      theta1<-coef(m.outcome)[2]
      
      theta2<-coef(m.outcome)[3]
     
      theta3<-coef(m.outcome)[43]
      
      beta1<-coef(m.mediator)[2]
      
      
      # Mean Values for covariates
      fJ<-data.frame(table(dat_boot$Jours))
      fJ<-fJ %>%  rename(Jours=Var1)
      
      new.data <- data.frame(
      
      heat_wave=0,
      no2=mean(dat_boot$no2,na.rm = TRUE),
      no2Lag1=mean(dat_boot$no2Lag1,na.rm = TRUE),
      no2Lag2=mean(dat_boot$no2Lag2,na.rm = TRUE),
      no2moy=mean(dat_boot$no2moy,na.rm=TRUE),
      time=median(dat_boot$time),
      Jours=c(fJ[which.max(fJ$Freq),][1]),
      Vacances=levels(dat_boot$Vacances)[1],
      hol=levels(dat_boot$hol)[1]
      
      )
      
      suppressWarnings(pred<-predict(m.mediator,newdata =new.data))
      
      # estimate CDE and NIE
      
      boot.cde[rep] <- exp(theta1) #NDE
      
      boot.ref_int[rep] <-exp(theta3*pred) # Reference Interaction 
      
      boot.med_int[rep] <-exp(theta3*beta1) # Mediated Interaction
      
      boot.pie[rep]<-exp(theta2*beta1) # Pure Indirect Effect
      
      boot.terr[rep] <- boot.cde[rep]+boot.ref_int[rep]+boot.med_int[rep]+boot.pie[rep] #
      
      boot.total[rep] <-boot.cde[rep]*boot.ref_int[rep]*boot.med_int[rep]*boot.pie[rep] # Total Effect Risk Ratio
                   
      boot.prop_cde[rep] <- boot.cde[rep]/boot.terr[rep] # Proportion CDE
  
      boot.prop_refint[rep] <- boot.ref_int[rep]/boot.terr[rep]   # Proportion Reference Interaction
      
      boot.prop_medint[rep] <- boot.med_int[rep]/boot.terr[rep]   # Proportion Mediated Interaction
      
      boot.prop_pie[rep] <- boot.pie[rep]/boot.terr[rep] # Proportion PIE
      
      boot.opm[rep]<- (boot.pie[rep]+boot.med_int[rep])/boot.terr[rep] #Overall Proportion Mediated
      
      boot.opai[rep]<- (boot.ref_int[rep]+boot.med_int[rep])/boot.terr[rep] # Overall Proportion Attributable to Interaction
        
          } #end bootstrap

  results[i,1]<-cities[[i]]
  results[i,2]<-median(boot.cde, na.rm=T)
  results[i,3]<-quantile(boot.cde, 0.025, na.rm=T)
  results[i,4]<-quantile(boot.cde, 0.975, na.rm=T)
  results[i,5]<-median(boot.ref_int, na.rm=T)
  results[i,6]<-quantile(boot.ref_int, 0.025, na.rm=T)
  results[i,7]<-quantile(boot.ref_int, 0.975, na.rm=T)
  results[i,8]<-median(boot.med_int, na.rm=T)
  results[i,9]<-quantile(boot.med_int, 0.025, na.rm=T)
  results[i,10]<-quantile(boot.med_int, 0.975, na.rm=T)
  results[i,11]<-median(boot.pie, na.rm=T)
  results[i,12]<-quantile(boot.pie, 0.025, na.rm=T)
  results[i,13]<-quantile(boot.pie, 0.975, na.rm=T)
  results[i,14]<-median(boot.opm, na.rm=T)
  results[i,15]<-quantile(boot.opm, 0.025, na.rm=T)
  results[i,16]<-quantile(boot.opm, 0.975, na.rm=T)
  results[i,17]<-median(boot.opai, na.rm=T)
  results[i,18]<-quantile(boot.opai, 0.025, na.rm=T)
  results[i,19]<-quantile(boot.opai, 0.975, na.rm=T)
  
     }
  
  # outputs and confidence intervals

Results:

res<-data.frame(results)

for (i in (2:19)){
res[,i]<-as.numeric(as.character(res[,i]))
}

res[,2:19]<-round(res[,2:19],digits = 4)

kable(res)%>%kable_styling("condensed")%>%
  column_spec(1, bold = T, border_right = T)%>%
  column_spec(4, border_right = T)%>%
  column_spec(7, border_right = T)%>%
  column_spec(10, border_right = T)%>%
  column_spec(13, border_right = T)%>%
  column_spec(16, border_right = T)
city CDE CDE_2.5 CDE_97.5 RefInt RefInt_2.5 RefInt_97.5 MedInt MedInt_2.5 MedInt_97.5 PIE PIE_2.5 PIE_97.5 OPM OPM_2.5 OPM_97.5 OPAI OPAI_2.5 OPAI_97.5
bm 1.2990 0.5823 2.9833 0.8229 0.4981 1.3168 0.8988 0.6664 1.1619 1.0302 0.9771 1.0886 0.4758 0.3270 0.5385 0.4248 0.2244 0.6068
bordeaux 0.9725 0.4440 2.1806 1.1354 0.6408 1.9819 1.0505 0.8398 1.3251 1.0176 0.9822 1.0581 0.4933 0.3973 0.5140 0.5245 0.3188 0.6928
clermont 1.1805 0.4341 3.5141 0.9681 0.4055 2.0120 0.9909 0.7797 1.2337 1.0142 0.9773 1.0535 0.4767 0.3145 0.5063 0.4708 0.2110 0.6912
dijon 0.8324 0.2367 2.0094 1.1487 0.6386 2.5067 1.0571 0.8404 1.4821 1.0284 0.9756 1.0810 0.5054 0.4109 0.5328 0.5432 0.3295 0.7615
grenoble 1.3141 0.6713 2.6491 0.8533 0.5034 1.3844 0.9474 0.7911 1.1205 1.0322 0.9977 1.0698 0.4778 0.3675 0.5133 0.4342 0.2638 0.5923
lehavre 0.8154 0.3389 1.7221 1.1400 0.7076 1.9075 1.0596 0.8577 1.3657 1.0653 1.0045 1.1299 0.5177 0.4386 0.5446 0.5394 0.3578 0.6957
lille 0.9590 0.6965 1.3273 1.0216 0.8372 1.2675 1.0105 0.9167 1.1253 1.0456 1.0158 1.0790 0.5088 0.4753 0.5279 0.5030 0.4280 0.5790
lyon 1.3537 0.7336 2.3671 0.8608 0.5955 1.2935 0.9494 0.8238 1.0949 1.0147 0.9923 1.0378 0.4691 0.3828 0.5113 0.4329 0.2924 0.5757
marseille 1.0527 0.5736 1.8412 1.0119 0.6330 1.6272 1.0014 0.9526 1.0628 1.0055 0.9969 1.0158 0.4907 0.4418 0.5047 0.4944 0.3551 0.6285
montpellier 1.0376 0.4281 2.0259 0.8836 0.4343 1.9959 1.0000 0.9834 1.0220 0.9998 0.9940 1.0033 0.4999 0.4231 0.5218 0.4812 0.3247 0.6798
nancy 1.2194 0.3510 2.4876 0.8711 0.5298 1.9305 0.9522 0.7992 1.2816 0.9907 0.9549 1.0287 0.4794 0.3751 0.5242 0.4509 0.2781 0.7081
nantes 0.5666 0.2372 1.2922 1.4637 0.8830 2.4452 1.1899 0.9467 1.5055 1.0316 0.9915 1.0755 0.5147 0.4668 0.5336 0.6234 0.4454 0.7579
nice 0.7796 0.5570 1.0634 1.2849 0.9722 1.8150 1.0052 0.9913 1.0279 0.9997 0.9933 1.0049 0.4918 0.4560 0.5064 0.5639 0.4894 0.6423
paris 1.1930 0.9677 1.4468 0.9267 0.8188 1.0590 0.9694 0.9190 1.0214 1.0319 1.0195 1.0453 0.4855 0.4623 0.5033 0.4599 0.4135 0.5091
rennes 1.9142 0.5456 5.7366 0.6819 0.3224 1.6027 0.8353 0.5898 1.2149 1.0444 0.9852 1.1129 0.4204 0.2123 0.5203 0.3387 0.1196 0.6435
rouen 1.1893 0.7280 1.8146 0.9454 0.7290 1.2532 0.9654 0.8228 1.1523 1.0559 1.0077 1.1054 0.4871 0.4293 0.5275 0.4597 0.3510 0.5754
strasbourg 0.7983 0.4656 1.3296 1.1959 0.8930 1.6565 1.0906 0.9459 1.2778 1.0262 0.9896 1.0668 0.5131 0.4721 0.5317 0.5567 0.4372 0.6631
toulouse 0.8809 0.3680 1.9945 1.1396 0.5743 2.3193 1.0339 0.8687 1.2504 1.0307 1.0013 1.0595 0.4972 0.4216 0.5154 0.5330 0.3225 0.7170

Forestplot

Controlled Direct Effect

Pure Indirect Effect

Reference Interaction

Reference Interaction