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.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
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:

## create dataframe without missing values
for (i in 1:length(villes_s)){
  villes_s[[i]]<-na.omit(villes_s[[i]][, c("cardiaque","o3","Jours","Vacances","hol","heat_wave","no2moy")])
  #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(cardiaque~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 1.0113 0.8811 1.1590 0.9861 0.9336 1.0371 -0.0502 -6.2040 5.7065
bordeaux 0.9794 0.9081 1.0531 1.0132 0.9838 1.0436 -0.1550 -6.1962 6.4993
clermont 1.0612 0.9670 1.1553 1.0001 0.9707 1.0314 0.0017 -1.4140 1.7256
grenoble 0.9912 0.9008 1.0856 1.0269 1.0007 1.0535 0.3693 -6.9107 9.0630
lehavre 0.9829 0.8741 1.1119 0.9966 0.9464 1.0477 -0.0365 -9.1960 6.1218
lille 0.9800 0.9074 1.0571 1.0053 0.9777 1.0376 -0.0834 -5.3120 4.3012
lyon 0.9721 0.9075 1.0288 0.9952 0.9741 1.0151 0.1064 -1.8714 3.2360
marseille 0.9691 0.9106 1.0298 1.0023 0.9944 1.0114 -0.0428 -1.2582 1.0739
montpellier 1.1522 1.0269 1.2754 1.0001 0.9973 1.0040 0.0006 -0.0280 0.0423
nancy 0.9184 0.8218 1.0183 0.9970 0.9631 1.0270 0.0288 -0.7384 1.2511
nantes 1.0370 0.9469 1.1286 1.0040 0.9724 1.0381 0.0540 -3.2350 3.2141
nice 0.9244 0.8358 1.0179 0.9993 0.9905 1.0034 0.0080 -0.0786 0.2041
paris 1.0174 0.9741 1.0601 0.9920 0.9773 1.0053 -0.1405 -8.1488 7.5663
rennes 0.9938 0.8321 1.1842 1.0033 0.9532 1.0584 -0.0070 -4.3666 3.2939
rouen 1.0498 0.9678 1.1447 0.9868 0.9503 1.0259 -0.1791 -4.9429 5.0951
strasbourg 1.0198 0.9400 1.1052 0.9955 0.9624 1.0291 -0.1160 -6.1755 4.4193
toulouse 1.0164 0.9469 1.0919 0.9979 0.9770 1.0167 -0.0346 -3.0517 2.8635

Forestplot

Natural Direct Effect

Natural Indirect Effect