D-R Curves for Ozone and mortality over the time with two-stage modeling - Linear and Non linear interaction Models
Anna Alari
28 février, 2021
Descriptive Statistics for O3
| Ville | Minimum | X1st.Qu | Median | Mean | X3rd.Qu | Max | NA.s |
|---|---|---|---|---|---|---|---|
| bordeaux | 1.9154062 | 51.33333 | 70.95833 | 70.70046 | 89.16369 | 169.5000 | NA |
| clermont | 1.7321429 | 54.12500 | 72.68750 | 72.06568 | 90.12500 | 178.5625 | 138 |
| grenoble | 0.4375000 | 35.51116 | 65.96875 | 65.82133 | 91.66518 | 185.8750 | NA |
| lehavre | 2.1250000 | 57.25000 | 70.62500 | 69.51529 | 82.12500 | 184.5000 | 59 |
| lille | 0.1666667 | 41.00000 | 59.00000 | 58.91898 | 74.75000 | 207.1250 | 239 |
| marseille | 3.1250000 | 54.75000 | 79.59742 | 77.44850 | 99.62500 | 185.6250 | 16 |
| nantes | 3.1250000 | 55.00000 | 70.37500 | 71.71581 | 87.62500 | 216.5000 | 71 |
| paris | 0.8592206 | 37.65048 | 57.72942 | 59.16113 | 77.16338 | 216.3857 | NA |
| rennes | 1.5625000 | 52.12500 | 65.93750 | 66.48971 | 80.25005 | 213.5625 | 17 |
| rouen | 1.2972346 | 46.96354 | 63.34375 | 63.70247 | 79.04688 | 178.6607 | NA |
| strasbourg | 0.0000000 | 36.75000 | 62.87500 | 64.30176 | 87.75000 | 206.5000 | 99 |
| toulouse | 1.6250000 | 57.31250 | 76.27500 | 76.08195 | 94.61250 | 180.9594 | NA |
Linear Interaction Model
Models For Lag0, Lag1 and La2 exposure
## change year variable with numeric ordinal number
for (i in 1:length(villes)){
villes[[i]]$year<-as.numeric(villes[[i]]$annee)
}
# define model objects
tot <- matrix(NA, nrow = 12, ncol = 7)
colnames(tot) <- c("Ville","Coefficient_B0","B0_CI95_Low","B0_CI95_High", "Coefficient_B1","LB_CI95_Low","B1_CI95_High")
modtot<-list()
st_errb0<-c()
st_errb1<-c()
# Linear Model for Lag0
for(i in 1:length(villes)) {
modtot[[i]]<-gam(nocc_tot~o3.x+o3.x:year+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]])
ll<-length(modtot[[i]]$coefficients)
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])
tot[i, 5] <- modtot[[i]]$coefficients[ll]
tot[i, 6] <- modtot[[i]]$coefficients[ll]-(1.96*est$se[ll])
tot[i, 7] <- modtot[[i]]$coefficients[ll]+(1.96*est$se[ll])
st_errb0[i]<- est$se[2]
st_errb1[i]<- est$se[ll]
}
# Linear Model for Lag1
# define model objects
tot_L1 <- matrix(NA, nrow = 12, ncol = 7)
colnames(tot_L1) <- c("Ville","Coefficient_B0","B0_CI95_Low","B0_CI95_High", "Coefficient_B1","LB_CI95_Low","B1_CI95_High")
modtot_L1<-list()
st_errb0_L1<-c()
st_errb1_L1<-c()
for(i in 1:length(villes)) {
modtot_L1[[i]]<-gam(nocc_tot~o3_L1+o3_L1:year+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_L1[[i]])
tot_L1[i, 1] <- as.character(villes_name[[i]])
ll<-length(modtot_L1[[i]]$coefficients)
tot_L1[i, 2] <- modtot_L1[[i]]$coefficients[2]
tot_L1[i, 3] <- modtot_L1[[i]]$coefficients[2]-(1.96*est$se[2])
tot_L1[i, 4] <- modtot_L1[[i]]$coefficients[2]+(1.96*est$se[2])
tot_L1[i, 5] <- modtot_L1[[i]]$coefficients[ll]
tot_L1[i, 6] <- modtot_L1[[i]]$coefficients[ll]-(1.96*est$se[ll])
tot_L1[i, 7] <- modtot_L1[[i]]$coefficients[ll]+(1.96*est$se[ll])
st_errb0_L1[i]<- est$se[2]
st_errb1_L1[i]<- est$se[ll]
}## Warning in gam.fit3(x = args$X, y = args$y, sp = lsp, Eb = args$Eb, UrS =
## args$UrS, : Algorithm did not converge# Linear Model for Lag2
# define model objects
tot_L2 <- matrix(NA, nrow = 12, ncol = 7)
colnames(tot_L2) <- c("Ville","Coefficient_B0","B0_CI95_Low","B0_CI95_High", "Coefficient_B1","LB_CI95_Low","B1_CI95_High")
modtot_L2<-list()
st_errb0_L2<-c()
st_errb1_L2<-c()
for(i in 1:length(villes)) {
modtot_L2[[i]]<-gam(nocc_tot~o3_L2+o3_L2:year+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_L2[[i]])
tot_L2[i, 1] <- as.character(villes_name[[i]])
ll<-length(modtot_L2[[i]]$coefficients)
tot_L2[i, 2] <- modtot_L2[[i]]$coefficients[2]
tot_L2[i, 3] <- modtot_L2[[i]]$coefficients[2]-(1.96*est$se[2])
tot_L2[i, 4] <- modtot_L2[[i]]$coefficients[2]+(1.96*est$se[2])
tot_L2[i, 5] <- modtot_L2[[i]]$coefficients[ll]
tot_L2[i, 6] <- modtot_L2[[i]]$coefficients[ll]-(1.96*est$se[ll])
tot_L2[i, 7] <- modtot_L2[[i]]$coefficients[ll]+(1.96*est$se[ll])
st_errb0_L2[i]<- est$se[2]
st_errb1_L2[i]<- est$se[ll]
}## Warning in gam.fit3(x = args$X, y = args$y, sp = lsp, Eb = args$Eb, UrS =
## args$UrS, : Algorithm did not convergeResults For Lag0 exposure
tot<-data.frame(tot)
for (i in (2:7)){
tot[,i]<-as.numeric(as.character(tot[,i]))
}
results<-data.frame("villes"=villes_name,"Increase_2002"=100*(exp(tot$Coefficient_B0*10 + tot$Coefficient_B1*2*10)-exp(10*tot$Coefficient_B0)),"Increase_2015"=100*(exp(tot$Coefficient_B0*10 + tot$Coefficient_B1*15*10)-exp(10*tot$Coefficient_B0)))
results$Temp_change<-results$Increase_2015-results$Increase_2002
tot_exp<-tot
tot_exp[,2:7]<-exp(tot_exp[,2:7])
tot_exp[,2:7]<-round(tot_exp[,2:7],digits = 4)
kable(tot_exp, align = c("r","c","c","c","c","c","c"))%>%kable_classic_2(full_width = F)%>%
column_spec(1, bold = T, border_right = T)%>%
column_spec(4, border_right = T)| Ville | Coefficient_B0 | B0_CI95_Low | B0_CI95_High | Coefficient_B1 | LB_CI95_Low | B1_CI95_High |
|---|---|---|---|---|---|---|
| bordeaux | 0.9995 | 0.9989 | 1.0001 | 1.0001 | 1.0000 | 1.0001 |
| clermont | 0.9999 | 0.9991 | 1.0006 | 1.0001 | 1.0001 | 1.0001 |
| grenoble | 0.9996 | 0.9990 | 1.0002 | 1.0001 | 1.0001 | 1.0002 |
| lehavre | 0.9997 | 0.9989 | 1.0005 | 1.0001 | 1.0000 | 1.0001 |
| lille | 1.0003 | 0.9999 | 1.0007 | 1.0000 | 1.0000 | 1.0000 |
| marseille | 1.0008 | 1.0004 | 1.0012 | 1.0000 | 1.0000 | 1.0000 |
| nantes | 0.9993 | 0.9986 | 0.9999 | 1.0001 | 1.0001 | 1.0002 |
| paris | 1.0006 | 1.0004 | 1.0009 | 0.9999 | 0.9999 | 1.0000 |
| rennes | 0.9994 | 0.9984 | 1.0003 | 1.0001 | 1.0001 | 1.0002 |
| rouen | 1.0002 | 0.9996 | 1.0009 | 1.0000 | 1.0000 | 1.0001 |
| strasbourg | 0.9995 | 0.9990 | 1.0001 | 1.0001 | 1.0000 | 1.0001 |
| toulouse | 0.9998 | 0.9992 | 1.0004 | 1.0001 | 1.0000 | 1.0001 |
results[,2:4]<-round(results[,2:4],digits = 3)
kable(results, align = c("r","c","c","c"))%>%kable_styling(full_width = F)%>%
column_spec(1, bold = T, border_right = T)%>%
column_spec(3, border_right = T)%>%
column_spec(4, bold = T)| villes | Increase_2002 | Increase_2015 | Temp_change |
|---|---|---|---|
| bordeaux | 0.104 | 0.785 | 0.681 |
| clermont | 0.184 | 1.391 | 1.206 |
| grenoble | 0.272 | 2.057 | 1.785 |
| lehavre | 0.110 | 0.825 | 0.716 |
| lille | -0.032 | -0.242 | -0.210 |
| marseille | -0.049 | -0.368 | -0.319 |
| nantes | 0.248 | 1.873 | 1.626 |
| paris | -0.119 | -0.889 | -0.770 |
| rennes | 0.205 | 1.547 | 1.342 |
| rouen | 0.091 | 0.687 | 0.595 |
| strasbourg | 0.154 | 1.162 | 1.008 |
| toulouse | 0.123 | 0.930 | 0.806 |
Results For Lag1 exposure
tot_L1<-data.frame(tot_L1)
for (i in (2:7)){
tot_L1[,i]<-as.numeric(as.character(tot_L1[,i]))
}
results_L1<-data.frame("villes"=villes_name,"Increase_2002"=100*(exp(tot_L1$Coefficient_B0*10 + tot_L1$Coefficient_B1*2*10)-exp(10*tot_L1$Coefficient_B0)),"Increase_2015"=100*(exp(tot_L1$Coefficient_B0*10 + tot_L1$Coefficient_B1*15*10)-exp(10*tot_L1$Coefficient_B0)))
results_L1$Temp_change<-results_L1$Increase_2015-results_L1$Increase_2002
tot_exp_L1<-tot
tot_exp_L1[,2:7]<-exp(tot_exp_L1[,2:7])
tot_exp_L1[,2:7]<-round(tot_exp_L1[,2:7],digits = 4)
kable(tot_exp_L1, align = c("r","c","c","c","c","c","c"))%>%kable_classic_2(full_width = F)%>%
column_spec(1, bold = T, border_right = T)%>%
column_spec(4, border_right = T)| Ville | Coefficient_B0 | B0_CI95_Low | B0_CI95_High | Coefficient_B1 | LB_CI95_Low | B1_CI95_High |
|---|---|---|---|---|---|---|
| bordeaux | 0.9995 | 0.9989 | 1.0001 | 1.0001 | 1.0000 | 1.0001 |
| clermont | 0.9999 | 0.9991 | 1.0006 | 1.0001 | 1.0001 | 1.0001 |
| grenoble | 0.9996 | 0.9990 | 1.0002 | 1.0001 | 1.0001 | 1.0002 |
| lehavre | 0.9997 | 0.9989 | 1.0005 | 1.0001 | 1.0000 | 1.0001 |
| lille | 1.0003 | 0.9999 | 1.0007 | 1.0000 | 1.0000 | 1.0000 |
| marseille | 1.0008 | 1.0004 | 1.0012 | 1.0000 | 1.0000 | 1.0000 |
| nantes | 0.9993 | 0.9986 | 0.9999 | 1.0001 | 1.0001 | 1.0002 |
| paris | 1.0006 | 1.0004 | 1.0009 | 0.9999 | 0.9999 | 1.0000 |
| rennes | 0.9994 | 0.9984 | 1.0003 | 1.0001 | 1.0001 | 1.0002 |
| rouen | 1.0002 | 0.9996 | 1.0009 | 1.0000 | 1.0000 | 1.0001 |
| strasbourg | 0.9995 | 0.9990 | 1.0001 | 1.0001 | 1.0000 | 1.0001 |
| toulouse | 0.9998 | 0.9992 | 1.0004 | 1.0001 | 1.0000 | 1.0001 |
results_L1[,2:4]<-round(results_L1[,2:4],digits = 3)
kable(results_L1, align = c("r","c","c","c"))%>%kable_styling(full_width = F)%>%
column_spec(1, bold = T, border_right = T)%>%
column_spec(3, border_right = T)%>%
column_spec(4, bold = T)| villes | Increase_2002 | Increase_2015 | Temp_change |
|---|---|---|---|
| bordeaux | 0.101 | 0.760 | 0.659 |
| clermont | 0.188 | 1.417 | 1.229 |
| grenoble | 0.272 | 2.059 | 1.787 |
| lehavre | 0.111 | 0.839 | 0.728 |
| lille | -0.036 | -0.266 | -0.230 |
| marseille | -0.052 | -0.391 | -0.338 |
| nantes | 0.235 | 1.779 | 1.544 |
| paris | -0.122 | -0.909 | -0.787 |
| rennes | 0.194 | 1.461 | 1.267 |
| rouen | 0.096 | 0.719 | 0.623 |
| strasbourg | 0.145 | 1.093 | 0.948 |
| toulouse | 0.121 | 0.914 | 0.792 |
Results For Lag2 exposure
tot_L2<-data.frame(tot_L2)
for (i in (2:7)){
tot_L2[,i]<-as.numeric(as.character(tot_L2[,i]))
}
results_L2<-data.frame("villes"=villes_name,"Increase_2002"=100*(exp(tot_L2$Coefficient_B0*10 + tot_L2$Coefficient_B1*2*10)-exp(10*tot_L2$Coefficient_B0)),"Increase_2015"=100*(exp(tot_L2$Coefficient_B0*10 + tot_L2$Coefficient_B1*15*10)-exp(10*tot_L2$Coefficient_B0)))
results_L2$Temp_change<-results_L2$Increase_2015-results_L2$Increase_2002
tot_exp_L2<-tot
tot_exp_L2[,2:7]<-exp(tot_exp_L2[,2:7])
tot_exp_L2[,2:7]<-round(tot_exp_L2[,2:7],digits = 4)
kable(tot_exp_L2, align = c("r","c","c","c","c","c","c"))%>%kable_classic_2(full_width = F)%>%
column_spec(1, bold = T, border_right = T)%>%
column_spec(4, border_right = T)| Ville | Coefficient_B0 | B0_CI95_Low | B0_CI95_High | Coefficient_B1 | LB_CI95_Low | B1_CI95_High |
|---|---|---|---|---|---|---|
| bordeaux | 0.9995 | 0.9989 | 1.0001 | 1.0001 | 1.0000 | 1.0001 |
| clermont | 0.9999 | 0.9991 | 1.0006 | 1.0001 | 1.0001 | 1.0001 |
| grenoble | 0.9996 | 0.9990 | 1.0002 | 1.0001 | 1.0001 | 1.0002 |
| lehavre | 0.9997 | 0.9989 | 1.0005 | 1.0001 | 1.0000 | 1.0001 |
| lille | 1.0003 | 0.9999 | 1.0007 | 1.0000 | 1.0000 | 1.0000 |
| marseille | 1.0008 | 1.0004 | 1.0012 | 1.0000 | 1.0000 | 1.0000 |
| nantes | 0.9993 | 0.9986 | 0.9999 | 1.0001 | 1.0001 | 1.0002 |
| paris | 1.0006 | 1.0004 | 1.0009 | 0.9999 | 0.9999 | 1.0000 |
| rennes | 0.9994 | 0.9984 | 1.0003 | 1.0001 | 1.0001 | 1.0002 |
| rouen | 1.0002 | 0.9996 | 1.0009 | 1.0000 | 1.0000 | 1.0001 |
| strasbourg | 0.9995 | 0.9990 | 1.0001 | 1.0001 | 1.0000 | 1.0001 |
| toulouse | 0.9998 | 0.9992 | 1.0004 | 1.0001 | 1.0000 | 1.0001 |
results_L2[,2:4]<-round(results_L2[,2:4],digits = 3)
kable(results_L2, align = c("r","c","c","c"))%>%kable_styling(full_width = F)%>%
column_spec(1, bold = T, border_right = T)%>%
column_spec(3, border_right = T)%>%
column_spec(4, bold = T)| villes | Increase_2002 | Increase_2015 | Temp_change |
|---|---|---|---|
| bordeaux | 0.104 | 0.784 | 0.680 |
| clermont | 0.188 | 1.419 | 1.231 |
| grenoble | 0.282 | 2.133 | 1.852 |
| lehavre | 0.111 | 0.837 | 0.726 |
| lille | -0.038 | -0.283 | -0.245 |
| marseille | -0.053 | -0.393 | -0.341 |
| nantes | 0.240 | 1.815 | 1.575 |
| paris | -0.123 | -0.917 | -0.794 |
| rennes | 0.203 | 1.534 | 1.331 |
| rouen | 0.101 | 0.761 | 0.660 |
| strasbourg | 0.148 | 1.112 | 0.965 |
| toulouse | 0.127 | 0.954 | 0.827 |
Metanalysis of the results
For Lag0 exposure
## For Lag0
resL0 <- data.frame("ville"=villes_name,"Coefficient_B0" = as.numeric(as.character(tot[, 2])),"St_err_B0" = st_errb0,"Coefficient_B1" = as.numeric(as.character(tot[, 5])),"St_err_B1" = st_errb1)
meta_L0_b0<-rma(Coefficient_B0,sei=St_err_B0,data=resL0[,1:3])
meta_L0_b1<-rma(Coefficient_B1,sei=St_err_B1,data=resL0)
results_meta<-data.frame("Increase_2002"=100*(exp(meta_L0_b0$b*10 + meta_L0_b1$b*2*10)-exp(10*meta_L0_b0$b)),"Increase_2015"=100*(exp(meta_L0_b0$b*10 + meta_L0_b1$b*15*10)-exp(10*meta_L0_b0$b)))
results_meta$Temp_change<-results_meta$Increase_2015-results_meta$Increase_2002
print(results_meta)## Increase_2002 Increase_2015 Temp_change
## intrcpt 0.1051638 0.7914312 0.6862673Results For Lag1 exposure
## For Lag1
resL1 <- data.frame("ville"=villes_name,"Coefficient_B0" = as.numeric(as.character(tot_L1[, 2])),"St_err_B0" = st_errb0_L1,"Coefficient_B1" = as.numeric(as.character(tot_L1[, 5])),"St_err_B1" = st_errb1_L1)
meta_L1_b0<-rma(Coefficient_B0,sei=St_err_B0,data=resL1[,1:3])
meta_L1_b1<-rma(Coefficient_B1,sei=St_err_B1,data=resL1)
results_meta_L1<-data.frame("Increase_2002"=100*(exp(meta_L1_b0$b*10 + meta_L1_b1$b*2*10)-exp(10*meta_L1_b0$b)),"Increase_2015"=100*(exp(meta_L1_b0$b*10 + meta_L1_b1$b*15*10)-exp(10*meta_L1_b0$b)))
results_meta_L1$Temp_change<-results_meta_L1$Increase_2015-results_meta_L1$Increase_2002
print(results_meta_L1)## Increase_2002 Increase_2015 Temp_change
## intrcpt 0.10206 0.7679939 0.6659339Results For Lag2 exposure
## For Lag2
resL2 <- data.frame("ville"=villes_name,"Coefficient_B0" = as.numeric(as.character(tot_L2[, 2])),"St_err_B0" = st_errb0_L2,"Coefficient_B1" = as.numeric(as.character(tot_L2[, 5])),"St_err_B1" = st_errb1_L2)
meta_L2_b0<-rma(Coefficient_B0,sei=St_err_B0,data=resL2[,1:3])
meta_L2_b1<-rma(Coefficient_B1,sei=St_err_B1,data=resL2)
results_meta_L2<-data.frame("Increase_2002"=100*(exp(meta_L2_b0$b*10 + meta_L2_b1$b*2*10)-exp(10*meta_L2_b0$b)),"Increase_2015"=100*(exp(meta_L2_b0$b*10 + meta_L2_b1$b*15*10)-exp(10*meta_L2_b0$b)))
results_meta_L2$Temp_change<-results_meta_L2$Increase_2015-results_meta_L2$Increase_2002
print(results_meta_L2)## Increase_2002 Increase_2015 Temp_change
## intrcpt 0.1051875 0.7916099 0.6864224Non-Linear Interaction Model
Models For Lag0, Lag1 and La2 exposure
# define model objects
nl_coeffB0 <- matrix(NA, nrow = 12, ncol = 4)
nl_coeffB1 <- matrix(NA, nrow = 12, ncol = 13)
colnames(nl_coeffB0) <- c("Ville","Coefficient_B0","B0_CI95_Low","B0_CI95_High")
colnames(nl_coeffB1) <- c("Ville","Coefficient1","C1_CI95_Low","C1_CI95_High","Coefficient2","C2_CI95_Low","C2_CI95_High","Coefficient3","C3_CI95_Low","C3_CI95_High","Coefficient4","C4_CI95_Low","C4_CI95_High")
nl_modtot<-list()
year<-list()
pred_nnlin<-list()
y_nl<-matrix(NA,nrow=length(villes),ncol=4)
S_nl<-list()
## FOR LAG1
nl_coeffB0_L1 <- matrix(NA, nrow = 12, ncol = 4)
nl_coeffB1_L1 <- matrix(NA, nrow = 12, ncol = 13)
colnames(nl_coeffB0_L1) <- c("Ville","Coefficient_B0","B0_CI95_Low","B0_CI95_High")
colnames(nl_coeffB1_L1) <- c("Ville","Coefficient1","C1_CI95_Low","C1_CI95_High","Coefficient2","C2_CI95_Low","C2_CI95_High","Coefficient3","C3_CI95_Low","C3_CI95_High","Coefficient4","C4_CI95_Low","C4_CI95_High")
nl_modtot_L1<-list()
pred_nnlin_L1<-list()
y_nl_L1<-matrix(NA,nrow=length(villes),ncol=4)
S_nl_L1<-list()
## FOR LAG2
nl_coeffB0_L2 <- matrix(NA, nrow = 12, ncol = 4)
nl_coeffB1_L2 <- matrix(NA, nrow = 12, ncol = 13)
colnames(nl_coeffB0_L2) <- c("Ville","Coefficient_B0","B0_CI95_Low","B0_CI95_High")
colnames(nl_coeffB1_L2) <- c("Ville","Coefficient1","C1_CI95_Low","C1_CI95_High","Coefficient2","C2_CI95_Low","C2_CI95_High","Coefficient3","C3_CI95_Low","C3_CI95_High","Coefficient4","C4_CI95_Low","C4_CI95_High")
nl_modtot_L2<-list()
pred_nnlin_L2<-list()
y_nl_L2<-matrix(NA,nrow=length(villes),ncol=4)
S_nl_L2<-list()# Non-Linear Model for Lag0
for(i in 1:length(villes)) {
year[[i]] <- onebasis(villes[[i]]$year,"ns",knots=c(4,7,11))
yr<-year[[i]]
nl_modtot[[i]]<-gam(nocc_tot~o3+o3:yr+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(nl_modtot[[i]])
nl_coeffB0[i,1] <- as.character(villes_name[[i]])
nl_coeffB1[i,1] <- as.character(villes_name[[i]])
ll<-length(nl_modtot[[i]]$coefficients)
nl_coeffB0[i, 2] <- nl_modtot[[i]]$coefficients[2]
nl_coeffB0[i, 3] <- nl_modtot[[i]]$coefficients[2]-(1.96*est$se[2])
nl_coeffB0[i, 4] <- nl_modtot[[i]]$coefficients[2]+(1.96*est$se[2])
nl_coeffB1[i, 2] <- nl_modtot[[i]]$coefficients[ll-3]
nl_coeffB1[i, 3] <- nl_modtot[[i]]$coefficients[ll-3]-(1.96*est$se[ll-3])
nl_coeffB1[i, 4] <- nl_modtot[[i]]$coefficients[ll-3]+(1.96*est$se[ll-3])
nl_coeffB1[i, 5] <- nl_modtot[[i]]$coefficients[ll-2]
nl_coeffB1[i, 6] <- nl_modtot[[i]]$coefficients[ll-2]-(1.96*est$se[ll-2])
nl_coeffB1[i, 7] <- nl_modtot[[i]]$coefficients[ll-2]+(1.96*est$se[ll-2])
nl_coeffB1[i, 8] <- nl_modtot[[i]]$coefficients[ll-1]
nl_coeffB1[i, 9] <- nl_modtot[[i]]$coefficients[ll-1]-(1.96*est$se[ll-1])
nl_coeffB1[i, 10] <- nl_modtot[[i]]$coefficients[ll]+(1.96*est$se[ll-1])
nl_coeffB1[i, 11] <- nl_modtot[[i]]$coefficients[ll]
nl_coeffB1[i, 12] <- nl_modtot[[i]]$coefficients[ll]-(1.96*est$se[ll])
nl_coeffB1[i, 13] <- nl_modtot[[i]]$coefficients[ll]+(1.96*est$se[ll])
pred_nnlin[[i]]<-crosspred(yr,nl_modtot[[i]],cen=0)
y_nl[i,] <- pred_nnlin[[i]]$coef
S_nl[[i]] <- pred_nnlin[[i]]$vcov
}
# Non-Linear Model for Lag1
for(i in 1:length(villes)) {
year[[i]] <- onebasis(villes[[i]]$year,"ns",knots=c(4,7,11))
yr<-year[[i]]
nl_modtot_L1[[i]]<-gam(nocc_tot~o3_L1+o3_L1:yr+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(nl_modtot_L1[[i]])
nl_coeffB0_L1[i,1] <- as.character(villes_name[[i]])
nl_coeffB1_L1[i,1] <- as.character(villes_name[[i]])
ll<-length(nl_modtot_L1[[i]]$coefficients)
nl_coeffB0_L1[i, 2] <- nl_modtot_L1[[i]]$coefficients[2]
nl_coeffB0_L1[i, 3] <- nl_modtot_L1[[i]]$coefficients[2]-(1.96*est$se[2])
nl_coeffB0_L1[i, 4] <- nl_modtot_L1[[i]]$coefficients[2]+(1.96*est$se[2])
nl_coeffB1_L1[i, 2] <- nl_modtot_L1[[i]]$coefficients[ll-3]
nl_coeffB1_L1[i, 3] <- nl_modtot_L1[[i]]$coefficients[ll-3]-(1.96*est$se[ll-3])
nl_coeffB1_L1[i, 4] <- nl_modtot_L1[[i]]$coefficients[ll-3]+(1.96*est$se[ll-3])
nl_coeffB1_L1[i, 5] <- nl_modtot_L1[[i]]$coefficients[ll-2]
nl_coeffB1_L1[i, 6] <- nl_modtot_L1[[i]]$coefficients[ll-2]-(1.96*est$se[ll-2])
nl_coeffB1_L1[i, 7] <- nl_modtot_L1[[i]]$coefficients[ll-2]+(1.96*est$se[ll-2])
nl_coeffB1_L1[i, 8] <- nl_modtot_L1[[i]]$coefficients[ll-1]
nl_coeffB1_L1[i, 9] <- nl_modtot_L1[[i]]$coefficients[ll-1]-(1.96*est$se[ll-1])
nl_coeffB1_L1[i, 10] <- nl_modtot_L1[[i]]$coefficients[ll]+(1.96*est$se[ll-1])
nl_coeffB1_L1[i, 11] <- nl_modtot_L1[[i]]$coefficients[ll]
nl_coeffB1_L1[i, 12] <- nl_modtot_L1[[i]]$coefficients[ll]-(1.96*est$se[ll])
nl_coeffB1_L1[i, 13] <- nl_modtot_L1[[i]]$coefficients[ll]+(1.96*est$se[ll])
pred_nnlin_L1[[i]]<-crosspred(yr,nl_modtot_L1[[i]],cen=0)
y_nl_L1[i,] <- pred_nnlin_L1[[i]]$coef
S_nl_L1[[i]] <- pred_nnlin_L1[[i]]$vcov
}
# Non-Linear Model for Lag2
for(i in 1:length(villes)) {
year[[i]] <- onebasis(villes[[i]]$year,"ns",knots=c(4,7,11))
yr<-year[[i]]
nl_modtot_L2[[i]]<-gam(nocc_tot~o3_L2+o3_L2:yr+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(nl_modtot_L2[[i]])
nl_coeffB0_L2[i,1] <- as.character(villes_name[[i]])
nl_coeffB1_L2[i,1] <- as.character(villes_name[[i]])
ll<-length(nl_modtot_L2[[i]]$coefficients)
nl_coeffB0_L2[i, 2] <- nl_modtot_L2[[i]]$coefficients[2]
nl_coeffB0_L2[i, 3] <- nl_modtot_L2[[i]]$coefficients[2]-(1.96*est$se[2])
nl_coeffB0_L2[i, 4] <- nl_modtot_L2[[i]]$coefficients[2]+(1.96*est$se[2])
nl_coeffB1_L2[i, 2] <- nl_modtot_L2[[i]]$coefficients[ll-3]
nl_coeffB1_L2[i, 3] <- nl_modtot_L2[[i]]$coefficients[ll-3]-(1.96*est$se[ll-3])
nl_coeffB1_L2[i, 4] <- nl_modtot_L2[[i]]$coefficients[ll-3]+(1.96*est$se[ll-3])
nl_coeffB1_L2[i, 5] <- nl_modtot_L2[[i]]$coefficients[ll-2]
nl_coeffB1_L2[i, 6] <- nl_modtot_L2[[i]]$coefficients[ll-2]-(1.96*est$se[ll-2])
nl_coeffB1_L2[i, 7] <- nl_modtot_L2[[i]]$coefficients[ll-2]+(1.96*est$se[ll-2])
nl_coeffB1_L2[i, 8] <- nl_modtot_L2[[i]]$coefficients[ll-1]
nl_coeffB1_L2[i, 9] <- nl_modtot_L2[[i]]$coefficients[ll-1]-(1.96*est$se[ll-1])
nl_coeffB1_L2[i, 10] <- nl_modtot_L2[[i]]$coefficients[ll]+(1.96*est$se[ll-1])
nl_coeffB1_L2[i, 11] <- nl_modtot_L2[[i]]$coefficients[ll]
nl_coeffB1_L2[i, 12] <- nl_modtot_L2[[i]]$coefficients[ll]-(1.96*est$se[ll])
nl_coeffB1_L2[i, 13] <- nl_modtot_L2[[i]]$coefficients[ll]+(1.96*est$se[ll])
pred_nnlin_L2[[i]]<-crosspred(yr,nl_modtot_L2[[i]],cen=0)
y_nl_L2[i,] <- pred_nnlin_L2[[i]]$coef
S_nl_L2[[i]] <- pred_nnlin_L2[[i]]$vcov
}Meta-analysis results
meta_nl <- mvmeta(y_nl,S_nl,method="ml")
sink("C:/Users/Anna/Dropbox/Leo/change DR curves/results/Ozone/meta_nl.txt")
summary(meta_nl)## Call: mvmeta(formula = y_nl ~ 1, S = S_nl, method = "ml")
##
## Multivariate random-effects meta-analysis
## Dimension: 4
## Estimation method: ML
##
## Fixed-effects coefficients
## Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
## y1 -0.0006 0.0002 -3.2474 0.0012 -0.0010 -0.0002 **
## y2 0.0000 0.0002 0.0424 0.9662 -0.0004 0.0004
## y3 -0.0008 0.0005 -1.5940 0.1109 -0.0017 0.0002
## y4 0.0016 0.0003 4.8023 0.0000 0.0009 0.0022 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Between-study random-effects (co)variance components
## Structure: General positive-definite
## Std. Dev Corr
## y1 0.0005 y1 y2 y3
## y2 0.0007 0.8431
## y3 0.0015 0.9817 0.9300
## y4 0.0010 0.8756 0.9980 0.9515
##
## Multivariate Cochran Q-test for heterogeneity:
## Q = 344.1948 (df = 44), p-value = 0.0000
## I-square statistic = 87.2%
##
## 12 studies, 48 observations, 4 fixed and 10 random-effects parameters
## logLik AIC BIC
## 287.0542 -546.1084 -519.9115sink()
#plot pooled RR - nonlinear
yearplot <- seq(0,14)
yrplt <- onebasis(yearplot,"ns",knots=c(4,7,11))
plot_nl <- crosspred(yrplt,coef=coef(meta_nl),vcov=vcov(meta_nl),cen=0,model.link="log",by=1)
plot(plot_nl,"overall",ci="lines",col=1,lwd=3,xaxt="n",ylab="% increase in risk",xlab="year",main="L0")
axis(1, at=1:14, labels=as.character(2002:2015))
meta_nl_L1 <- mvmeta(y_nl_L1,S_nl_L1,method="ml")
sink("C:/Users/Anna/Dropbox/Leo/change DR curves/results/Ozone/meta_nl_L1.txt")
summary(meta_nl_L1)## Call: mvmeta(formula = y_nl_L1 ~ 1, S = S_nl_L1, method = "ml")
##
## Multivariate random-effects meta-analysis
## Dimension: 4
## Estimation method: ML
##
## Fixed-effects coefficients
## Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
## y1 -0.0004 0.0001 -3.5589 0.0004 -0.0007 -0.0002 ***
## y2 -0.0000 0.0002 -0.1764 0.8599 -0.0003 0.0003
## y3 -0.0006 0.0003 -1.8115 0.0701 -0.0012 0.0000 .
## y4 0.0012 0.0002 4.8836 0.0000 0.0007 0.0016 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Between-study random-effects (co)variance components
## Structure: General positive-definite
## Std. Dev Corr
## y1 0.0003 y1 y2 y3
## y2 0.0004 0.8558
## y3 0.0010 0.9897 0.9210
## y4 0.0008 0.8840 0.9984 0.9418
##
## Multivariate Cochran Q-test for heterogeneity:
## Q = 370.9699 (df = 44), p-value = 0.0000
## I-square statistic = 88.1%
##
## 12 studies, 48 observations, 4 fixed and 10 random-effects parameters
## logLik AIC BIC
## 304.5128 -581.0256 -554.8288sink()
plot_nl_L1 <- crosspred(yrplt,coef=coef(meta_nl_L1),vcov=vcov(meta_nl_L1),cen=0,model.link="log",by=1)
plot(plot_nl_L1,"overall",ci="lines",col=1,lwd=3,xaxt="n",ylab="% increase in risk",xlab="year",main="L1")
axis(1, at=1:14, labels=as.character(2002:2015))
meta_nl_L2 <- mvmeta(y_nl_L2,S_nl_L2,method="ml")
sink("C:/Users/Anna/Dropbox/Leo/change DR curves/results/Ozone/meta_nl_L2.txt")
summary(meta_nl_L2)## Call: mvmeta(formula = y_nl_L2 ~ 1, S = S_nl_L2, method = "ml")
##
## Multivariate random-effects meta-analysis
## Dimension: 4
## Estimation method: ML
##
## Fixed-effects coefficients
## Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
## y1 -0.0004 0.0001 -3.0633 0.0022 -0.0007 -0.0002 **
## y2 -0.0000 0.0002 -0.1002 0.9202 -0.0003 0.0003
## y3 -0.0005 0.0003 -1.6841 0.0922 -0.0012 0.0001 .
## y4 0.0012 0.0002 4.7868 0.0000 0.0007 0.0017 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Between-study random-effects (co)variance components
## Structure: General positive-definite
## Std. Dev Corr
## y1 0.0004 y1 y2 y3
## y2 0.0004 0.7723
## y3 0.0010 0.9854 0.8693
## y4 0.0008 0.8810 0.9809 0.9488
##
## Multivariate Cochran Q-test for heterogeneity:
## Q = 388.4143 (df = 44), p-value = 0.0000
## I-square statistic = 88.7%
##
## 12 studies, 48 observations, 4 fixed and 10 random-effects parameters
## logLik AIC BIC
## 304.5209 -581.0418 -554.8450sink()
plot_nl_L2 <- crosspred(yrplt,coef=coef(meta_nl_L2),vcov=vcov(meta_nl_L2),cen=0,model.link="log",by=1)
plot(plot_nl_L2,"overall",ci="lines",col=1,lwd=3,xaxt="n",ylab="% increase in risk",xlab="year",main="L2")
axis(1, at=1:14, labels=as.character(2002:2015))