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("respi","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(respi~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.0072 0.8596 1.1717 1.0514 0.9897 1.1161 0.4908 -6.7967 6.0548
bordeaux 1.0816 0.9882 1.1802 1.0167 0.9832 1.0483 0.1772 -0.3103 1.1872
clermont 1.0759 0.9277 1.2570 1.0017 0.9494 1.0517 0.0161 -1.8876 3.2967
grenoble 1.0861 0.9686 1.2067 1.0174 0.9840 1.0508 0.1707 -0.5268 1.2597
lehavre 1.0913 0.9377 1.2441 0.9934 0.9410 1.0462 -0.1030 -3.0525 2.7846
lille 1.0226 0.9455 1.1082 1.0104 0.9791 1.0399 0.1846 -3.6692 5.2497
lyon 1.0319 0.9445 1.1170 0.9833 0.9558 1.0089 -0.2158 -7.5167 7.3884
marseille 0.9791 0.9101 1.0438 1.0079 1.0000 1.0200 -0.1200 -3.7823 3.3794
montpellier 0.9794 0.8699 1.0965 0.9999 0.9957 1.0029 0.0001 -0.3717 0.2641
nancy 1.0120 0.9005 1.1232 1.0182 0.9823 1.0592 0.1712 -6.0290 5.5112
nantes 1.0478 0.9271 1.1621 1.0167 0.9776 1.0622 0.1906 -3.9127 3.0593
nice 1.1002 1.0106 1.1978 0.9988 0.9884 1.0043 -0.0135 -0.2152 0.0862
paris 1.0458 0.9954 1.0956 1.0015 0.9834 1.0198 0.0319 -0.9029 0.8586
rennes 1.0624 0.8984 1.2339 1.0292 0.9731 1.0981 0.2529 -2.3993 3.8309
rouen 0.9432 0.8534 1.0395 0.9815 0.9341 1.0312 0.2100 -1.7175 2.2413
strasbourg 1.0280 0.9285 1.1389 1.0523 1.0124 1.0957 0.5691 -2.8958 5.0313
toulouse 1.1650 1.0618 1.2769 0.9963 0.9712 1.0188 -0.0274 -0.2626 0.1582

Forestplot

Natural Direct Effect

Natural Indirect Effect