Diff in Diff

Seuil d’Information et d’Alerte selon les années et le polluant

Les seuils réglementés pour NO2 n’ont pas changé au cours de la période observée.
Le seuil d’alerte pour l’ozone a changé en 2006.
La réglementation a inclut des seuils pour les PM10 seulement à partir de la fin de l’année 2007 et leur seuil d’alerte a été ensuite modifié en 2012

Seuil d’Information et Alerte par année et Polluant
	O3		NO2		PM10
annee	Seuil_Info	Seuil_Alerte	Seuil_Inf	Seuil_Alert	S_Info	S_Alert
2000	180	360	200	400	NA	NA
2001	180	360	200	400	NA	NA
2002	180	360	200	400	NA	NA
2003	180	360	200	400	NA	NA
2004	180	360	200	400	NA	NA
2005	180	360	200	400	NA	NA
2006	180	240	200	400	NA	NA
2007	180	240	200	400	NA	NA
2008	180	240	200	400	80	125
2009	180	240	200	400	80	125
2010	180	240	200	400	80	125
2011	180	240	200	400	80	125
2012	180	240	200	400	50	80
2013	180	240	200	400	50	80
2014	180	240	200	400	50	80
2015	180	240	200	400	50	80
2016	180	240	200	400	50	80
2017	180	240	200	400	50	80
Note:
Les couleurs indiquent les différentes périodes selon les changements de réglementation

Les crithères d’éligibilité aux mesures anti-pic de pollution ont donc changé dans le temps.

Création variable éligibilité aux mesures anti-pic de pollution

data$elig<-NA

# Periode 2000-2006
for (i in 1 : 2562) {
data$elig[i]<- if (data$o3[i]>data$O3_Inf[i]) 1 else 0
}
# Année 2007
for (i in 2563 : 2927) {
  data$elig[i]<- if (data$o3[i]>data$O3_Inf[i] | data$no2[i]>data$NO2_Inf[i]) 1 else 0
}
# Periode 2008-fin
for (i in 2928 : nrow(data)) {
  data$elig[i]<- if (data$pm10full[i]>data$PM10_Inf[i] | data$o3[i]>data$O3_Inf[i] | data$no2[i]>data$NO2_Inf[i]) 1 else 0
}

Nombre de jours éligibles sur toute la période d’étude
Eligibilité	Nombre de Jours
0	6476
1	99

Nombre de jours éligibles par année
	2000	2001	2002	2003	2004	2005	2006	2007	2008	2009	2010	2011	2012	2013	2014	2015	2016	2017
0	366	365	365	360	366	363	364	365	366	362	364	362	339	343	357	357	353	359
1	0	0	0	5	0	2	1	0	0	3	1	3	27	22	8	8	13	6

Dans le but de notre analyse, on va considérer 3 périodes: - de 2000 à 2005 : NO2 et O3 réglementés - de 2006 à 2007 : NO2 et O3 réglementés, nouveau seuil pour O3 - à partir de 2008 : PM10 inclus dans la réglementation

Nombre de jours éligibles pendant la première période (2000-2005)
Eligibilité	Nombre de Jours
0	2190
1	7

Jours éligibles au cours de la première périodes et concentrations moyennes journalières de O3 correspondantes:

Date	O3_Concentration
2003-07-15	193.0676
2003-08-03	187.7442
2003-08-07	193.4258
2003-08-08	216.3857
2003-08-12	189.4658
2005-07-14	187.7017
2005-07-15	180.6635

Nombre de jours éligibles pendant la deuxième période (2006-2007)
Eligibilité	Nombre de Jours
0	729
1	1

Jours éligibles au cours de la deuxième période et concentrations moyennes journalières de O3 et NO2 correspondantes:

Date	O3_Concentration	NO2_Concentration
2006-07-26	188.0962	NA

Nombre de jours éligibles pendant la trosième période (à partir de 2008)
Eligibilité	Nombre de Jours
0	3557
1	91

Creation d’une variable indicatrice de période:

Creation variable Période

data$annee<-as.numeric(as.character(data$annee))
data$period<-NA
for (i in 1:nrow(data)) {
  data$period[i]<-if (data$annee[i]<2006) "1" 
  else if (data$annee[i]>2005 & data$annee[i]<2008) "2" 
  else if (data$annee[i]>2007) "3"
  }

Calcul Score de Propension

Deux modèles sont comparés pour calculer le score de propension:

un modèle incluant un seul lag pour la variable ozone (Mod1)
un modèle incluant 3 lag pour la variable ozone (Mod2)

mod1<-glm(elig~tempmin+tempmax+Jours+mois+jours_fériés+Vacances+annee+lag_O3,family="binomial",data=data)
mod2<-glm(elig~tempmin+tempmax+Jours+mois+jours_fériés+Vacances+annee+lag_O3+lag2_O3+lag3_O3,family="binomial",data=data)

data$SP1[1]<-NA
data$SP1[2:nrow(data)]<-mod1$fitted
data$SP2<-NA
data$SP2[4:nrow(data)]<-mod2$fitted

Résultats Modèle 1 avec un seul LAG pour Ozone

Observations	6574 (1 missing obs. deleted)
Dependent variable	elig
Type	Generalized linear model
Family	binomial
Link	logit

χ²(39)	443.375
Pseudo-R² (Cragg-Uhler)	0.451
Pseudo-R² (McFadden)	0.432
AIC	663.887
BIC	935.522

	Est.	2.5%	97.5%	z val.	p
(Intercept)	-24.193	-4302.833	4254.448	-0.011	0.991
tempmin	-0.483	-0.596	-0.369	-8.360	0.000
tempmax	0.245	0.159	0.330	5.615	0.000
Joursjeudi	1.572	0.566	2.578	3.063	0.002
Jourslundi	0.017	-1.162	1.196	0.028	0.977
Joursmardi	1.200	0.155	2.245	2.252	0.024
Joursmercredi	1.268	0.233	2.303	2.402	0.016
Jourssamedi	1.014	-0.048	2.075	1.872	0.061
Joursvendredi	1.454	0.437	2.470	2.803	0.005
mois10	-17.640	-3398.031	3362.751	-0.010	0.992
mois11	-1.080	-2.412	0.252	-1.589	0.112
mois12	-0.029	-0.876	0.819	-0.067	0.947
mois2	-0.451	-1.294	0.392	-1.048	0.295
mois3	0.294	-0.511	1.098	0.716	0.474
mois4	-1.953	-3.170	-0.736	-3.146	0.002
mois5	-2.931	-4.685	-1.177	-3.276	0.001
mois6	-18.895	-3453.653	3415.864	-0.011	0.991
mois7	-0.409	-2.264	1.445	-0.433	0.665
mois8	-0.216	-2.103	1.672	-0.224	0.823
mois9	-18.630	-3379.179	3341.919	-0.011	0.991
jours_fériés	0.895	-0.666	2.455	1.123	0.261
Vacances	-1.360	-2.160	-0.559	-3.328	0.001
annee2001	0.098	-6132.310	6132.506	0.000	1.000
annee2002	-0.087	-6075.926	6075.751	-0.000	1.000
annee2003	17.358	-4261.283	4295.998	0.008	0.994
annee2004	-0.394	-6001.263	6000.475	-0.000	1.000
annee2005	16.817	-4261.824	4295.458	0.008	0.994
annee2006	15.734	-4262.907	4294.375	0.007	0.994
annee2007	-0.125	-6080.325	6080.075	-0.000	1.000
annee2008	-0.420	-6041.351	6040.512	-0.000	1.000
annee2009	16.717	-4261.923	4295.358	0.008	0.994
annee2010	15.679	-4262.962	4294.320	0.007	0.994
annee2011	16.982	-4261.658	4295.623	0.008	0.994
annee2012	19.528	-4259.113	4298.168	0.009	0.993
annee2013	19.309	-4259.331	4297.950	0.009	0.993
annee2014	18.064	-4260.577	4296.704	0.008	0.993
annee2015	18.010	-4260.631	4296.650	0.008	0.993
annee2016	18.743	-4259.898	4297.383	0.009	0.993
annee2017	17.782	-4260.859	4296.423	0.008	0.994
lag_O3	0.033	0.021	0.045	5.442	0.000
Standard errors: MLE

Résultats Modèle 2 avec 3 LAG pour Ozone

Observations	6572 (3 missing obs. deleted)
Dependent variable	elig
Type	Generalized linear model
Family	binomial
Link	logit

χ²(41)	443.868
Pseudo-R² (Cragg-Uhler)	0.451
Pseudo-R² (McFadden)	0.432
AIC	667.333
BIC	952.537

	Est.	2.5%	97.5%	z val.	p
(Intercept)	-24.277	-4313.208	4264.654	-0.011	0.991
tempmin	-0.485	-0.598	-0.371	-8.369	0.000
tempmax	0.245	0.160	0.331	5.609	0.000
Joursjeudi	1.570	0.564	2.577	3.058	0.002
Jourslundi	0.044	-1.136	1.224	0.072	0.942
Joursmardi	1.169	0.121	2.217	2.187	0.029
Joursmercredi	1.258	0.220	2.296	2.376	0.017
Jourssamedi	1.015	-0.047	2.076	1.874	0.061
Joursvendredi	1.456	0.440	2.473	2.809	0.005
mois10	-17.644	-3395.029	3359.741	-0.010	0.992
mois11	-1.067	-2.402	0.268	-1.566	0.117
mois12	-0.033	-0.884	0.817	-0.077	0.939
mois2	-0.453	-1.297	0.391	-1.052	0.293
mois3	0.269	-0.540	1.079	0.652	0.514
mois4	-2.021	-3.256	-0.785	-3.206	0.001
mois5	-3.020	-4.796	-1.243	-3.332	0.001
mois6	-18.947	-3464.281	3426.387	-0.011	0.991
mois7	-0.512	-2.394	1.370	-0.533	0.594
mois8	-0.279	-2.176	1.619	-0.288	0.773
mois9	-18.670	-3380.528	3343.188	-0.011	0.991
jours_fériés	0.889	-0.675	2.453	1.114	0.265
Vacances	-1.353	-2.153	-0.553	-3.315	0.001
annee2001	0.101	-6139.850	6140.053	0.000	1.000
annee2002	-0.082	-6086.391	6086.227	-0.000	1.000
annee2003	17.335	-4271.596	4306.266	0.008	0.994
annee2004	-0.399	-6009.281	6008.483	-0.000	1.000
annee2005	16.813	-4272.118	4305.745	0.008	0.994
annee2006	15.726	-4273.206	4304.657	0.007	0.994
annee2007	-0.140	-6084.999	6084.719	-0.000	1.000
annee2008	-0.426	-6048.421	6047.569	-0.000	1.000
annee2009	16.718	-4272.214	4305.649	0.008	0.994
annee2010	15.671	-4273.261	4304.602	0.007	0.994
annee2011	16.977	-4271.954	4305.908	0.008	0.994
annee2012	19.528	-4269.403	4308.459	0.009	0.993
annee2013	19.303	-4269.628	4308.234	0.009	0.993
annee2014	18.058	-4270.873	4306.989	0.008	0.993
annee2015	17.998	-4270.933	4306.929	0.008	0.993
annee2016	18.761	-4270.170	4307.692	0.009	0.993
annee2017	17.785	-4271.146	4306.717	0.008	0.994
lag_O3	0.029	0.012	0.046	3.374	0.001
lag2_O3	0.005	-0.016	0.026	0.481	0.630
lag3_O3	0.002	-0.016	0.019	0.187	0.852
Standard errors: MLE

Statistiques Descriptives pour le Score de Propension:

Modèle 1

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.      NA's 
## 0.0000000 0.0000000 0.0000125 0.0150593 0.0045602 0.8229836         1

Modèle 2

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.      NA's 
## 0.0000000 0.0000000 0.0000151 0.0150639 0.0045756 0.8009219         3

Le modèle 1 présente un AIC inferieur et un seule prédiction manquante (a cause des observations correspondantes au premier lag de l’ozone qui ont été supprimées dans le modèle). Le SP predit par le modèle 2 présente 3 estimations manquantes.

Test moyenne Score de Propension Test de Student pour difference moyenne entre variable Score de Propension (x) et variable éligibilité (y):

Modèle 1:

## 
##  Welch Two Sample t-test
## 
## data:  data$SP1 and data$elig
## t = 0.0013819, df = 9308.2, p-value = 0.9989
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.003246773  0.003251354
## sample estimates:
##  mean of x  mean of y 
## 0.01505932 0.01505703

Modèle 2:

## 
##  Welch Two Sample t-test
## 
## data:  data$SP2 and data$elig
## t = 0.0041457, df = 9315.5, p-value = 0.9967
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.003243051  0.003256798
## sample estimates:
##  mean of x  mean of y 
## 0.01506391 0.01505703

La moyenne du score de propension predit par le modèle 1 est plus proche de la moyenne de la variable éligibilité que la moyenne du SP predit par le modèle 2.

Comparaison de la moyenne du score de propension entre jours éligibles et pas éligibles avec test de Student.

Modèle 1:

## 
##  Welch Two Sample t-test
## 
## data:  SP1 by elig
## t = -10.975, df = 98.165, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.2508522 -0.1740275
## sample estimates:
## mean in group 0 mean in group 1 
##      0.01186012      0.22429993

Modèle 2:

## 
##  Welch Two Sample t-test
## 
## data:  SP2 by elig
## t = -11.013, df = 98.166, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.2513655 -0.1746089
## sample estimates:
## mean in group 0 mean in group 1 
##      0.01185549      0.22484265

Standard Deviation de la variable SP et Eligibilité:

Modèle 1:

## data$elig: 0
## [1] 0.04515635
## ------------------------------------------------------------ 
## data$elig: 1
## [1] 0.1925178

Modèle 2:

## data$elig: 0
## [1] 0.045221
## ------------------------------------------------------------ 
## data$elig: 1
## [1] 0.1923468

Appariement par le Score de Propension

Modèle 1

data1<-data[-1,]
var.match<-cbind(data1$SP1,data1$period)
App1<-Match(Tr = data1$elig,X=var.match, ties=FALSE)

summary(App1)

## 
## Estimate...  0 
## SE.........  0 
## T-stat.....  NaN 
## p.val......  NA 
## 
## Original number of observations..............  6574 
## Original number of treated obs...............  99 
## Matched number of observations...............  99 
## Matched number of observations  (unweighted).  99

control=data.frame(App1$index.control)
treated=data.frame(App1$index.treated)

data1$nrow<-c(1:nrow(data1))

colnames(control) <- c("nrow")
colnames(treated) <- c("nrow")

cas<- cbind(treated, PAIR = c(1:99))
control <- cbind(control, PAIR = c(1:99))

casm<-merge(cas,data1,by="nrow")
controlm<- merge(control,data1, by ="nrow")

datam<-rbind(casm,controlm)
datam<-datam[order(datam[,"PAIR"]),]

Modèle 2:

data2<-data[-c(1:3),]
var.match<-cbind(data2$SP1,data2$period)
App2<-Match(Tr = data2$elig,X=var.match, ties=FALSE)

summary(App2)

## 
## Estimate...  0 
## SE.........  0 
## T-stat.....  NaN 
## p.val......  NA 
## 
## Original number of observations..............  6572 
## Original number of treated obs...............  99 
## Matched number of observations...............  99 
## Matched number of observations  (unweighted).  99

control=data.frame(App2$index.control)
treated=data.frame(App2$index.treated)

data2$nrow<-c(1:nrow(data2))

colnames(control) <- c("nrow")
colnames(treated) <- c("nrow")

cas<- cbind(treated, PAIR = c(1:99))
control <- cbind(control, PAIR = c(1:99))

casm<-merge(cas,data2,by="nrow")
controlm<- merge(control,data2, by ="nrow")

datam2<-rbind(casm,controlm)
datam2<-datam2[order(datam2[,"PAIR"]),]

Comparaison SP entre jours éligibles et jours pas éligibles Après et Avant Appariement:

Modèle 1:

Modèle 2:

Comparaison de la moyenne du score de propension entre jours éligibles et pas éligibles après appariement avec test de Student:

Modèle 1:

## 
##  Welch Two Sample t-test
## 
## data:  SP1 by elig
## t = -0.0052813, df = 196, p-value = 0.9958
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.05399377  0.05370535
## sample estimates:
## mean in group 0 mean in group 1 
##       0.2241557       0.2242999

Modèle 2:

## 
##  Welch Two Sample t-test
## 
## data:  SP2 by elig
## t = -0.049393, df = 196, p-value = 0.9607
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.05514030  0.05244579
## sample estimates:
## mean in group 0 mean in group 1 
##       0.2234954       0.2248427

La valeur des moyennes du SP entre jours éligibles et jours non éligibles après appariement sont plus proches avec appariement par SP du modèle 1.

Balance statistics for matching

Standardized mean differences for continuos variables after matching with SP1 and SP2:

	temp min	temp max	lag_O3
SP mod1	0.0663546	0.0603000	0.0832185
SP mod2	0.0791104	0.0659674	0.0677056

Love Plot for standardized mean difference

Differences in Differences

Après création de 2 dataframes par type d’appariement incluant uniquement la période 1 et la période 3 (pour la période 2 on a uniquement 1 jour éligible), analyse diff in diff via regression linéaire pour évaluer l’effet des mesures sur le niveau de 03 et via regression de Poisson pour évaluer effet sur mortalité non accidentelle:

library(geepack)

dfdf_o3.1<-geeglm(formula=o3~period+elig+period*elig, family=gaussian, id=PAIR, data=data1_Ap1)
dfdf_o3.2<-geeglm(formula=o3~period+elig+period*elig, family=gaussian, id=PAIR, data=data1_Ap2)


dfdf_mort1<-geeglm(formula=nocc_tot~period+elig+period*elig, family=poisson(link = "log"), id=PAIR, data=data1_Ap1)
dfdf_mort2<-geeglm(formula=nocc_tot~period+elig+period*elig, family=poisson(link = "log"), id=PAIR, data=data1_Ap2)

Effet sur niveau de concentration d’ozone période 1 vs. Période post mésures, avec appariement par SP1 :

## 
## Call:
## geeglm(formula = o3 ~ period + elig + period * elig, family = gaussian, 
##     data = data1_Ap1, id = PAIR)
## 
##  Coefficients:
##              Estimate  Std.err   Wald Pr(>|W|)    
## (Intercept)    62.494   21.731  8.270  0.00403 ** 
## period3        -7.672   21.890  0.123  0.72597    
## elig          130.142   23.742 30.048 4.21e-08 ***
## period3:elig -140.335   24.072 33.986 5.55e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = independence 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)      760   91.44
## Number of clusters:   98  Maximum cluster size: 2

Effet sur niveau de concentration d’ozone période 1 vs. Période post mésures, avec appariement par SP2 :

## 
## Call:
## geeglm(formula = o3 ~ period + elig + period * elig, family = gaussian, 
##     data = data1_Ap2, id = PAIR)
## 
##  Coefficients:
##              Estimate Std.err   Wald Pr(>|W|)    
## (Intercept)      35.1    10.2  11.88  0.00057 ***
## period3          19.9    10.5   3.54  0.05979 .  
## elig            157.5    10.5 223.09  < 2e-16 ***
## period3:elig   -167.9    11.2 224.27  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = independence 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)      691    68.9
## Number of clusters:   98  Maximum cluster size: 2

Effet sur mortalité non accidentelle période 1 vs. Période post mésures, avec appariement par SP1 :

## 
## Call:
## geeglm(formula = nocc_tot ~ period + elig + period * elig, family = poisson(link = "log"), 
##     data = data1_Ap1, id = PAIR)
## 
##  Coefficients:
##              Estimate Std.err    Wald Pr(>|W|)    
## (Intercept)    4.8708  0.0798 3727.37   <2e-16 ***
## period3       -0.1620  0.0812    3.99    0.046 *  
## elig           0.6676  0.3424    3.80    0.051 .  
## period3:elig  -0.6463  0.3429    3.55    0.059 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = independence 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)     9.38    5.39
## Number of clusters:   91  Maximum cluster size: 2

Effet sur mortalité non accidentelle période 2 vs. Période post mésures, avec appariement par SP2 :

## 
## Call:
## geeglm(formula = nocc_tot ~ period + elig + period * elig, family = poisson(link = "log"), 
##     data = data1_Ap2, id = PAIR)
## 
##  Coefficients:
##              Estimate Std.err     Wald Pr(>|W|)    
## (Intercept)    4.8329  0.0253 36552.17  < 2e-16 ***
## period3       -0.1270  0.0289    19.39  1.1e-05 ***
## elig           0.7056  0.3303     4.56    0.033 *  
## period3:elig  -0.6813  0.3309     4.24    0.040 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = independence 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)     9.12    5.44
## Number of clusters:   89  Maximum cluster size: 2

Diff in Diff

Anna Alari

06/02/2020

Seuil d’Information et d’Alerte selon les années et le polluant

Création variable éligibilité aux mesures anti-pic de pollution

Creation variable Période

Calcul Score de Propension

Résultats Modèle 1 avec un seul LAG pour Ozone

Résultats Modèle 2 avec 3 LAG pour Ozone

Statistiques Descriptives pour le Score de Propension:

Appariement par le Score de Propension

Modèle 1

Modèle 2:

Balance statistics for matching

Love Plot for standardized mean difference

Differences in Differences

	2000	2001	2002	2003	2004	2005	2006	2007	2008	2009	2010	2011	2012	2013	2014	2015	2016	2017
0	366	365	365	360	366	363	364	365	366	362	364	362	339	343	357	357	353	359
1	0	0	0	5	0	2	1	0	0	3	1	3	27	22	8	8	13	6

	2000	2001	2002	2003	2004	2005	2006	2007	2008	2009	2010	2011	2012	2013	2014	2015	2016	2017
0	366	365	365	360	366	363	364	365	366	362	364	362	339	343	357	357	353	359
1	0	0	0	5	0	2	1	0	0	3	1	3	27	22	8	8	13	6

	2000	2001	2002	2003	2004	2005	2006	2007	2008	2009	2010	2011	2012	2013	2014	2015	2016	2017
0	366	365	365	360	366	363	364	365	366	362	364	362	339	343	357	357	353	359
1	0	0	0	5	0	2	1	0	0	3	1	3	27	22	8	8	13	6