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
}

Exclusion canicule 2008

for (i in (1:nrow(data))) {
data$can[i]<- if (year(data$Dates[i])==2003 & month(data$Dates[i])==8) 1 else 0
}

## Warning: Unknown or uninitialised column: 'can'.

data<-data[which(data$can==0),]

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

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	333	366	363	364	365	366	362	364	362	339	343	357	357	353	359
1	0	0	0	1	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	2158
1	3

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
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	760
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	3562
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	6543 (1 missing obs. deleted)
Dependent variable	elig
Type	Generalized linear model
Family	binomial
Link	logit

χ²(39)	453.761
Pseudo-R² (Cragg-Uhler)	0.476
Pseudo-R² (McFadden)	0.457
AIC	618.985
BIC	890.431

	Est.	2.5%	97.5%	z val.	p
(Intercept)	-24.034	-4229.162	4181.095	-0.011	0.991
tempmin	-0.500	-0.617	-0.384	-8.402	0.000
tempmax	0.245	0.156	0.334	5.403	0.000
Joursjeudi	1.719	0.629	2.809	3.092	0.002
Jourslundi	0.280	-0.966	1.525	0.440	0.660
Joursmardi	1.292	0.165	2.419	2.247	0.025
Joursmercredi	1.448	0.336	2.560	2.553	0.011
Jourssamedi	1.253	0.115	2.391	2.158	0.031
Joursvendredi	1.622	0.521	2.723	2.887	0.004
mois10	-17.448	-3378.296	3343.400	-0.010	0.992
mois11	-1.058	-2.396	0.280	-1.550	0.121
mois12	0.001	-0.848	0.851	0.003	0.997
mois2	-0.481	-1.339	0.377	-1.099	0.272
mois3	0.565	-0.251	1.380	1.357	0.175
mois4	-1.515	-2.742	-0.288	-2.420	0.016
mois5	-2.388	-4.162	-0.615	-2.640	0.008
mois6	-18.226	-3451.917	3415.464	-0.010	0.992
mois7	0.491	-1.355	2.337	0.521	0.602
mois8	-16.808	-3433.855	3400.238	-0.010	0.992
mois9	-18.065	-3405.266	3369.136	-0.010	0.992
jours_fériés	1.046	-0.550	2.642	1.285	0.199
Vacances	-1.343	-2.147	-0.538	-3.270	0.001
annee2001	0.199	-6044.425	6044.823	0.000	1.000
annee2002	0.025	-5995.891	5995.940	0.000	1.000
annee2003	15.847	-4189.282	4220.976	0.007	0.994
annee2004	-0.313	-5926.723	5926.097	-0.000	1.000
annee2005	16.891	-4188.238	4222.020	0.008	0.994
annee2006	15.752	-4189.377	4220.881	0.007	0.994
annee2007	-0.035	-5991.146	5991.075	-0.000	1.000
annee2008	-0.329	-5971.841	5971.182	-0.000	1.000
annee2009	16.763	-4188.366	4221.892	0.008	0.994
annee2010	15.675	-4189.454	4220.804	0.007	0.994
annee2011	17.054	-4188.075	4222.182	0.008	0.994
annee2012	19.627	-4185.502	4224.755	0.009	0.993
annee2013	19.406	-4185.723	4224.534	0.009	0.993
annee2014	18.212	-4186.917	4223.340	0.008	0.993
annee2015	18.079	-4187.049	4223.208	0.008	0.993
annee2016	18.811	-4186.317	4223.940	0.009	0.993
annee2017	17.827	-4187.301	4222.956	0.008	0.993
lag_O3	0.024	0.011	0.037	3.679	0.000
Standard errors: MLE

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

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

χ²(41)	453.871
Pseudo-R² (Cragg-Uhler)	0.476
Pseudo-R² (McFadden)	0.457
AIC	622.816
BIC	907.822

	Est.	2.5%	97.5%	z val.	p
(Intercept)	-24.073	-4239.714	4191.569	-0.011	0.991
tempmin	-0.501	-0.618	-0.384	-8.390	0.000
tempmax	0.244	0.155	0.333	5.364	0.000
Joursjeudi	1.712	0.623	2.801	3.082	0.002
Jourslundi	0.278	-0.969	1.525	0.436	0.663
Joursmardi	1.300	0.171	2.428	2.257	0.024
Joursmercredi	1.432	0.318	2.545	2.521	0.012
Jourssamedi	1.250	0.113	2.387	2.155	0.031
Joursvendredi	1.618	0.518	2.718	2.883	0.004
mois10	-17.441	-3376.562	3341.680	-0.010	0.992
mois11	-1.047	-2.387	0.292	-1.532	0.125
mois12	0.011	-0.840	0.862	0.025	0.980
mois2	-0.477	-1.335	0.382	-1.088	0.277
mois3	0.559	-0.260	1.378	1.338	0.181
mois4	-1.536	-2.780	-0.292	-2.420	0.016
mois5	-2.409	-4.202	-0.616	-2.633	0.008
mois6	-18.250	-3449.633	3413.133	-0.010	0.992
mois7	0.472	-1.394	2.339	0.496	0.620
mois8	-16.816	-3435.151	3401.518	-0.010	0.992
mois9	-18.074	-3405.465	3369.317	-0.010	0.992
jours_fériés	1.040	-0.558	2.639	1.276	0.202
Vacances	-1.340	-2.147	-0.534	-3.257	0.001
annee2001	0.213	-6052.612	6053.037	0.000	1.000
annee2002	0.030	-6003.142	6003.203	0.000	1.000
annee2003	15.856	-4199.786	4231.498	0.007	0.994
annee2004	-0.313	-5933.238	5932.613	-0.000	1.000
annee2005	16.893	-4198.749	4232.535	0.008	0.994
annee2006	15.754	-4199.888	4231.397	0.007	0.994
annee2007	-0.032	-5999.158	5999.093	-0.000	1.000
annee2008	-0.330	-5978.574	5977.913	-0.000	1.000
annee2009	16.765	-4198.877	4232.407	0.008	0.994
annee2010	15.673	-4199.969	4231.315	0.007	0.994
annee2011	17.059	-4198.583	4232.701	0.008	0.994
annee2012	19.629	-4196.013	4235.271	0.009	0.993
annee2013	19.410	-4196.232	4235.052	0.009	0.993
annee2014	18.217	-4197.425	4233.858	0.008	0.993
annee2015	18.082	-4197.560	4233.724	0.008	0.993
annee2016	18.823	-4196.818	4234.465	0.009	0.993
annee2017	17.839	-4197.803	4233.481	0.008	0.993
lag_O3	0.025	0.006	0.043	2.600	0.009
lag2_O3	-0.003	-0.026	0.019	-0.281	0.779
lag3_O3	0.004	-0.015	0.023	0.408	0.683
Standard errors: MLE

Statistiques Descriptives pour le Score de Propension:

Modèle 1

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
## 0.000000 0.000000 0.000000 0.014519 0.003184 0.821557        1

Modèle 2

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
## 0.000000 0.000000 0.000000 0.014524 0.003192 0.820011        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.0013488, df = 9486.4, p-value = 0.9989
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.003222387  0.003226825
## sample estimates:
##  mean of x  mean of y 
## 0.01451933 0.01451711

Modèle 2:

## 
##  Welch Two Sample t-test
## 
## data:  data$SP2 and data$elig
## t = 0.0040469, df = 9490.5, p-value = 0.9968
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.003218439  0.003231756
## sample estimates:
##  mean of x  mean of y 
## 0.01452377 0.01451711

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 = -11.59, df = 94.154, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.2717117 -0.1922308
## sample estimates:
## mean in group 0 mean in group 1 
##      0.01115127      0.24312249

Modèle 2:

## 
##  Welch Two Sample t-test
## 
## data:  SP2 by elig
## t = -11.602, df = 94.155, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.2717354 -0.1923187
## sample estimates:
## mean in group 0 mean in group 1 
##      0.01115387      0.24318090

Standard Deviation de la variable SP et Eligibilité:

Modèle 1:

## data$elig: 0
## [1] 0.04597785
## ------------------------------------------------------------ 
## data$elig: 1
## [1] 0.1950073

Modèle 2:

## data$elig: 0
## [1] 0.04602732
## ------------------------------------------------------------ 
## data$elig: 1
## [1] 0.1948495

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..............  6543 
## Original number of treated obs...............  95 
## Matched number of observations...............  95 
## Matched number of observations  (unweighted).  95

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:95))
control <- cbind(control, PAIR = c(1:95))

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..............  6541 
## Original number of treated obs...............  95 
## Matched number of observations...............  95 
## Matched number of observations  (unweighted).  95

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:95))
control <- cbind(control, PAIR = c(1:95))

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.017072, df = 187.99, p-value = 0.9864
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.05607278  0.05511058
## sample estimates:
## mean in group 0 mean in group 1 
##       0.2426414       0.2431225

Modèle 2:

## 
##  Welch Two Sample t-test
## 
## data:  SP2 by elig
## t = -0.067465, df = 187.98, p-value = 0.9463
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.05739923  0.05360295
## sample estimates:
## mean in group 0 mean in group 1 
##       0.2412828       0.2431809

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.0532614	0.0583192	0.0397870
SP mod2	0.0509367	0.0623877	0.0365915

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)    57.429   11.919  23.216 1.45e-06 ***
## period3        -7.636   12.246   0.389    0.533    
## elig          129.715    9.032 206.242  < 2e-16 ***
## period3:elig -134.879    9.951 183.730  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = independence 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)    714.4   90.25
## Number of clusters:   94  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)      66.5    18.0 13.59  0.00023 ***
## period3         -16.8    18.2  0.85  0.35723    
## elig            120.6    15.1 63.50  1.6e-15 ***
## period3:elig   -125.7    15.7 64.46  1.0e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = independence 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)      695    88.1
## Number of clusters:   94  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.8013  0.0396 14717.96   <2e-16 ***
## period3       -0.0747  0.0427     3.06     0.08 .  
## elig          -0.0622  0.1880     0.11     0.74    
## period3:elig   0.0657  0.1892     0.12     0.73    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = independence 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)     1.93   0.227
## Number of clusters:   90  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.7536  0.0812 3429.25   <2e-16 ***
## period3       -0.0311  0.0826    0.14     0.71    
## elig          -0.0145  0.2298    0.00     0.95    
## period3:elig   0.0221  0.2307    0.01     0.92    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = independence 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)     1.89   0.225
## Number of clusters:   91  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	333	366	363	364	365	366	362	364	362	339	343	357	357	353	359
1	0	0	0	1	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	333	366	363	364	365	366	362	364	362	339	343	357	357	353	359
1	0	0	0	1	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	333	366	363	364	365	366	362	364	362	339	343	357	357	353	359
1	0	0	0	1	0	2	1	0	0	3	1	3	27	22	8	8	13	6