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)
  }

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

Creating Heat Wave variable

trshld95<-c()

for (i in 1:length(villes_s)){
  trshld95[i]<-quantile(villes_s[[i]]$tempmax, probs=c(0.975),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 2143 53
bordeaux 2143 53
clermont 2142 54
grenoble 2141 55
lehavre 2142 54
lille 2143 53
lyon 2142 54
marseille 2142 54
montpellier 2145 51
nancy 2142 54
nantes 2141 55
nice 2142 54
paris 2144 52
rennes 2142 54
rouen 2143 53
strasbourg 2143 53
toulouse 2141 55 
# 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:

## 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",]
  villes_s[[i]]$time<-1:nrow(villes_s[[i]])
} 
   

    
  # start bootstrap
     
  nreps=1000 #number of bootstrap reps
  results<-matrix(NA,nrow = 17,ncol=10)
  colnames(results)<-c("city","NDE","NDE_2.5","NDE_97.5","NIE","NIE_2.5","NIE_97.5","PMM","PMM_2.5","PMM_97.5")
  
for (i in (1:length(villes_s))){
  boot.nde<- rep(NA,nreps)
  boot.nie <- rep(NA,nreps)
  boot.pmm <- 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+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+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]
     
      beta1<-coef(m.mediator)[2]
      
      # estimate CDE and NIE
      
      boot.nde[rep] <- exp(theta1) #NDE
      
      boot.nie[rep] <-exp(theta2*beta1) #NIE 
      
      boot.pmm[rep] <-boot.nde[rep]*(boot.nie[rep]-1)/(boot.nde[rep]*boot.nie[rep]-1) #NIE 
                   
      
    } #end bootstrap

  results[i,1]<-cities[[i]]
  results[i,2]<-median(boot.nde, na.rm=T)
  results[i,3]<-quantile(boot.nde, 0.025, na.rm=T)
  results[i,4]<-quantile(boot.nde, 0.975, na.rm=T)
  results[i,5]<-median(boot.nie, na.rm=T)
  results[i,6]<-quantile(boot.nie, 0.025, na.rm=T)
  results[i,7]<-quantile(boot.nie, 0.975, na.rm=T)
  results[i,8]<-median(boot.pmm, na.rm=T)
  results[i,9]<-quantile(boot.pmm, 0.025, na.rm=T)
  results[i,10]<-quantile(boot.pmm, 0.975, na.rm=T)
  
     }
  
  # output CDE and NIE and confidence intervals

Results:

res<-data.frame(results)

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

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

kable(res)%>%kable_styling()%>%
  column_spec(1, bold = T, border_right = T)%>%
  column_spec(4, border_right = T)%>%
  column_spec(7, border_right = T)
city NDE NDE_2.5 NDE_97.5 NIE NIE_2.5 NIE_97.5 PMM PMM_2.5 PMM_97.5
bm 0.8798 0.7533 1.0320 1.0379 0.9897 1.0897 -0.2849 -4.0327 3.3879
bordeaux 1.2016 1.0471 1.3769 1.0310 0.9915 1.0754 0.1582 -0.0411 0.5046
clermont 1.2079 1.0502 1.4044 1.0136 0.9742 1.0510 0.0706 -0.1501 0.3630
grenoble 0.9727 0.8296 1.1328 1.0443 1.0054 1.0867 0.2775 -9.5640 9.5318
lehavre 1.0771 0.9051 1.2740 1.0659 1.0096 1.1359 0.4346 -1.8672 3.2290
lille 0.9906 0.8870 1.0974 1.0548 1.0243 1.0894 0.7334 -11.3172 10.4998
lyon 1.0717 0.9572 1.1813 1.0177 0.9981 1.0394 0.1961 -0.6175 1.6263
marseille 1.1312 1.0247 1.2441 1.0043 0.9998 1.0129 0.0361 -0.0035 0.1835
montpellier 0.9460 0.8144 1.1083 0.9997 0.9933 1.0034 0.0018 -0.3581 0.2450
nancy 0.9486 0.7972 1.1146 0.9943 0.9623 1.0292 0.0280 -1.8903 1.7853
nantes 1.0400 0.8429 1.2488 1.0419 1.0047 1.0878 0.3155 -3.7076 5.5010
nice 0.9870 0.8897 1.1113 1.0004 0.9940 1.0078 0.0011 -1.0567 0.8452
paris 1.0479 0.9897 1.1148 1.0376 1.0251 1.0534 0.4530 0.2315 1.2776
rennes 1.1250 0.7755 1.4903 1.0407 0.9802 1.1029 0.1660 -2.1066 2.4103
rouen 1.1513 0.9782 1.3183 1.0636 1.0162 1.1182 0.3252 0.0680 1.2485
strasbourg 1.1396 0.9813 1.2955 1.0398 1.0063 1.0784 0.2456 0.0054 1.1205
toulouse 1.1057 0.9549 1.2751 1.0290 1.0057 1.0562 0.2172 -0.8547 1.7425

Forestplot

Natural Direct Effect

Natural Indirect Effect