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

seuil_IBMax<-c(33,35,34,34,33,33,34,35,35,34,34,31,31,34,33,34,36)
seuil_IBMin<-c(18,21,19,19,18,18,20,24,22,18,20,24,21,19,19,19,21)

## imputation NA pour Marseille
villes_s[[8]]$tempmax[which(is.na(villes_s[[8]]$tempmax))] <- villes_s[[8]]$tempmaxmoy7j[which(is.na(villes_s[[8]]$tempmax))]
villes_s[[8]]$tempmin[which(is.na(villes_s[[8]]$tempmin))]<-mean(villes_s[[8]]$tempmin)

## imputation NA pour rouen
villes_s[[15]]$tempmax[which(is.na(villes_s[[15]]$tempmax))] <- villes_s[[15]]$tempmaxmoy7j[which(is.na(villes_s[[15]]$tempmax))]


## Creation moyenne glissante sur 3 jours
for (i in 1:length(villes_s)){
villes_s[[i]]$tmax3<-as.vector(filter(villes_s[[i]]$tempmax,sides = 1, filter = c(1,1,1)/3))
villes_s[[i]]$tmin3<-as.vector(filter(villes_s[[i]]$tempmin,sides = 1, filter = c(1,1,1)/3))
}

## Creation Varibale Heat wave
for (i in 1:length(villes_s)){
  for (r in 3:nrow(villes_s[[i]])){
  villes_s[[i]]$heat_wave[r]<-  if (villes_s[[i]]$tmax3[r] > seuil_IBMax[i] & villes_s[[i]]$tmin3[r] > seuil_IBMin[i]) 1 else 0 
  }}
)

 

Number of Heat Wave for each city:

Heat wave=0 Heat wave=1
bm 2190 4
bordeaux 2180 14
clermont 2163 31
grenoble 2193 1
lehavre 2192 2
lille 2190 4
lyon 2133 61
marseille 2190 4
montpellier 2191 3
nancy 2179 15
nantes 2188 6
nice 2173 21
paris 2177 17
rennes 2189 5
rouen 2190 4
strasbourg 2185 9
toulouse 2181 13 
# 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 On Non-Accidental Mortality

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 NA NA NA NA NA NA NA NA NA
bordeaux NA NA NA NA NA NA NA NA NA
clermont 1.4037 1.1026 1.6904 1.0137 0.9794 1.0464 0.0450 -0.0834 0.2170
grenoble 1.1751 1.0363 1.3665 1.0372 1.0048 1.0775 0.1990 0.0185 0.5671
lehavre NA NA NA NA NA NA NA NA NA
lille NA NA NA NA NA NA NA NA NA
lyon 1.0966 0.9828 1.2092 1.0113 0.9992 1.0273 0.1067 -0.2470 0.6370
marseille NA NA NA 0.9907 0.9809 0.9994 NA NA NA
montpellier 0.6503 0.4798 0.7925 0.9999 0.9912 1.0087 0.0000 -0.0147 0.0151
nancy 1.2331 0.9132 1.5811 0.9924 0.9591 1.0240 -0.0372 -0.6176 0.7326
nantes NA NA NA NA NA NA NA NA NA
nice 1.1098 0.9064 1.3163 1.0008 0.9870 1.0173 0.0044 -0.3411 0.3990
paris 1.3261 1.2586 1.3889 1.0050 0.9863 1.0342 0.0403 -0.0136 0.1390
rennes NA NA NA NA NA NA NA NA NA
rouen NA NA NA NA NA NA NA NA NA
strasbourg 1.5499 1.1788 1.9308 1.0339 1.0064 1.0842 0.0925 0.0203 0.2073
toulouse 1.1521 0.8331 1.4758 1.0405 1.0106 1.1201 0.1519 -2.6395 1.2595

Forestplot

Natural Direct Effect

Natural Indirect Effect