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

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[[6]]$tempmax[which(is.na(villes_s[[6]]$tempmax))] <- villes_s[[6]]$tempmaxmoy7j[which(is.na(villes_s[[6]]$tempmax))]
villes_s[[6]]$tempmin[which(is.na(villes_s[[6]]$tempmin))]<-mean(villes_s[[6]]$tempmin)

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


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

for (i in 1:length(villes_s)){
  trshld95[i]<-quantile(villes_s[[i]]$tempmax, probs=c(0.95),na.rm=TRUE)
  for (r in 2:nrow(villes_s[[i]])){
  villes_s[[i]]$heat_wave[r]<-  if (villes_s[[i]]$tempmax[r] >= trshld95[i] & villes_s[[i]]$tempmax[r-1] >= trshld95[i]) 1 else 0
  }}

## reducing database only for years from 2000 to 2015
for (i in 1:length(villes_s)){
  villes_s[[i]]<-villes_s[[i]][villes_s[[i]]$Dates<="2016-01-01",]
  }

 

Number of Heat Wave for each city, for the period from 2000 to 2015:

Heat wave=0 Heat wave=1
bordeaux 1918 33
clermont 1908 43
grenoble 1896 55
lehavre 1922 29
lyon 1897 54
marseille 1914 37
montpellier 1919 32
nancy 1904 47
nantes 1917 34
nice 1910 41
paris 1914 37
rennes 1918 33
rouen 1916 35
strasbourg 1905 46
toulouse 1917 34

Creating Ozone binary variable

Number of O3>90 perc episodes per city:

O3>90perc=0 O3>90perc=1
bordeaux 1756 196
clermont 1756 196
grenoble 1756 196
lehavre 1756 196
lyon 1756 196
marseille 1757 195
montpellier 1756 196
nancy 1756 196
nantes 1756 196
nice 1757 195
paris 1756 196
rennes 1756 196
rouen 1756 196
strasbourg 1757 195
toulouse 1756 196

Analysis On Non-Accidental Mortality

Continuous Ozone

Code:

# some data management
 for (i in 1:length(villes_s)){
   villes_s[[i]]$time<-1:nrow(villes_s[[i]]) # Create a new variable for time 
 } 
   

    
  # start bootstrap
     
  nreps=800 #number of bootstrap reps
  results<-matrix(NA,nrow = 15,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(3*length(time)/122))+Jours+Vacances+hol,data=dat_boot,family=quasipoisson)
      
      # mediator model - multinomial linear regression
      m.mediator<-lm(o3~heat_wave+no2moy+ns(time,df=round(3*length(time)/122))+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) #PMM 
                   
      
    } #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

Binary Ozone

Code:

  # start bootstrap
     
  nreps=800 #number of bootstrap reps
  results_bin<-matrix(NA,nrow = 15,ncol=10)
  colnames(results_bin)<-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_90+ns(time,df=round(3*length(time)/122))+Jours+Vacances+hol,data=dat_boot,family=quasipoisson)
      
      # mediator model - multinomial linear regression
      m.mediator<-lm(o3~heat_wave+no2moy+ns(time,df=round(3*length(time)/122))+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) #PMM 
                   
      
    } #end bootstrap

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

Results

Results For Continuous Ozone:

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
bordeaux 1.3702 1.2006 1.5496 1.0283 1.0057 1.0564 0.0964 0.0175 0.2079
clermont 1.1967 1.0192 1.3971 1.0207 0.9924 1.0517 0.1122 -0.0866 0.4743
grenoble 1.0659 0.9652 1.1766 1.0373 1.0106 1.0652 0.3721 -0.0768 2.2157
lehavre 1.2033 0.9950 1.4297 1.0359 0.9866 1.0859 0.1741 -0.0848 0.9293
lyon 1.2010 1.0772 1.3482 1.0362 1.0202 1.0559 0.1785 0.0900 0.3748
marseille 1.0809 0.9985 1.1681 1.0045 0.9985 1.0137 0.0583 -0.0629 0.3570
montpellier 1.0616 0.9091 1.1997 1.0040 0.9966 1.0159 0.0430 -0.3535 0.6719
nancy 1.0816 0.9236 1.2502 1.0022 0.9740 1.0323 0.0209 -1.1772 1.3248
nantes 1.2275 1.0578 1.4100 1.0630 1.0240 1.1092 0.2552 0.0953 0.5908
nice 1.0803 0.9670 1.2030 1.0076 1.0001 1.0204 0.0896 -0.3230 0.6701
paris 1.5090 1.2747 1.7860 1.0540 1.0334 1.0835 0.1405 0.0901 0.2105
rennes 1.0941 0.9126 1.3202 1.0833 1.0204 1.1536 0.4703 -1.2430 2.6851
rouen 1.3343 1.1909 1.4943 1.0627 1.0226 1.1061 0.2020 0.0786 0.3505
strasbourg 1.3256 1.1942 1.4845 1.0221 1.0009 1.0459 0.0823 0.0033 0.1839
toulouse 1.1182 0.9945 1.2593 1.0397 1.0186 1.0665 0.2655 0.1126 0.9686 
write_xlsx(res,"C:/Users/Anna/Dropbox/Heat wave and O3/95 perc/mortality/nonlinearity/results_continuous.xlsx")

Results For Binary Ozone:

res_bin<-data.frame(results_bin)

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

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

kable(res_bin)%>%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
bordeaux 1.3420 1.1933 1.5418 9.4435 2.4144 43.3369 0.9709 0.8390 0.9951
clermont 1.1488 0.9859 1.3195 11.1484 2.6048 76.0637 0.9888 0.9085 1.0006
grenoble 1.0866 0.9798 1.1983 3.8783 0.7263 21.9007 0.9782 0.4689 1.2497
lehavre 1.2222 1.0277 1.4617 2.9794 0.3880 22.8593 0.9462 -0.9066 3.2414
lyon 1.2273 1.0986 1.3620 2.7301 0.8683 8.7584 0.9081 -0.0369 0.9884
marseille 1.0796 0.9992 1.1572 1.2078 0.9409 1.8339 0.7607 -0.4090 1.1529
montpellier 1.0650 0.9239 1.2124 1.0966 0.7500 1.8226 0.7841 -3.8483 3.7360
nancy 1.0793 0.9247 1.2691 1.6065 0.3394 8.1996 0.9703 -0.5560 2.8744
nantes 1.2384 1.0722 1.4041 17.7760 2.0582 194.6154 0.9895 0.8471 0.9993
nice 1.0665 0.9601 1.1826 1.9172 1.1626 3.8908 0.9446 0.6585 1.0578
paris 1.5039 1.2716 1.8484 18.3586 8.0459 54.9329 0.9815 0.9568 0.9933
rennes 1.0944 0.9039 1.2923 167.3239 6.2881 4905.3791 0.9996 0.9756 1.0006
rouen 1.3367 1.2053 1.5075 24.9632 2.2676 283.4419 0.9897 0.8441 0.9992
strasbourg 1.2881 1.1637 1.4357 8.9491 2.4775 38.6887 0.9727 0.8593 0.9951
toulouse 1.1361 1.0051 1.2752 3.8333 1.4378 14.0706 0.9629 0.7633 0.9985 
write_xlsx(res_bin,"C:/Users/Anna/Dropbox/Heat wave and O3/95 perc/mortality/nonlinearity/results_binary.xlsx")

Forestplot

Natural Direct Effect

Natural Indirect Effect

Analysis On Cardiovascual Mortality

Continuous Ozone

Code:

# some data management
 for (i in 1:length(villes_s)){
   villes_s[[i]]$time<-1:nrow(villes_s[[i]]) # Create a new variable for time 
 } 
   

    
  # start bootstrap
     
  nreps=800 #number of bootstrap reps
  results<-matrix(NA,nrow = 15,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(cv_tot~heat_wave+o3+ns(time,df=round(3*length(time)/122))+Jours+Vacances+hol,data=dat_boot,family=quasipoisson)
      
      # mediator model - multinomial linear regression
      m.mediator<-lm(o3~heat_wave+no2moy+ns(time,df=round(3*length(time)/122))+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) #PMM 
                   
      
    } #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

Binary Ozone

Code:

  # start bootstrap
     
  nreps=800 #number of bootstrap reps
  results_bin<-matrix(NA,nrow = 15,ncol=10)
  colnames(results_bin)<-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_90+ns(time,df=round(3*length(time)/122))+Jours+Vacances+hol,data=dat_boot,family=quasipoisson)
      
      # mediator model - multinomial linear regression
      m.mediator<-lm(o3~heat_wave+no2moy+ns(time,df=round(3*length(time)/122))+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) #PMM 
                   
      
    } #end bootstrap

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

Results

Results For Continuous Ozone:

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
bordeaux 1.1805 0.9834 1.4122 1.0394 0.9977 1.0810 0.1991 -0.0609 0.9903
clermont 1.2642 0.9993 1.5720 1.0027 0.9469 1.0547 0.0116 -0.4733 0.5567
grenoble 1.1760 0.9485 1.4136 1.0626 1.0141 1.1143 0.2885 0.0187 1.7133
lehavre 1.0103 0.6769 1.4242 1.1315 1.0264 1.2370 0.5079 -5.8677 11.1371
lyon 1.3128 1.1028 1.5516 1.0350 1.0017 1.0686 0.1278 0.0064 0.3323
marseille 1.0471 0.9107 1.1931 1.0054 0.9977 1.0200 0.0556 -0.9545 0.9519
montpellier 0.8265 0.5738 1.0734 1.0057 0.9917 1.0287 -0.0218 -0.4196 0.2077
nancy 1.0220 0.7817 1.3316 1.0038 0.9473 1.0612 0.0183 -2.0764 2.1590
nantes 1.2669 1.0330 1.5491 1.0713 0.9973 1.1574 0.2551 -0.0084 0.7314
nice 1.0348 0.8523 1.2491 1.0172 1.0020 1.0427 0.1150 -2.8257 2.7074
paris 1.5661 1.2879 1.9031 1.0508 1.0251 1.0884 0.1273 0.0705 0.2027
rennes 1.2515 0.8683 1.7779 1.0440 0.9363 1.1733 0.1617 -0.7664 2.3572
rouen 1.1475 0.9002 1.4342 1.0878 1.0084 1.1739 0.3900 -0.9800 1.9698
strasbourg 1.3251 1.0759 1.5663 0.9888 0.9448 1.0298 -0.0542 -0.3668 0.1540
toulouse 1.1202 0.9170 1.3379 1.0629 1.0221 1.1086 0.3504 -1.7777 1.6377 
write_xlsx(res,"C:/Users/Anna/Dropbox/Heat wave and O3/95 perc/mortality/nonlinearity/results_cv_continuous.xlsx")

Results For Binary Ozone:

res_bin<-data.frame(results_bin)

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

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

kable(res_bin)%>%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
bordeaux 1.3533 1.1891 1.5207 9.3289 2.3915 40.9740 0.9701 0.8380 0.9942
clermont 1.1491 0.9872 1.3053 11.6916 2.3390 70.3759 0.9895 0.8922 1.0008
grenoble 1.0875 0.9792 1.2004 4.2208 0.7737 22.6770 0.9808 0.6842 1.4558
lehavre 1.2280 1.0238 1.4522 3.3136 0.3458 26.7763 0.9553 -0.8958 3.0124
lyon 1.2228 1.0844 1.3614 2.8155 0.8946 9.2258 0.9141 0.0775 0.9906
marseille 1.0797 0.9974 1.1637 1.1800 0.9393 1.8067 0.7390 -0.5087 1.3177
montpellier 1.0637 0.9198 1.2079 1.0984 0.6775 1.8126 0.7971 -3.0240 4.1783
nancy 1.0693 0.9033 1.2433 1.6762 0.3376 8.2908 0.9765 -0.7752 2.3436
nantes 1.2341 1.0730 1.4170 17.2521 1.4197 220.2446 0.9887 0.7407 0.9993
nice 1.0641 0.9414 1.1893 1.8767 1.1730 4.0518 0.9440 0.6760 1.0877
paris 1.5075 1.2614 1.8101 18.3054 7.1884 52.8454 0.9811 0.9530 0.9934
rennes 1.0895 0.9134 1.3113 148.8960 4.6891 6343.3455 0.9996 0.9733 1.0011
rouen 1.3441 1.2049 1.4987 24.7026 2.6714 278.8204 0.9895 0.8533 0.9991
strasbourg 1.2901 1.1681 1.4387 8.1450 1.9369 37.7110 0.9698 0.7751 0.9948
toulouse 1.1311 1.0093 1.2866 3.8985 1.2894 13.4205 0.9640 0.6982 0.9987 
write_xlsx(res_bin,"C:/Users/Anna/Dropbox/Heat wave and O3/95 perc/mortality/nonlinearity/results_cv_binary.xlsx")

Forestplot

Natural Direct Effect

Natural Indirect Effect

Analysis On Respiratory Mortality

Continuous Ozone

Code:

# some data management
 for (i in 1:length(villes_s)){
   villes_s[[i]]$time<-1:nrow(villes_s[[i]]) # Create a new variable for time 
 } 
   

    
  # start bootstrap
     
  nreps=800 #number of bootstrap reps
  results<-matrix(NA,nrow = 15,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(3*length(time)/122))+Jours+Vacances+hol,data=dat_boot,family=quasipoisson)
      
      # mediator model - multinomial linear regression
      m.mediator<-lm(o3~heat_wave+no2moy+ns(time,df=round(3*length(time)/122))+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) #PMM 
                   
      
    } #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

Binary Ozone

Code:

  # start bootstrap
     
  nreps=800 #number of bootstrap reps
  results_bin<-matrix(NA,nrow = 15,ncol=10)
  colnames(results_bin)<-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_90+ns(time,df=round(3*length(time)/122))+Jours+Vacances+hol,data=dat_boot,family=quasipoisson)
      
      # mediator model - multinomial linear regression
      m.mediator<-lm(o3~heat_wave+no2moy+ns(time,df=round(3*length(time)/122))+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) #PMM 
                   
      
    } #end bootstrap

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

Results

Results For Continuous Ozone:

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
bordeaux 1.3781 1.2106 1.5562 1.0277 1.0032 1.0529 0.0936 0.0108 0.1986
clermont 1.1896 1.0327 1.3661 1.0226 0.9932 1.0528 0.1205 -0.0498 0.5092
grenoble 1.0669 0.9664 1.1859 1.0377 1.0136 1.0681 0.3669 0.1023 2.2203
lehavre 1.2008 0.9921 1.4269 1.0351 0.9824 1.0881 0.1682 -0.1246 0.8268
lyon 1.1964 1.0673 1.3440 1.0367 1.0197 1.0559 0.1822 0.0947 0.3970
marseille 1.0817 0.9987 1.1751 1.0046 0.9986 1.0133 0.0550 -0.0456 0.3438
montpellier 1.0667 0.9190 1.2004 1.0038 0.9968 1.0159 0.0390 -0.3859 0.6752
nancy 1.0775 0.9213 1.2534 1.0024 0.9752 1.0296 0.0120 -1.2266 1.1804
nantes 1.2277 1.0703 1.3993 1.0627 1.0219 1.1063 0.2574 0.0909 0.5600
nice 1.0738 0.9620 1.2006 1.0080 1.0003 1.0199 0.0934 -0.5800 0.8159
paris 1.5000 1.2509 1.8075 1.0541 1.0326 1.0845 0.1411 0.0905 0.2054
rennes 1.0990 0.9055 1.3477 1.0830 1.0201 1.1535 0.4578 -0.8879 2.9449
rouen 1.3354 1.1975 1.4988 1.0627 1.0231 1.1108 0.2008 0.0732 0.3459
strasbourg 1.3260 1.1959 1.4582 1.0236 1.0006 1.0484 0.0885 0.0027 0.1966
toulouse 1.1257 0.9926 1.2910 1.0401 1.0162 1.0681 0.2662 0.0948 1.0090 
write_xlsx(res,"C:/Users/Anna/Dropbox/Heat wave and O3/95 perc/mortality/nonlinearity/results_respi_continuous.xlsx")

Results For Binary Ozone:

res_bin<-data.frame(results_bin)

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

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

kable(res_bin)%>%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
bordeaux 1.3491 1.1911 1.5197 9.0932 2.5691 40.8435 0.9697 0.8492 0.9945
clermont 1.1399 0.9803 1.3202 11.5198 2.2758 55.9332 0.9896 0.8985 1.0017
grenoble 1.0888 0.9822 1.1930 4.0500 0.7673 20.5685 0.9794 0.5587 1.4068
lehavre 1.2243 1.0046 1.4510 2.9230 0.3575 29.1886 0.9497 -1.9935 2.5076
lyon 1.2218 1.1019 1.3708 2.7138 0.9394 8.6337 0.9115 0.3711 0.9890
marseille 1.0756 0.9970 1.1673 1.1905 0.9441 1.8660 0.7593 -0.7839 1.8365
montpellier 1.0634 0.9294 1.2016 1.1182 0.7543 1.7685 0.8067 -1.7927 4.8911
nancy 1.0731 0.9035 1.2513 1.7144 0.3004 9.0761 0.9761 -0.3275 2.3438
nantes 1.2424 1.0776 1.4143 14.7298 1.4829 202.6463 0.9872 0.6969 0.9994
nice 1.0656 0.9569 1.1754 1.8797 1.1859 4.1339 0.9420 0.6876 1.0666
paris 1.5163 1.2742 1.8083 18.1280 7.8118 55.6802 0.9812 0.9553 0.9936
rennes 1.0922 0.9264 1.3178 183.3340 6.5026 4203.8244 0.9996 0.9738 1.0006
rouen 1.3400 1.1999 1.5048 25.7431 2.9641 296.4434 0.9901 0.8748 0.9993
strasbourg 1.2940 1.1633 1.4296 8.4227 2.0242 40.1522 0.9703 0.8063 0.9948
toulouse 1.1410 1.0050 1.2765 3.9950 1.2861 13.5976 0.9642 0.6917 0.9992 
write_xlsx(res_bin,"C:/Users/Anna/Dropbox/Heat wave and O3/95 perc/mortality/nonlinearity/results_respi_binary.xlsx")

Forestplot

Natural Direct Effect

Natural Indirect Effect