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),nocc_tot=scale(nocc_tot),o3=scale(o3))
}

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.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<-lm(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)
      
      # 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] <- theta1 #NDE
      
      boot.ref_int[rep] <-theta3*pred # Reference Interaction 
      
      boot.med_int[rep] <-theta3*beta1 # Mediated Interaction
      
      boot.pie[rep]<-theta2*beta1 # Pure Indirect Effect
      
      boot.total[rep] <- boot.cde[rep]+boot.ref_int[rep]+boot.med_int[rep]+boot.pie[rep] # Total Effect 
                   
      boot.prop_cde[rep] <- boot.cde[rep]/boot.total[rep] # Proportion CDE
  
      boot.prop_refint[rep] <- boot.ref_int[rep]/boot.total[rep]   # Proportion Reference Interaction
      
      boot.prop_medint[rep] <- boot.med_int[rep]/boot.total[rep]   # Proportion Mediated Interaction
      
      boot.prop_pie[rep] <- boot.pie[rep]/boot.total[rep] # Proportion PIE
      
      boot.opm[rep]<- (boot.pie[rep]+boot.med_int[rep])/boot.total[rep] #Overall Proportion Mediated
      
      boot.opai[rep]<- (boot.ref_int[rep]+boot.med_int[rep])/boot.total[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 0.1412 -0.7018 1.0240 0.0194 -0.0864 0.2705 -0.2763 -1.0538 0.3618 0.0808 -0.0628 0.2359 -0.6686 -26.1905 20.1397 -0.3755 -31.1243 22.5345
bordeaux 0.2607 -0.6217 1.1210 0.0016 -0.1187 0.1583 0.1919 -0.5918 1.0442 0.0552 -0.0571 0.1770 0.4838 -1.0734 2.8524 0.3821 -1.2629 2.4590
clermont 0.3400 -0.3752 1.0915 0.0003 -0.1546 0.1554 -0.0155 -0.6531 0.5672 0.0326 -0.0531 0.1204 0.0303 -4.7061 4.9687 -0.0521 -3.9143 3.6069
dijon -0.1415 -1.0810 0.6002 -0.0029 -0.1694 0.1111 0.1319 -0.3966 0.8957 0.0599 -0.0531 0.1669 -0.2436 -27.4860 16.5546 -0.3467 -21.9906 14.8206
grenoble 0.2760 -0.2379 0.8767 -0.0002 -0.1238 0.1255 -0.1214 -0.6497 0.3701 0.0825 -0.0083 0.1774 -0.1190 -6.3430 5.5803 -0.3634 -8.1836 5.4100
lehavre -0.1871 -0.9068 0.4783 -0.0256 -0.1839 0.0874 0.1614 -0.3608 0.7952 0.1444 0.0085 0.2823 0.3065 -25.8867 21.4164 -0.1526 -16.0917 11.9523
lille -0.0813 -0.6699 0.5904 0.0021 -0.0839 0.1182 0.0573 -0.3838 0.5435 0.1907 0.0677 0.3282 0.9583 -15.2135 13.4889 0.2201 -9.5463 8.4310
lyon 0.5562 -0.2036 1.2175 0.0133 -0.1219 0.2544 -0.1875 -0.7212 0.3540 0.0499 -0.0266 0.1270 -0.3443 -1.8354 2.7581 -0.3352 -2.1356 1.8961
marseille 0.2318 -0.3609 0.7151 -0.0005 -0.1854 0.1537 0.0077 -0.2160 0.2770 0.0231 -0.0127 0.0661 0.0242 -4.5597 3.9724 -0.0016 -3.2108 3.6511
montpellier -0.2280 -0.5140 0.0716 -0.0111 -0.2145 0.1217 -0.0003 -0.0406 0.0497 -0.0005 -0.0148 0.0085 0.0055 -0.3117 0.3040 0.0697 -1.6118 1.3448
nancy 0.1181 -0.8811 0.7543 0.0122 -0.1470 0.2126 -0.1294 -0.6187 0.5841 -0.0222 -0.1098 0.0695 -0.8604 -14.8097 14.9397 -0.6076 -12.3969 12.3216
nantes -0.5317 -1.5223 0.4626 -0.0264 -0.2378 0.1001 0.5714 -0.1665 1.2975 0.0941 -0.0252 0.2218 1.1740 -37.7056 30.1182 0.8016 -33.1956 24.4721
nice 0.0356 -0.2189 0.2942 -0.0120 -0.1745 0.0826 0.0160 -0.0273 0.0856 -0.0009 -0.0202 0.0146 0.0261 -1.6473 1.9673 0.1223 -3.9662 3.6055
paris 0.3847 0.0412 0.7062 -0.0051 -0.0911 0.0446 -0.1070 -0.3123 0.1026 0.1145 0.0702 0.1609 0.0174 -0.4447 0.8033 -0.2977 -0.9758 0.4130
rennes 0.5019 -0.2924 1.4682 0.0310 -0.1257 0.2520 -0.3785 -1.1492 0.4282 0.0874 -0.0279 0.2188 -0.8990 -10.7549 11.4079 -0.9208 -11.3317 13.0179
rouen 0.2881 -0.3494 0.7971 0.0037 -0.0865 0.0987 -0.0610 -0.5837 0.5376 0.1611 0.0202 0.2984 0.2406 -1.1748 3.7754 -0.1352 -1.7612 2.5042
strasbourg -0.1851 -0.9156 0.4942 -0.0025 -0.1357 0.1295 0.3104 -0.1582 0.8403 0.0695 -0.0301 0.1770 1.2069 -22.0654 22.0152 0.9438 -16.2076 17.0276
toulouse -0.0198 -0.6370 0.6738 0.0040 -0.1323 0.1944 0.1078 -0.4509 0.7389 0.0955 0.0034 0.1829 0.8376 -9.4552 12.3973 0.5584 -9.1326 13.7519

Forestplot

Controlled Direct Effect

Pure Indirect Effect

Reference Interaction

Reference Interaction