D-R Curves for Ozone and nocc mortality over the time with two-stage modeling - Basic and Stratification Models

Descriptive Statistics for O3

Descriptive Statistics for o3

Ville	Minimum	X1st.Qu	Median	Mean	X3rd.Qu	Max	NA.s
bordeaux	1.9154062	49.55208	69.45833	69.41858	88.12500	169.5000	NA
clermont	1.7321429	54.31250	72.48661	72.01823	89.87500	178.5625	144
grenoble	0.4375000	35.29129	65.78125	65.64074	91.66964	185.8750	NA
lehavre	2.1250000	56.68750	69.87500	69.16447	81.75000	184.5000	76
lille	0.1666667	41.00000	58.75000	58.65181	74.21250	207.1250	285
lyon	0.4375000	39.43799	64.94705	65.15421	89.04688	229.3750	8
marseille	3.1250000	54.31250	79.00000	77.03170	99.10938	195.0625	16
montpellier	7.5458333	61.66667	81.16250	81.22084	100.05417	200.8083	9
nancy	0.5000000	52.00000	68.25000	70.65748	87.37500	221.0000	170
nantes	1.1250000	52.93750	68.75000	70.20679	86.12500	216.5000	75
nice	14.5625000	60.50000	89.25000	87.41249	110.50000	206.3750	60
paris	0.8592206	36.77541	56.96864	58.49854	76.43874	216.3857	NA
rennes	1.5625000	50.75000	64.87500	65.52426	79.50000	213.5625	18
rouen	1.2972346	47.00000	63.15625	63.59236	78.75000	178.6607	16
strasbourg	0.0000000	36.62500	61.75000	63.87351	87.25000	206.5000	113
toulouse	1.6250000	56.33750	75.54063	75.36355	93.80938	180.9594	NA

Dose-Response Curves for o3 and Non-accidental Mortality, Basic Model

For Lag0, Lag1 and Lag2 exposure

# SET MODELS OBJECTS#
tot <- matrix(NA, nrow = 16, ncol = 10)
colnames(tot) <- c("Ville","Coefficient_L0","L0_CI95_Low","L0_CI95_High", "Coefficient_L1","L1_CI95_Low","L1_CI95_High","Coefficient_L2","L2_CI95_Low","L2_CI95_High")
poll<-list()
modtot<-list() 
pred<-list()
y<-matrix(NA,nrow=16,ncol=1)
S<-list()

# MODEL #
for(i in 1:length(villes)) {
  poll[[i]] <- onebasis(villes[[i]]$o3.x,"lin")
  pol<-poll[[i]]
  modtot[[i]]<-gam(nocc_tot~pol+ns(time,df=round(8*length(time)/365))+ns(tempmoy,df=6)+ns(temp3,df=6)+Jours+Vacances,data=villes[[i]],family=quasipoisson)
  est<-summary(modtot[[i]])
  tot[i, 1] <- as.character(villes_name[[i]])
  tot[i, 2] <- modtot[[i]]$coefficients[2]
  tot[i, 3] <- modtot[[i]]$coefficients[2]-(1.96*est$se[2])
  tot[i, 4] <- modtot[[i]]$coefficients[2]+(1.96*est$se[2])
  pred[[i]]<-crosspred(pol,modtot[[i]])
  y[i,] <- pred[[i]]$coef
  S[[i]] <- pred[[i]]$vcov
}


# SET MODELS OBJECTS#
poll_L1<-list()
modtot_L1<-list() 
pred_L1<-list()
y_L1<-matrix(NA,nrow=16,ncol=1)
S_L1<-list()

for(i in 1:length(villes)) {
  poll_L1[[i]] <- onebasis(villes[[i]]$o3_L1,"lin")
  pol_L1<-poll_L1[[i]]
  modtot_L1[[i]]<-gam(nocc_tot~pol_L1+ns(time,df=round(8*length(time)/365))+ns(tempmoy,df=6)+ns(temp3,df=6)+Jours+Vacances,data=villes[[i]],family=quasipoisson)
  est_L1<-summary(modtot_L1[[i]])
  tot[i, 5] <- modtot_L1[[i]]$coefficients[2]
  tot[i, 6] <- modtot_L1[[i]]$coefficients[2]-(1.96*est_L1$se[2])
  tot[i, 7] <- modtot_L1[[i]]$coefficients[2]+(1.96*est_L1$se[2])
  pred_L1[[i]]<-crosspred(pol_L1,modtot_L1[[i]])
  y_L1[i,] <- pred_L1[[i]]$coef
  S_L1[[i]] <- pred_L1[[i]]$vcov
  }

# SET MODELS OBJECTS#
poll_L2<-list()
modtot_L2<-list() 
pred_L2<-list()
y_L2<-matrix(NA,nrow=16,ncol=1)
S_L2<-list()

# MODEL #
for(i in 1:length(villes)) {
  poll_L2[[i]] <- onebasis(villes[[i]]$o3_L2,"lin")
  pol_L2<-poll_L2[[i]]
  modtot_L2[[i]]<-gam(nocc_tot~pol_L2+ns(time,df=round(8*length(time)/365))+ns(tempmoy,df=6)+ns(temp3,df=6)+Jours+Vacances,data=villes[[i]],family=quasipoisson)
  est_L2<-summary(modtot_L2[[i]])
  tot[i, 8] <- modtot_L2[[i]]$coefficients[2]
  tot[i, 9] <- modtot_L2[[i]]$coefficients[2]-(1.96*est_L2$se[2])
  tot[i, 10] <- modtot_L2[[i]]$coefficients[2]+(1.96*est_L2$se[2])
  pred_L2[[i]]<-crosspred(pol_L2,modtot_L2[[i]])
  y_L2[i,] <- pred_L2[[i]]$coef
  S_L2[[i]] <- pred_L2[[i]]$vcov
  }

## Warning in gam.fit3(x = args$X, y = args$y, sp = lsp, Eb = args$Eb, UrS =
## args$UrS, : Algorithm did not converge

Ville	Coefficient_L0	L0_CI95_Low	L0_CI95_High	Coefficient_L1	L1_CI95_Low	L1_CI95_High	Coefficient_L2	L2_CI95_Low	L2_CI95_High
bordeaux	1.00000	0.99951	1.00048	1.00039	0.99991	1.00087	1.00026	0.99979	1.00074
clermont	1.00039	0.99972	1.00106	1.00024	0.99957	1.00092	1.00077	1.00010	1.00144
grenoble	1.00003	0.99950	1.00056	0.99964	0.99910	1.00017	1.00003	0.99950	1.00057
lehavre	1.00023	0.99954	1.00092	1.00036	0.99967	1.00105	1.00028	0.99961	1.00095
lille	1.00023	0.99989	1.00058	1.00032	0.99998	1.00066	1.00013	0.99979	1.00047
lyon	1.00017	0.99981	1.00053	1.00049	1.00013	1.00084	1.00048	1.00012	1.00083
marseille	1.00060	1.00024	1.00096	1.00037	1.00001	1.00073	1.00048	1.00012	1.00084
montpellier	1.00058	0.99994	1.00123	1.00105	1.00040	1.00170	1.00081	1.00017	1.00147
nancy	1.00020	0.99960	1.00079	1.00041	0.99981	1.00101	1.00093	1.00033	1.00153
nantes	1.00073	1.00020	1.00126	1.00129	1.00077	1.00182	1.00109	1.00057	1.00161
nice	1.00070	1.00021	1.00119	1.00097	1.00048	1.00147	1.00093	1.00044	1.00142
paris	1.00016	0.99995	1.00036	1.00010	0.99990	1.00030	0.99994	0.99974	1.00014
rennes	0.99995	0.99911	1.00078	0.99928	0.99844	1.00011	0.99921	0.99838	1.00003
rouen	1.00058	1.00004	1.00111	1.00083	1.00030	1.00137	1.00072	1.00020	1.00124
strasbourg	1.00002	0.99953	1.00051	1.00012	0.99963	1.00060	1.00020	0.99971	1.00068
toulouse	1.00050	0.99995	1.00104	1.00075	1.00022	1.00129	1.00087	1.00033	1.00141

Meta-Analysis results

meta <- mvmeta(y,S,method="ml")
summary(meta)

## Call:  mvmeta(formula = y ~ 1, S = S, method = "ml")
## 
## Univariate random-effects meta-analysis
## Dimension: 1
## Estimation method: ML
## 
## Fixed-effects coefficients
##              Estimate  Std. Error       z  Pr(>|z|)  95%ci.lb  95%ci.ub     
## (Intercept)    0.0003      0.0001  4.8083    0.0000    0.0002    0.0004  ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Between-study random-effects (co)variance components
##   Std. Dev
##     0.0001
## 
## Univariate Cochran Q-test for heterogeneity:
## Q = 17.1344 (df = 15), p-value = 0.3109
## I-square statistic = 12.5%
## 
## 16 studies, 16 observations, 1 fixed and 1 random-effects parameters
##    logLik        AIC        BIC  
##  109.4054  -214.8109  -213.2657

meta_L1 <- mvmeta(y_L1,S_L1,method="ml")
summary(meta_L1)

## Call:  mvmeta(formula = y_L1 ~ 1, S = S_L1, method = "ml")
## 
## Univariate random-effects meta-analysis
## Dimension: 1
## Estimation method: ML
## 
## Fixed-effects coefficients
##              Estimate  Std. Error       z  Pr(>|z|)  95%ci.lb  95%ci.ub     
## (Intercept)    0.0004      0.0001  3.9235    0.0001    0.0002    0.0006  ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Between-study random-effects (co)variance components
##   Std. Dev
##     0.0003
## 
## Univariate Cochran Q-test for heterogeneity:
## Q = 49.0712 (df = 15), p-value = 0.0000
## I-square statistic = 69.4%
## 
## 16 studies, 16 observations, 1 fixed and 1 random-effects parameters
##    logLik        AIC        BIC  
##  100.1221  -196.2442  -194.6991

meta_L2 <- mvmeta(y_L2,S_L2,method="ml")
summary(meta_L2)

## Call:  mvmeta(formula = y_L2 ~ 1, S = S_L2, method = "ml")
## 
## Univariate random-effects meta-analysis
## Dimension: 1
## Estimation method: ML
## 
## Fixed-effects coefficients
##              Estimate  Std. Error       z  Pr(>|z|)  95%ci.lb  95%ci.ub     
## (Intercept)    0.0004      0.0001  4.2788    0.0000    0.0002    0.0006  ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Between-study random-effects (co)variance components
##   Std. Dev
##     0.0003
## 
## Univariate Cochran Q-test for heterogeneity:
## Q = 53.0972 (df = 15), p-value = 0.0000
## I-square statistic = 71.7%
## 
## 16 studies, 16 observations, 1 fixed and 1 random-effects parameters
##    logLik        AIC        BIC  
##  101.0112  -198.0223  -196.4772

Pooled Dose-Responses Curves For Basic Model

pol <- seq(bound_o3[1],bound_o3[2],length=220)  # BASIS USED TO PREDICT PM10, EQUAL TO THAT USED FOR ESTIMATION
pol <- onebasis(pol,"lin")

plot <- crosspred(pol,coef=coef(meta),vcov=vcov(meta),model.link="log",by=1)
plot(plot,"overall",ci="lines",col=1,lwd=3,ylab="RR",xlab="Ozone",main="Pooled DR curves for 2000-2015")

plot <- crosspred(pol,coef=coef(meta_L1),vcov=vcov(meta_L1),model.link="log",by=1)
plot(plot,"overall",ci="lines",col=1,lwd=3,ylab="RR",xlab="Ozone_LAG1",main="")

plot <- crosspred(pol,coef=coef(meta_L2),vcov=vcov(meta_L2),model.link="log",by=1)
plot(plot,"overall",ci="lines",col=1,lwd=3,ylab="RR",xlab="Ozone_LAG2",main="")

3 Periods Stratified analysis

For Lag0, Lag1 and Lag2 exposure

######################### DEFINE DIFFERENT DATA PERIOD #####################################

y2001_y2005<-list()
for(i in 1:length(villesG1)) {
  y2001_y2005[[i]]<-villesG1[[i]][which (villesG1[[i]]$annee %in% c("2001","2002","2003","2004","2005")),]
  y2001_y2005[[i]]$time<-1:nrow(y2001_y2005[[i]])
  y2001_y2005[[i]]<-transform(y2001_y2005[[i]], aug2003=format(Dates, "%Y-%m")=="2003-08")
  y2001_y2005[[i]]<-y2001_y2005[[i]][y2001_y2005[[i]]$aug2003==FALSE,]  ### take off the heat wave of 2003
}

filter<-dplyr::filter
y2006_y2010<-list()
for(i in 1:length(villes)) {
  y2006_y2010[[i]]<- villes[[i]] %>% filter(annee %in% c("2006","2007","2008","2009","2010"))
  y2006_y2010[[i]]$time<-1:nrow(y2006_y2010[[i]])
}

y2011_y2015<-list()
for(i in 1:length(villes)) {
  y2011_y2015[[i]]<- villes[[i]] %>% filter(annee %in% c("2011","2012","2013","2014","2015"))
  y2011_y2015[[i]]$time<-1:nrow(y2011_y2015[[i]])
}

## o3 limit values for each period and each city
bound_o3_0105<-as.data.frame(matrix(NA,nrow=length(villesG1),ncol=2))
for(i in 1:length(villesG1)) {
  bound_o3_0105[i,]<-c(min(y2001_y2005[[i]]$o3.x, na.rm=TRUE),max(y2001_y2005[[i]]$o3.x,na.rm=TRUE))
}  

bound_o3_0610<-as.data.frame(matrix(NA,nrow=length(villes),ncol=2))
for(i in 1:length(villes)) {
  bound_o3_0610[i,]<-c(min(y2006_y2010[[i]]$o3.x, na.rm=TRUE),max(y2006_y2010[[i]]$o3.x,na.rm=TRUE))
}  

bound_o3_1115<-as.data.frame(matrix(NA,nrow=length(villes),ncol=2))
for(i in 1:length(villes)) {
  bound_o3_1115[i,]<-c(min(y2011_y2015[[i]]$o3.x, na.rm=TRUE),max(y2011_y2015[[i]]$o3.x,na.rm=TRUE))
}  

############################### SET MODELS OBJECTS############################################

tot0105 <- matrix(NA, nrow = 13, ncol = 10)
colnames(tot0105) <- c("Ville","Coefficient_L0","L0_CI95_Low","L0_CI95_High", "Coefficient_L1","L1_CI95_Low","L1_CI95_High","Coefficient_L2","L2_CI95_Low","L2_CI95_High")
poll_0105<-list()
modtot0105<-list() 
pred0105<-list()
y0105<-matrix(NA,nrow=13,ncol=1)
S0105<-list()
st_err0105<-c()

poll_0105_L1<-list()
modtot0105_L1<-list() 
pred0105_L1<-list()
y0105_L1<-matrix(NA,nrow=13,ncol=1)
S0105_L1<-list()
st_err0105_L1<-c()

poll_0105_L2<-list()
modtot0105_L2<-list() 
pred0105_L2<-list()
y0105_L2<-matrix(NA,nrow=13,ncol=1)
S0105_L2<-list()
st_err0105_L2<-c()


tot0610 <- matrix(NA, nrow = 16, ncol = 10)
colnames(tot0610) <- c("Ville","Coefficient_L0","L0_CI95_Low","L0_CI95_High", "Coefficient_L1","L1_CI95_Low","L1_CI95_High","Coefficient_L2","L2_CI95_Low","L2_CI95_High")
poll_0610<-list()
modtot0610<-list() 
pred0610<-list()
y0610<-matrix(NA,nrow=16,ncol=1)
S0610<-list()
st_err0610<-c()

poll_0610_L1<-list()
modtot0610_L1<-list() 
pred0610_L1<-list()
y0610_L1<-matrix(NA,nrow=16,ncol=1)
S0610_L1<-list()
st_err0610_L1<-c()


poll_0610_L2<-list()
modtot0610_L2<-list() 
pred0610_L2<-list()
y0610_L2<-matrix(NA,nrow=16,ncol=1)
S0610_L2<-list()
st_err0610_L2<-c()


tot1115 <- matrix(NA, nrow = 16, ncol = 10)
colnames(tot1115) <- c("Ville","Coefficient_L0","L0_CI95_Low","L0_CI95_High", "Coefficient_L1","L1_CI95_Low","L1_CI95_High","Coefficient_L2","L2_CI95_Low","L2_CI95_High")
poll_1115<-list()
modtot1115<-list() 
pred1115<-list()
y1115<-matrix(NA,nrow=16,ncol=1)
S1115<-list()
st_err1115<-c()


poll_1115_L1<-list()
modtot1115_L1<-list() 
pred1115_L1<-list()
y1115_L1<-matrix(NA,nrow=16,ncol=1)
S1115_L1<-list()
st_err1115_L1<-c()

poll_1115_L2<-list()
modtot1115_L2<-list() 
pred1115_L2<-list()
y1115_L2<-matrix(NA,nrow=16,ncol=1)
S1115_L2<-list()
st_err1115_L2<-c()




#################################### MODELS ################################################
##### For First Period
# LAG0
for(i in 1:length(villesG1)) {
  poll_0105[[i]] <- onebasis(y2001_y2005[[i]]$o3.x,"lin")
  pol<-poll_0105[[i]]
  modtot0105[[i]]<-gam(nocc_tot~pol+ns(time,df=round(8*length(time)/365))+ns(tempmoy,df=6)+ns(temp3,df=6)+Jours+Vacances,data=y2001_y2005[[i]],family=quasipoisson)
  est<-summary(modtot0105[[i]])
  st_err0105[i]<-est$se[2]
  tot0105[i, 1] <- as.character(villes_name[[i]])
  tot0105[i, 2] <- modtot0105[[i]]$coefficients[2]
  tot0105[i, 3] <- modtot0105[[i]]$coefficients[2]-(1.96*est$se[2])
  tot0105[i, 4] <- modtot0105[[i]]$coefficients[2]+(1.96*est$se[2])
  pred0105[[i]]<-crosspred(pol,modtot0105[[i]])
  y0105[i,] <- pred0105[[i]]$coef
  S0105[[i]] <- pred0105[[i]]$vcov
  }

# LAG1
for(i in 1:length(villesG1)) {
  poll_0105_L1[[i]] <- onebasis(y2001_y2005[[i]]$o3_L1,"lin")
  pol_L1<-poll_0105_L1[[i]]
  modtot0105_L1[[i]]<-gam(nocc_tot~pol_L1+ns(time,df=round(8*length(time)/365))+ns(tempmoy,df=6)+ns(temp3,df=6)+Jours+Vacances,data=y2001_y2005[[i]],family=quasipoisson)
  est_L1<-summary(modtot0105_L1[[i]])
  st_err0105_L1[i]<-est_L1$se[2]
  tot0105[i, 5] <- modtot0105_L1[[i]]$coefficients[2]
  tot0105[i, 6] <- modtot0105_L1[[i]]$coefficients[2]-(1.96*est_L1$se[2])
  tot0105[i, 7] <- modtot0105_L1[[i]]$coefficients[2]+(1.96*est_L1$se[2])
  pred0105_L1[[i]]<-crosspred(pol_L1,modtot0105_L1[[i]])
  y0105_L1[i,] <- pred0105_L1[[i]]$coef
  S0105_L1[[i]] <- pred0105_L1[[i]]$vcov
  }

# LAG2
for(i in 1:length(villesG1)) {
  poll_0105_L2[[i]] <- onebasis(y2001_y2005[[i]]$o3_L2,"lin")
  pol_L2<-poll_0105_L2[[i]]
  modtot0105_L2[[i]]<-gam(nocc_tot~pol_L2+ns(time,df=round(8*length(time)/365))+ns(tempmoy,df=6)+ns(temp3,df=6)+Jours+Vacances,data=y2001_y2005[[i]],family=quasipoisson)
  est_L2<-summary(modtot0105_L2[[i]])
  st_err0105_L2[i]<-est_L2$se[2]
  tot0105[i, 8] <- modtot0105_L2[[i]]$coefficients[2]
  tot0105[i, 9] <- modtot0105_L2[[i]]$coefficients[2]-(1.96*est_L2$se[2])
  tot0105[i, 10] <- modtot0105_L2[[i]]$coefficients[2]+(1.96*est_L2$se[2])
  pred0105_L2[[i]]<-crosspred(pol_L2,modtot0105_L2[[i]])
  y0105_L2[i,] <- pred0105_L2[[i]]$coef
  S0105_L2[[i]] <- pred0105_L2[[i]]$vcov
  }

### Second Period
# L0
for(i in 1:length(villes)) {
  poll_0610[[i]] <- onebasis(y2006_y2010[[i]]$o3.x,"lin")
  pol<-poll_0610[[i]]
  modtot0610[[i]]<-gam(nocc_tot~pol+ns(time,df=round(8*length(time)/365))+ns(tempmoy,df=6)+ns(temp3,df=6)+Jours+Vacances,data=y2006_y2010[[i]],family=quasipoisson)
  est<-summary(modtot0610[[i]])
  st_err0610[i]<-est$se[2]
  tot0610[i, 1] <- as.character(villes_name[[i]])
  tot0610[i, 2] <- modtot0610[[i]]$coefficients[2]
  tot0610[i, 3] <- modtot0610[[i]]$coefficients[2]-(1.96*est$se[2])
  tot0610[i, 4] <- modtot0610[[i]]$coefficients[2]+(1.96*est$se[2])
  pred0610[[i]]<-crosspred(pol,modtot0610[[i]])
  y0610[i,] <- pred0610[[i]]$coef
  S0610[[i]] <- pred0610[[i]]$vcov
}

# L1
for(i in 1:length(villes)) {
  poll_0610_L1[[i]] <- onebasis(y2006_y2010[[i]]$o3_L1,"lin")
  pol_L1<-poll_0610_L1[[i]]
  modtot0610_L1[[i]]<-gam(nocc_tot~pol_L1+ns(time,df=round(8*length(time)/365))+ns(tempmoy,df=6)+ns(temp3,df=6)+Jours+Vacances,data=y2006_y2010[[i]],family=quasipoisson)
  est_L1<-summary(modtot0610_L1[[i]])
  st_err0610_L1[i]<-est_L1$se[2]
  tot0610[i, 5] <- modtot0610_L1[[i]]$coefficients[2]
  tot0610[i, 6] <- modtot0610_L1[[i]]$coefficients[2]-(1.96*est_L1$se[2])
  tot0610[i, 7] <- modtot0610_L1[[i]]$coefficients[2]+(1.96*est_L1$se[2])
  pred0610_L1[[i]]<-crosspred(pol_L1,modtot0610_L1[[i]])
  y0610_L1[i,] <- pred0610_L1[[i]]$coef
  S0610_L1[[i]] <- pred0610_L1[[i]]$vcov
}

for(i in 1:length(villes)) {
  poll_0610_L2[[i]] <- onebasis(y2006_y2010[[i]]$o3_L2,"lin")
  pol_L2<-poll_0610_L2[[i]]
  modtot0610_L2[[i]]<-gam(nocc_tot~pol_L2+ns(time,df=round(8*length(time)/365))+ns(tempmoy,df=6)+ns(temp3,df=6)+Jours+Vacances,data=y2006_y2010[[i]],family=quasipoisson)
  est_L2<-summary(modtot0610_L2[[i]])
  st_err0610_L2[i]<-est_L2$se[2]
  tot0610[i, 8] <- modtot0610_L2[[i]]$coefficients[2]
  tot0610[i, 9] <- modtot0610_L2[[i]]$coefficients[2]-(1.96*est_L2$se[2])
  tot0610[i, 10] <- modtot0610_L2[[i]]$coefficients[2]+(1.96*est_L2$se[2])
  pred0610_L2[[i]]<-crosspred(pol_L2,modtot0610_L2[[i]])
  y0610_L2[i,] <- pred0610_L2[[i]]$coef
  S0610_L2[[i]] <- pred0610_L2[[i]]$vcov
}

### Third Period
# L0
for(i in 1:length(villes)) {
  poll_1115[[i]] <- onebasis(y2011_y2015[[i]]$o3.x,"lin")
  pol<-poll_1115[[i]]
  modtot1115[[i]]<-gam(nocc_tot~pol+ns(time,df=round(8*length(time)/365))+ns(tempmoy,df=6)+ns(temp3,df=6)+Jours+Vacances,data=y2011_y2015[[i]],family=quasipoisson)
  est<-summary(modtot1115[[i]])
  st_err1115[i]<-est$se[2]
  tot1115[i, 1] <- as.character(villes_name[[i]])
  tot1115[i, 2] <- modtot1115[[i]]$coefficients[2]
  tot1115[i, 3] <- modtot1115[[i]]$coefficients[2]-(1.96*est$se[2])
  tot1115[i, 4] <- modtot1115[[i]]$coefficients[2]+(1.96*est$se[2])
  pred1115[[i]]<-crosspred(pol,modtot1115[[i]])
  y1115[i,] <- pred1115[[i]]$coef
  S1115[[i]] <- pred1115[[i]]$vcov
}

# L1
for(i in 1:length(villes)) {
  poll_1115_L1[[i]] <- onebasis(y2011_y2015[[i]]$o3_L1,"lin")
  pol_L1<-poll_1115_L1[[i]]
  modtot1115_L1[[i]]<-gam(nocc_tot~pol_L1+ns(time,df=round(8*length(time)/365))+ns(tempmoy,df=6)+ns(temp3,df=6)+Jours+Vacances,data=y2011_y2015[[i]],family=quasipoisson)
  est_L1<-summary(modtot1115_L1[[i]])
  st_err1115_L1[i]<-est_L1$se[2]
  tot1115[i, 5] <- modtot1115_L1[[i]]$coefficients[2]
  tot1115[i, 6] <- modtot1115_L1[[i]]$coefficients[2]-(1.96*est_L1$se[2])
  tot1115[i, 7] <- modtot1115_L1[[i]]$coefficients[2]+(1.96*est_L1$se[2])
  pred1115_L1[[i]]<-crosspred(pol_L1,modtot1115_L1[[i]])
  y1115_L1[i,] <- pred1115_L1[[i]]$coef
  S1115_L1[[i]] <- pred1115_L1[[i]]$vcov
}


# L2
for(i in 1:length(villes)) {
  poll_1115_L2[[i]] <- onebasis(y2011_y2015[[i]]$o3_L2,"lin")
  pol_L2<-poll_1115_L2[[i]]
  modtot1115_L2[[i]]<-gam(nocc_tot~pol_L2+ns(time,df=round(8*length(time)/365))+ns(tempmoy,df=6)+ns(temp3,df=6)+Jours+Vacances,data=y2011_y2015[[i]],family=quasipoisson)
  est_L2<-summary(modtot1115_L2[[i]])
  st_err1115_L2[i]<-est_L2$se[2]
  tot1115[i, 8] <- modtot1115_L2[[i]]$coefficients[2]
  tot1115[i, 9] <- modtot1115_L2[[i]]$coefficients[2]-(1.96*est_L2$se[2])
  tot1115[i, 10] <- modtot1115_L2[[i]]$coefficients[2]+(1.96*est_L2$se[2])
  pred1115_L2[[i]]<-crosspred(pol_L2,modtot1115_L2[[i]])
  y1115_L2[i,] <- pred1115_L2[[i]]$coef
  S1115_L2[[i]] <- pred1115_L2[[i]]$vcov
}

Meta-Analysis results

For Lag0

First Period:

resL0_1 <- data.frame("ville"=G1_names,"coef" = as.numeric(as.character(tot0105[, 2])),"St_err" = st_err0105)

meta_L0_1<-rma(coef,sei=St_err,data=resL0_1)

#Forestplot
forest(meta_L0_1, slab = resL0_1$ville, refline=1,xlab="Lag 0",transf = exp, digits=4L)

meta_0105 <- mvmeta(y0105,S0105,method="ml")
summary(meta_0105)

## Call:  mvmeta(formula = y0105 ~ 1, S = S0105, method = "ml")
## 
## Univariate random-effects meta-analysis
## Dimension: 1
## Estimation method: ML
## 
## Fixed-effects coefficients
##              Estimate  Std. Error       z  Pr(>|z|)  95%ci.lb  95%ci.ub   
## (Intercept)    0.0001      0.0001  0.5971    0.5505   -0.0001    0.0003   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Between-study random-effects (co)variance components
##   Std. Dev
##     0.0000
## 
## Univariate Cochran Q-test for heterogeneity:
## Q = 8.2349 (df = 12), p-value = 0.7665
## I-square statistic = 1.0%
## 
## 13 studies, 13 observations, 1 fixed and 1 random-effects parameters
##    logLik        AIC        BIC  
##   83.7379  -163.4759  -162.3460

Second Period:

resL0_2 <- data.frame("ville"=villes_name,"coef" = as.numeric(as.character(tot0610[, 2])),"St_err" = st_err0610)

meta_L0_2<-rma(coef,sei=St_err,data=resL0_2)

#Forestplot
forest(meta_L0_2, slab = resL0_2$ville, refline=1,xlab="Lag 0",transf = exp, digits=4L)

meta_0610 <- mvmeta(y0610,S0610,method="ml")
summary(meta_0610)

## Call:  mvmeta(formula = y0610 ~ 1, S = S0610, method = "ml")
## 
## Univariate random-effects meta-analysis
## Dimension: 1
## Estimation method: ML
## 
## Fixed-effects coefficients
##              Estimate  Std. Error       z  Pr(>|z|)  95%ci.lb  95%ci.ub   
## (Intercept)    0.0002      0.0001  1.6024    0.1091   -0.0000    0.0004   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Between-study random-effects (co)variance components
##   Std. Dev
##     0.0000
## 
## Univariate Cochran Q-test for heterogeneity:
## Q = 13.9945 (df = 15), p-value = 0.5259
## I-square statistic = 1.0%
## 
## 16 studies, 16 observations, 1 fixed and 1 random-effects parameters
##    logLik        AIC        BIC  
##  100.4299  -196.8597  -195.3146

Third Period:

resL0_3 <- data.frame("ville"=villes_name,"coef" = as.numeric(as.character(tot1115[, 2])),"St_err" = st_err1115)

meta_L0_3<-rma(coef,sei=St_err,data=resL0_3)

#Forestplot
forest(meta_L0_3, slab = resL0_3$ville, refline=1,xlab="Lag 0",transf = exp, digits=4L)

meta_1115 <- mvmeta(y1115,S1115,method="ml")
summary(meta_1115)

## Call:  mvmeta(formula = y1115 ~ 1, S = S1115, method = "ml")
## 
## Univariate random-effects meta-analysis
## Dimension: 1
## Estimation method: ML
## 
## Fixed-effects coefficients
##              Estimate  Std. Error       z  Pr(>|z|)  95%ci.lb  95%ci.ub    
## (Intercept)    0.0003      0.0001  2.9355    0.0033    0.0001    0.0005  **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Between-study random-effects (co)variance components
##   Std. Dev
##     0.0000
## 
## Univariate Cochran Q-test for heterogeneity:
## Q = 6.7365 (df = 15), p-value = 0.9646
## I-square statistic = 1.0%
## 
## 16 studies, 16 observations, 1 fixed and 1 random-effects parameters
##    logLik        AIC        BIC  
##  104.6798  -205.3597  -203.8145

Results for Lag1

First Period:

resL1_1 <- data.frame("ville"=G1_names,"coef" = as.numeric(as.character(tot0105[, 5])),"St_err" = st_err0105_L1)

meta_L1_1<-rma(coef,sei=St_err,data=resL1_1)

#Forestplot
forest(meta_L1_1, slab = resL1_1$ville, refline=1,xlab="Lag 1", transf = exp, digits=4L)

meta_0105_L1 <- mvmeta(y0105_L1,S0105_L1,method="ml")
summary(meta_0105_L1)

## Call:  mvmeta(formula = y0105_L1 ~ 1, S = S0105_L1, method = "ml")
## 
## Univariate random-effects meta-analysis
## Dimension: 1
## Estimation method: ML
## 
## Fixed-effects coefficients
##              Estimate  Std. Error       z  Pr(>|z|)  95%ci.lb  95%ci.ub   
## (Intercept)    0.0002      0.0001  2.1487    0.0317    0.0000    0.0004  *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Between-study random-effects (co)variance components
##   Std. Dev
##     0.0000
## 
## Univariate Cochran Q-test for heterogeneity:
## Q = 9.0010 (df = 12), p-value = 0.7028
## I-square statistic = 1.0%
## 
## 13 studies, 13 observations, 1 fixed and 1 random-effects parameters
##    logLik        AIC        BIC  
##   83.4108  -162.8216  -161.6917

Second Period:

resL1_2 <- data.frame("ville"=villes_name,"coef" = as.numeric(as.character(tot0610[, 5])),"St_err" = st_err0610_L1)

meta_L1_2<-rma(coef,sei=St_err,data=resL1_2)

#Forestplot
forest(meta_L1_2, slab = resL1_2$ville, refline=1,xlab="Lag 1",transf = exp, digits=4L)

meta_0610_L1 <- mvmeta(y0610_L1,S0610_L1,method="ml")
summary(meta_0610_L1)

## Call:  mvmeta(formula = y0610_L1 ~ 1, S = S0610_L1, method = "ml")
## 
## Univariate random-effects meta-analysis
## Dimension: 1
## Estimation method: ML
## 
## Fixed-effects coefficients
##              Estimate  Std. Error       z  Pr(>|z|)  95%ci.lb  95%ci.ub    
## (Intercept)    0.0003      0.0001  2.8342    0.0046    0.0001    0.0005  **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Between-study random-effects (co)variance components
##   Std. Dev
##     0.0000
## 
## Univariate Cochran Q-test for heterogeneity:
## Q = 23.6661 (df = 15), p-value = 0.0710
## I-square statistic = 36.6%
## 
## 16 studies, 16 observations, 1 fixed and 1 random-effects parameters
##    logLik        AIC        BIC  
##   95.8630  -187.7260  -186.1808

Third Period:

resL1_3 <- data.frame("ville"=villes_name,"coef" = as.numeric(as.character(tot1115[, 5])),"St_err" = st_err1115_L1)

meta_L1_3<-rma(coef,sei=St_err,data=resL1_3)

#Forestplot
forest(meta_L1_3, slab = resL1_3$ville, refline=1,xlab="Lag 1",transf = exp, digits=4L)

meta_1115_L1 <- mvmeta(y1115_L1,S1115_L1,method="ml")
summary(meta_1115_L1)

## Call:  mvmeta(formula = y1115_L1 ~ 1, S = S1115_L1, method = "ml")
## 
## Univariate random-effects meta-analysis
## Dimension: 1
## Estimation method: ML
## 
## Fixed-effects coefficients
##              Estimate  Std. Error       z  Pr(>|z|)  95%ci.lb  95%ci.ub     
## (Intercept)    0.0004      0.0001  4.1511    0.0000    0.0002    0.0006  ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Between-study random-effects (co)variance components
##   Std. Dev
##     0.0000
## 
## Univariate Cochran Q-test for heterogeneity:
## Q = 9.2588 (df = 15), p-value = 0.8636
## I-square statistic = 1.0%
## 
## 16 studies, 16 observations, 1 fixed and 1 random-effects parameters
##    logLik        AIC        BIC  
##  103.5027  -203.0055  -201.4603

Results for Lag2

First Period:

resL2_1 <- data.frame("ville"=G1_names,"coef" = as.numeric(as.character(tot0105[, 8])),"St_err" = st_err0105_L2)

meta_L2_1<-rma(coef,sei=St_err,data=resL2_1)

#Forestplot
forest(meta_L2_1, slab = resL1_1$ville, refline=1,xlab="Lag 2", transf = exp, digits=4L)

meta_0105_L2 <- mvmeta(y0105_L2,S0105_L2,method="ml")
summary(meta_0105_L2)

## Call:  mvmeta(formula = y0105_L2 ~ 1, S = S0105_L2, method = "ml")
## 
## Univariate random-effects meta-analysis
## Dimension: 1
## Estimation method: ML
## 
## Fixed-effects coefficients
##              Estimate  Std. Error       z  Pr(>|z|)  95%ci.lb  95%ci.ub   
## (Intercept)    0.0002      0.0001  1.8595    0.0630   -0.0000    0.0004  .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Between-study random-effects (co)variance components
##   Std. Dev
##     0.0000
## 
## Univariate Cochran Q-test for heterogeneity:
## Q = 11.4551 (df = 12), p-value = 0.4904
## I-square statistic = 1.0%
## 
## 13 studies, 13 observations, 1 fixed and 1 random-effects parameters
##    logLik        AIC        BIC  
##   82.3687  -160.7374  -159.6076

Second Period:

resL2_2 <- data.frame("ville"=villes_name,"coef" = as.numeric(as.character(tot0610[, 8])),"St_err" = st_err0610_L2)

meta_L2_2<-rma(coef,sei=St_err,data=resL2_2)

#Forestplot
forest(meta_L2_2, slab = resL2_2$ville, refline=1,xlab="Lag 2",transf = exp, digits=4L)

meta_0610_L2 <- mvmeta(y0610_L2,S0610_L2,method="ml")
summary(meta_0610_L2)

## Call:  mvmeta(formula = y0610_L2 ~ 1, S = S0610_L2, method = "ml")
## 
## Univariate random-effects meta-analysis
## Dimension: 1
## Estimation method: ML
## 
## Fixed-effects coefficients
##              Estimate  Std. Error       z  Pr(>|z|)  95%ci.lb  95%ci.ub   
## (Intercept)    0.0002      0.0001  1.3253    0.1851   -0.0001    0.0004   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Between-study random-effects (co)variance components
##   Std. Dev
##     0.0002
## 
## Univariate Cochran Q-test for heterogeneity:
## Q = 19.9531 (df = 15), p-value = 0.1737
## I-square statistic = 24.8%
## 
## 16 studies, 16 observations, 1 fixed and 1 random-effects parameters
##    logLik        AIC        BIC  
##   97.8582  -191.7164  -190.1712

Third Period:

resL2_3 <- data.frame("ville"=villes_name,"coef" = as.numeric(as.character(tot1115[, 8])),"St_err" = st_err1115_L2)

meta_L2_3<-rma(coef,sei=St_err,data=resL1_3)

#Forestplot
forest(meta_L2_3, slab = resL2_3$ville, refline=1,xlab="Lag 2",transf = exp, digits=4L)

meta_1115_L2 <- mvmeta(y1115_L2,S1115_L2,method="ml")
summary(meta_1115_L2)

## Call:  mvmeta(formula = y1115_L2 ~ 1, S = S1115_L2, method = "ml")
## 
## Univariate random-effects meta-analysis
## Dimension: 1
## Estimation method: ML
## 
## Fixed-effects coefficients
##              Estimate  Std. Error       z  Pr(>|z|)  95%ci.lb  95%ci.ub     
## (Intercept)    0.0004      0.0001  4.2984    0.0000    0.0002    0.0006  ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Between-study random-effects (co)variance components
##   Std. Dev
##     0.0000
## 
## Univariate Cochran Q-test for heterogeneity:
## Q = 10.9673 (df = 15), p-value = 0.7549
## I-square statistic = 1.0%
## 
## 16 studies, 16 observations, 1 fixed and 1 random-effects parameters
##    logLik        AIC        BIC  
##  102.7181  -201.4362  -199.8910

Table of results: change in mortality risk per 10 μg/m3 increase in PM10 for 3 time periods, with Temporal Change expressed in percentage points

####Results for Lag0

L0<-matrix(NA,nrow=16,ncol=12)
colnames(L0)<-c("villes","Period1","P1_CI95_Low","P1_CI95_High","Period2","P2_CI95_Low","P2_CI95_High","Period3","P3_CI95_Low","P3_CI95_High","Temp_Change1_2","Temp_Change2_3")
L0<-as.data.frame(L0)
L0[,1]<-names(villes)


L0[which(L0$villes %in% names(villesG1)),]$Period1<-exp(tot0105$Coefficient_L0*10)
L0$Period2<-exp(tot0610$Coefficient_L0*10)
L0$Period3<-exp(tot1115$Coefficient_L0*10)

L0[which(L0$villes %in% names(villesG1)),]$P1_CI95_Low<-exp(tot0105$L0_CI95_Low*10)
L0$P2_CI95_Low<-exp(tot0610$L0_CI95_Low*10)
L0$P3_CI95_Low<-exp(tot1115$L0_CI95_Low*10)

L0[which(L0$villes %in% names(villesG1)),]$P1_CI95_High<-exp(tot0105$L0_CI95_High*10)
L0$P2_CI95_High<-exp(tot0610$L0_CI95_High*10)
L0$P3_CI95_High<-exp(tot1115$L0_CI95_High*10)


L0[which(L0$villes %in% names(villesG1)),]$Temp_Change1_2<-(exp(tot0610[which(tot0610$Ville %in% names(villesG1)),]$Coefficient_L0*10)-exp(tot0105$Coefficient_L0*10))*100

L0$Temp_Change2_3<-(exp(tot1115$Coefficient_L0*10)-exp(tot0610$Coefficient_L0*10))*100


L0[,2:12]<-round(L0[,2:12],digits = 3)


kable(L0) %>% kable_styling(bootstrap_options = c("hover", "condensed"),full_width = F, position = "center")%>% add_header_above(c(" ", "Period 1" = 3, "Period 2" = 3, "Period 3"=3, "Temporal Change (%)"=2))%>%column_spec(1, bold = T, border_right = T)%>%
  column_spec(4, border_right = T)%>%
  column_spec(7, border_right = T)%>%
  column_spec(10, border_right = T)

	Period 1			Period 2			Period 3			Temporal Change (%)
villes	Period1	P1_CI95_Low	P1_CI95_High	Period2	P2_CI95_Low	P2_CI95_High	Period3	P3_CI95_Low	P3_CI95_High	Temp_Change1_2	Temp_Change2_3
bordeaux	0.997	0.988	1.007	1.004	0.994	1.014	0.999	0.991	1.008	0.700	-0.501
clermont	1.005	0.992	1.017	1.000	0.987	1.013	1.001	0.989	1.013	-0.501	0.100
grenoble	1.002	0.992	1.011	1.006	0.996	1.016	0.999	0.990	1.009	0.402	-0.702
lehavre	0.996	0.983	1.009	1.000	0.987	1.014	1.005	0.992	1.018	0.399	0.501
lille	1.000	0.994	1.006	1.002	0.995	1.009	1.004	0.997	1.011	0.200	0.201
lyon	NA	NA	NA	1.002	0.995	1.009	1.003	0.997	1.010	NA	0.100
marseille	1.006	1.000	1.013	1.004	0.998	1.011	1.005	0.999	1.012	-0.201	0.100
montpellier	1.005	0.994	1.016	1.004	0.992	1.017	0.999	0.987	1.012	-0.100	-0.501
nancy	NA	NA	NA	0.994	0.983	1.006	1.000	0.989	1.011	NA	0.598
nantes	0.994	0.984	1.004	0.996	0.986	1.007	1.004	0.995	1.014	0.199	0.800
nice	NA	NA	NA	1.006	0.997	1.016	0.994	0.984	1.004	NA	-1.200
paris	1.000	0.997	1.003	1.002	0.999	1.006	1.003	1.000	1.006	0.200	0.100
rennes	1.002	0.986	1.018	0.998	0.981	1.015	1.004	0.989	1.020	-0.400	0.601
rouen	1.003	0.993	1.013	1.005	0.995	1.016	1.005	0.996	1.015	0.201	0.000
strasbourg	0.996	0.987	1.005	0.996	0.987	1.005	1.004	0.995	1.013	0.000	0.800
toulouse	1.004	0.994	1.014	0.988	0.978	0.999	1.008	0.998	1.018	-1.594	1.996

#write_xlsx(L0,"C:/Users/Anna/Dropbox/Leo/change DR curves/results/Ozone/L0_periods.xlsx")


resL0 <- data.frame("Period1" = as.numeric(as.character(meta_L0_1$b)),"Period2" = as.numeric(as.character(meta_L0_2$b)),"Period3"=as.numeric(as.character(meta_L0_3$b)))

resL0$change12<-(exp(resL0$Period2*10)-exp(resL0$Period1*10))*100
resL0$change23<-(exp(resL0$Period3*10)-exp(resL0$Period2*10))*100

resL0[1:3]<-exp(resL0[1:3])
kable(resL0)%>% kable_styling(bootstrap_options = c("hover", "condensed"),full_width = F, position = "center")

Period1	Period2	Period3	change12	change23
1.000061	1.000162	1.000289	0.1003545	0.1270233

####Results for Lag1

	Period 1			Period 2			Period 3			Temporal Change (%)
villes	Period1	P1_CI95_Low	P1_CI95_High	Period2	P2_CI95_Low	P2_CI95_High	Period3	P3_CI95_Low	P3_CI95_High	Temp_Change1_2	Temp_Change2_3
bordeaux	1.001	0.992	1.010	1.009	1.000	1.019	1.003	0.995	1.012	0.804	-0.604
clermont	1.005	0.992	1.018	0.991	0.978	1.004	1.005	0.993	1.017	-1.397	1.397
grenoble	0.996	0.986	1.006	0.998	0.988	1.008	0.999	0.989	1.009	0.199	0.100
lehavre	1.000	0.987	1.013	0.997	0.984	1.010	1.013	1.000	1.026	-0.300	1.608
lille	1.004	0.998	1.011	1.001	0.995	1.008	1.002	0.996	1.009	-0.301	0.100
lyon	NA	NA	NA	1.006	0.999	1.013	1.003	0.997	1.010	NA	-0.301
marseille	1.002	0.995	1.008	1.005	0.998	1.011	1.000	0.993	1.007	0.301	-0.501
montpellier	1.007	0.995	1.018	1.005	0.993	1.018	1.010	0.998	1.022	-0.201	0.504
nancy	NA	NA	NA	0.999	0.988	1.011	1.004	0.993	1.015	NA	0.501
nantes	0.999	0.989	1.009	1.011	1.000	1.021	1.008	0.998	1.017	1.206	-0.303
nice	NA	NA	NA	1.011	1.001	1.020	1.004	0.994	1.015	NA	-0.705
paris	1.002	0.999	1.005	1.003	1.000	1.006	1.004	1.001	1.007	0.100	0.100
rennes	0.995	0.979	1.011	0.996	0.979	1.013	0.993	0.978	1.009	0.100	-0.298
rouen	1.003	0.993	1.013	1.009	0.999	1.019	1.007	0.998	1.017	0.604	-0.202
strasbourg	1.002	0.994	1.011	0.991	0.982	1.000	1.007	0.999	1.016	-1.096	1.598
toulouse	1.013	1.003	1.023	0.997	0.987	1.007	1.007	0.997	1.017	-1.608	1.002

Period1	Period2	Period3	change12	change23
1.000219	1.00025	1.000404	0.0311357	0.1539444

####Results for Lag2

	Period 1			Period 2			Period 3			Temporal Change (%)
villes	Period1	P1_CI95_Low	P1_CI95_High	Period2	P2_CI95_Low	P2_CI95_High	Period3	P3_CI95_Low	P3_CI95_High	Temp_Change1_2	Temp_Change2_3
bordeaux	1.000	0.991	1.009	1.005	0.995	1.014	1.005	0.996	1.013	0.501	0.000
clermont	1.012	0.999	1.025	0.999	0.986	1.011	1.006	0.994	1.018	-1.307	0.702
grenoble	0.999	0.990	1.009	0.996	0.986	1.007	1.007	0.997	1.017	-0.299	1.102
lehavre	1.001	0.989	1.013	0.999	0.986	1.012	1.008	0.996	1.021	-0.200	0.903
lille	1.003	0.996	1.009	0.999	0.992	1.005	1.002	0.995	1.008	-0.400	0.300
lyon	NA	NA	NA	1.001	0.994	1.007	1.006	0.999	1.013	NA	0.502
marseille	1.006	0.999	1.012	1.006	1.000	1.013	0.999	0.992	1.006	0.000	-0.702
montpellier	1.013	1.002	1.024	0.994	0.981	1.006	1.008	0.996	1.020	-1.907	1.401
nancy	NA	NA	NA	1.010	0.998	1.021	1.012	1.001	1.023	NA	0.202
nantes	1.001	0.991	1.010	1.009	0.999	1.020	1.005	0.996	1.014	0.804	-0.403
nice	NA	NA	NA	1.007	0.997	1.017	1.005	0.994	1.015	NA	-0.201
paris	1.000	0.997	1.003	0.999	0.996	1.003	1.003	1.000	1.006	-0.100	0.400
rennes	0.992	0.976	1.007	0.992	0.976	1.009	0.992	0.977	1.007	0.000	0.000
rouen	0.999	0.990	1.009	1.010	1.000	1.020	1.006	0.996	1.015	1.105	-0.403
strasbourg	1.003	0.994	1.012	0.994	0.985	1.003	1.007	0.998	1.016	-0.899	1.301
toulouse	1.006	0.996	1.016	1.005	0.995	1.015	1.010	1.000	1.020	-0.101	0.504

Period1	Period2	Period3	change12	change23
1.000213	1.000165	1.000404	-0.0475805	0.2390857

Metanalysis of the results

Pooled Dose-Responses Curves

L0

pol0105 <- seq(min(bound_o3_0105[,1]),max(bound_o3_0105[,2]),length=220)
pol0105 <- onebasis(pol0105,"lin")

pol0610 <- seq(min(bound_o3_0610[,1]),max(bound_o3_0610[,2]),length=220)
pol0610 <- onebasis(pol0610,"lin")

pol1115 <- seq(min(bound_o3_1115[,1]),max(bound_o3_1115[,2]),length=220)
pol1115 <- onebasis(pol1115,"lin")


plot_0105 <- crosspred(pol0105,coef=coef(meta_0105),vcov=vcov(meta_0105),model.link="log",by=1)
plot(plot_0105,"overall",ci="lines",col=1,lwd=3,ylab="RR",xlab="PM10",main="2001-2005")

plot_0610 <- crosspred(pol0610,coef=coef(meta_0610),vcov=vcov(meta_0610),model.link="log",by=1)
plot(plot_0610,"overall",ci="lines",col=1,lwd=3,ylab="RR",xlab="PM10",main="2006-2010")

plot_1115 <- crosspred(pol1115,coef=coef(meta_1115),vcov=vcov(meta_1115),model.link="log",by=1)
plot(plot_1115,"overall",ci="lines",col=1,lwd=3,ylab="RR",xlab="PM10",main="2011-2015")

L1

plot_0105_L1 <- crosspred(pol0105,coef=coef(meta_0105_L1),vcov=vcov(meta_0105_L1),model.link="log",by=1)
plot(plot_0105_L1,"overall",ci="lines",col=1,lwd=3,ylab="RR",xlab="PM10_L1",main="2001-2005")

plot_0610_L1 <- crosspred(pol0610,coef=coef(meta_0610_L1),vcov=vcov(meta_0610_L1),model.link="log",by=1)
plot(plot_0610_L1,"overall",ci="lines",col=1,lwd=3,ylab="RR",xlab="PM10_L1",main="2006-2010")

plot_1115_L1 <- crosspred(pol1115,coef=coef(meta_1115_L1),vcov=vcov(meta_1115_L1),model.link="log",by=1)
plot(plot_1115_L1,"overall",ci="lines",col=1,lwd=3,ylab="RR",xlab="PM10_L1",main="2011-2015")

L2

plot_0105_L2 <- crosspred(pol0105,coef=coef(meta_0105_L2),vcov=vcov(meta_0105_L2),model.link="log",by=1)
plot(plot_0105_L2,"overall",ci="lines",col=1,lwd=3,ylab="RR",xlab="PM10_L2",main="2001-2005")

plot_0610_L2 <- crosspred(pol0610,coef=coef(meta_0610_L2),vcov=vcov(meta_0610_L2),model.link="log",by=1)
plot(plot_0610_L2,"overall",ci="lines",col=1,lwd=3,ylab="RR",xlab="PM10_L2",main="2006-2010")

plot_1115_L2 <- crosspred(pol1115,coef=coef(meta_1115_L2),vcov=vcov(meta_1115_L2),model.link="log",by=1)
plot(plot_1115_L2,"overall",ci="lines",col=1,lwd=3,ylab="RR",xlab="PM10",main="2011-2015")

Doses-Response Curves by cities

Bordeaux

Lag 0

## Plot per city, exemple
i=1
par (mfrow=c(2,2),mar=c(1,1,3,1))
plot(pred0105[[i]],ci="lines",col=1,lwd=3,ylab="RR",xlab="o3",main="2001-2005")
plot(pred0610[[i]],ci="lines",col=1,lwd=3,ylab="RR",xlab="o3",main="2006-2010")
plot(pred1115[[i]],ci="lines",col=1,lwd=3,ylab="RR",xlab="o3",main="2011-2015")

Lag1

Lag2

Clermont-Ferrand

Lag 0

Lag 1

Lag 2

Grenoble

Lag 0

Lag 1

Lag 2

Le Havre

Lag 0

Lag 1

Lag 2

Lille

Lag 0

Lag 1

Lag 2

Lyon

Lag 0

Lag 1

Lag 2

Marseille

Lag 0

Lag 1

Lag 2

Montpellier

Lag 0

Lag 1

Lag 2

Nancy

Lag 0

Lag 1

Lag 2

Nantes

Lag 0

Lag 1

Lag 2

Nice

Lag 0

Lag 1

Lag 2

Paris

Lag 0

Lag 1

Lag 2

Rennes

Lag 0

Lag 1

Lag 2

Rouen

Lag 0

Lag 1

Lag 2

Strasbourg

Lag 0

Lag 1

Lag 2

Toulouse

Lag 0

Lag 1

Lag 2