B.4 Appariement sur le score de propension avec une étape de régression pour réduire les biais de mauvais équilibrage

L’objectif de cette technique est de corriger les biais de mauvais équilibrage qui viennent du fait que lorsque l’on procède à l’appariement dans un échantillon de taille finie, on apparie chaque individu qui a reçu l’appel téléphonique à des individus qui eux n’en ont pas fait l’objet, et qui ont un score de propension très proche mais légèrement différent du premier. Ainsi, les caractéristiques observables pertinentes des seconds, définies par les variables de conditionnement, sont légèrement différentes de celles du premier, ce qui peut lorsque ces différences deviennent trop importantes générer des biais dans l’estimation des effets causaux moyens de l’appel téléphonique incitant à voter sur la participation électorale.

Pour limiter ces biais, Abadie et Imbens (2011) proposent de commencer par estimer deux régressions linéaires séparément, dans lesquelles on régresse la participation électorale \(Y\) sur les variables de conditionnement à l’intérieur de chaque groupe défini par l’intervention : tous les individus qui ont reçu l’appel téléphonique d’un côté, et tous ceux qui ne l’ont pas reçu de l’autre. Cela permet de construire deux jeux de valeurs prédites de la participation électorale \(\hat{Y^1}\) à partir de la régression dans le groupe de ceux qui ont reçu l’appel téléphonique, et \(\hat{Y^0}\) à partir de la régression dans le groupe de ceux qui ne l’ont pas reçu, ces deux jeux de valeurs prédites pouvant être considérés quelque soit le statut réel d’un individu vis-à-vis de l’appel téléphonique incitant à voter.

On peut ensuite appliquer la technique d’appariement. La différence viendra alors du fait que lorsque l’on veut construire les pseudo-effets causaux individuels, au lieu de considérer la différence entre les valeurs observées de la participation électorale \(Y_i-Y_j\) entre un individu \(i\) qui a reçu l’appel téléphonique, et l’individu \(j\) qui ne l’a pas reçu et qu’on lui a apparié, on va plutôt considérer la quantité \(Y_i - \left\{Y_j+\left(\hat{Y_i^0}-\hat{Y_j^0}\right)\right\}\) (et inversement lorsque \(i\) n’a pas reçu l’appel téléphonique). Le terme \(\hat{Y_i^0}-\hat{Y_j^0}\) est la différence entre les valeurs prédites pour \(i\) et pour \(j\) de la régression effectuée exclusivement dans le groupe de ceux qui ont reçu l’appel téléphonique : c’est une fonction linéaire des variables de conditionnement, qui doit donc être nul lorsque \(i\) et \(j\) ont exactement les mêmes valeurs des variables de conditionnement.

Le fragment de code suivant propose d’appliquer cette approche aux données tirées de Gerber et Green (2000) et Imai (2005). Avec cette correction, la valeur estimée des effets causaux moyens de l’appel téléphonique incitant à voter sur la participation électorale atteint 9.1 points de pourcentage. Cette valeur est plus élevée que celle obtenue sans la correction, et plus proche de celles obtenues avec d’autres techniques.

library(data.table)
library(ggplot2)
library(Matching)

#On récupère les données et on convertit en data.table
data(GerberGreenImai)
GerberGreenImai<-data.table(GerberGreenImai)[WARD!="22"]
#pour s'enlever le problème de support commun

#On estime le score de propension avec un modèle logit
pscore_estimates<-glm(PHN.C1 ~ PERSONS + 
                        WARD +
                        AGE +
                        MAJORPTY +
                        VOTE96.0 +
                        VOTE96.1 +
                        NEW + 
                        AGE2 +
                        PERSONS*VOTE96.0 + 
                        PERSONS*NEW,
                      data=GerberGreenImai,
                      family=binomial(link=logit))$fitted.values

#On remet ce score de propension estimé dans la table initiale
GerberGreenImai<-cbind(GerberGreenImai,
                       pscore_estimates)

#On estime les régressions séparément dans chacun des deux groupes
reg_imai_appel<-
  lm(VOTED98 ~
       PERSONS +
       WARD +
       AGE +
       MAJORPTY +
       VOTE96.0 +
       VOTE96.1 +
       AGE2 +
       PERSONS*VOTE96.0 + 
       PERSONS*NEW -
       NEW,
     data=GerberGreenImai[PHN.C1==1],
     x=TRUE)
reg_imai_sans_appel<-
  lm(VOTED98 ~
       PERSONS +
       WARD +
       AGE +
       MAJORPTY +
       VOTE96.0 +
       VOTE96.1 +
       AGE2 +
       PERSONS*VOTE96.0 + 
       PERSONS*NEW -
       NEW,
     data=GerberGreenImai[PHN.C1==0],
     x=TRUE)

#On veut récupèrer le vecteur beta_d de chacun de ces groupes pour calculer
# la valeur prédite par la régression dans toute la populaiton X'beta_d
#Il faut donc aussi construire la matrice X pour toute la population
#Le plus simple pour ça
matrice_population<-
  lm(VOTED98 ~
       PERSONS +
       WARD +
       AGE +
       MAJORPTY +
       VOTE96.0 +
       VOTE96.1 +
       AGE2 +
       PERSONS*VOTE96.0 +
       PERSONS*NEW -
       NEW,
     data=GerberGreenImai,
     x=TRUE)$x

#On n'a plus qu'à calculer les valeurs prédites
valeur_predite_appel<-
  matrice_population%*%as.matrix(reg_imai_appel$coefficients)
valeur_predite_sans_appel<-
  matrice_population%*%as.matrix(reg_imai_sans_appel$coefficients)

#On remet ces valeurs prédites dans la table de départ
GerberGreenImai<-cbind(GerberGreenImai,
                       valeur_predite_appel,
                       valeur_predite_sans_appel)
colnames(GerberGreenImai)[ncol(GerberGreenImai)+(-1:0)]<-
  c("valeur_predite_appel",
    "valeur_predite_sans_appel")

#On commence par apparier chaque individu observé à ses 5 plus proches voisins
# de l'autre groupe, avec remplacement (on n'élimine donc pas au fur et à 
# mesure les individus que l'on a réussi à apparier)
nb_match<-5

#Première étape : la petite table avec les individus de départ, la variable
# qui dit s'ils ont reçu ou non l'appel téléphonique, et la valeur estimée
# du score de propension
GerberGreenImai[,
                indiv_id:=as.character(.I)]
individus_de_depart<-GerberGreenImai[,
                                     c("indiv_id",
                                       "PHN.C1",
                                       "pscore_estimates")]

#On apparie dans un premier temps tous les individus possibles à tous les
# individus de l'autre groupe
individus_a_apparier<-individus_de_depart
individus_a_apparier[,
                     autre_groupe:=1-PHN.C1]
colnames(individus_a_apparier)[1]<-c("indiv_id_a_apparier")
appariement<-merge(individus_de_depart,
                   individus_a_apparier[,
                                        -c("PHN.C1")], 
                   by.x="PHN.C1",
                   by.y="autre_groupe",
                   allow.cartesian = TRUE)[,
                                           -c("autre_groupe")]

#On classe pour chaque individu de départ les couples potentiels dans l'ordre
# croissant de la valeur absolue de la différence des valeurs estimées du 
# score de propension
appariement[,
            abs_diff_pscore:=abs(pscore_estimates.x-
                                   pscore_estimates.y)]
setorder(appariement,
         indiv_id,
         abs_diff_pscore)

#On garde seulement pour chaque individu à apparier les 5 premières 
# propositions : ce sont les individus avec les plus petites différentes de 
# score de propension estimé
appariement[,
            rang_couple:=cumsum(PHN.C1+1-PHN.C1),
            by=indiv_id]
appariement<-appariement[rang_couple<=nb_match,
                         c("indiv_id",
                           "PHN.C1",
                           "indiv_id_a_apparier")]

#Avec appariement
#Pour récupérer les différences moyennes après appariement il faut un peu
# plus de travail
#On crée une table qui stocke (1) les données des individus de départ et (2) les
# données des individus qu'on leur apparie
GerberGreenImai_matching<-rbind(GerberGreenImai,
                                merge(appariement,
                                      GerberGreenImai,
                                      by.x=c("indiv_id_a_apparier"),
                                      by.y="indiv_id"),
                                idcol="Matching",
                                fill=TRUE)
#Pour bien faire la différence entre l'appel reçu par l'individu de départ
# et celui reçu par l'éventuel individu qu'on lui apparie
GerberGreenImai_matching[,
                         c("appel_indiv_dep",
                           "appel_ego"):=list(
                             fcase(
                               Matching==1,
                               PHN.C1,
                               Matching==2,
                               PHN.C1.x),
                             fcase(
                               Matching==1,
                               PHN.C1,
                               Matching==2,
                               PHN.C1.y)
                           )]

#Quand on regarde côté apparié, chaque individu compte 1/5 (le nombre de
# couples que l'on forme pour chaque individu de départ), mais les poids valent
# toujours 1 quand on regarde côté individus de départ
GerberGreenImai_matching[,
                         poids_matching:=fcase(
                           Matching==1,
                           1,
                           Matching==2,
                           1/nb_match
                         )]

#Si on réarrange les termes, on voit que le pseudo-effet individuel que l'on 
# veut estimer a la forme Y_i - \hat{Y_i^0} - (Y_j - \hat{Y_j^0}) quand i
# a reçu l'appel téléphonique (et avec des ^1 dans le cas contraire)
#Il y a donc un intérêt à considérer les résidus, c'est-à-dire les différences
# entre valeur observée et valeur prédite, **dans toute la population** pour
# chacune des régressions estimées dans seulement un des deux groupes
GerberGreenImai_matching[,
                         c("residu_appel",
                           "residu_sans_appel"):=
                           list(VOTED98-valeur_predite_appel,
                                VOTED98-valeur_predite_sans_appel)]

#Le résidu intéressant : pour les individus de départ c'est celui pour
# la régression estimée dans l'autre groupe, alors que pour les individus 
# appariés c'est celui de la régression estimée dans leur groupe
GerberGreenImai_matching[,
                         residu_matching:=
                           fcase(appel_indiv_dep==0,
                                 residu_appel,
                                 appel_indiv_dep==1,
                                 residu_sans_appel)]

#Fonction qui renvoie la différence moyenne entre le groupe qui fait
# l'objet de l'intervention et celui qui n'en fait pas l'objet
diff_moy<-function(variable,
                            groupe_intervention,
                            poids=NULL){
  
  if(is.null(poids)){
    
    (sum(groupe_intervention*variable)/
       sum(groupe_intervention)-
       (sum((1-groupe_intervention)*variable))
     /sum(1-groupe_intervention))
    
  }
  
  else{
    
    (sum(poids*groupe_intervention*variable)/
       sum(poids*groupe_intervention)-
       (sum(poids*(1-groupe_intervention)*variable))
     /sum(poids*(1-groupe_intervention)))
    
  }
  
}

#On n'a plus maintenant qu'à faire la différence entre les individus de
# départ et les individus appariés du point de vue de ces residus
#On commence par estimer la différence de taux de participation électorale
# à l'intérieur de chaque groupe défini par l'intervention, entre les individus
# de départ et ceux auxquels on les a appariés
effet_par_groupe<-
  GerberGreenImai_matching[,
                           lapply(X=.SD,
                                  FUN=function(x){
                                    diff_moy(
                                      groupe_intervention = appel_ego,
                                      variable = x, 
                                      poids=poids_matching)
                                    }),
                           .SDcols="residu_matching",
                           by=c("appel_indiv_dep")]

#On n'a plus qu'à agréger ces effets spécifiques à chaque groupe
effet_par_groupe[,
                 poids_groupe:=
                   fcase(
                     appel_indiv_dep==1,
                     mean(GerberGreenImai$PHN.C1),
                     appel_indiv_dep==0,
                     1-mean(GerberGreenImai$PHN.C1)
                     )]
ATE_matching_corrige<-
  effet_par_groupe[,
                   list(ATE=sum(residu_matching*poids_groupe)/
                          sum(poids_groupe))]
ATE_matching_corrige
##           ATE
## 1: 0.09065878