J'utilise le package Matching pour faire un matching de donénes sur le score de propension. La fonction Match.out avec mdata permet de récupérer les données appariées avec les variables utilisées dans le matching. J'ai néanmoins deux questions pour lesquelles j'apprécierai beaucoup votre aide :
Quand j'ai une variable factor comme Année ou Secteur ou encore Race dans le fichiers de données citées par le package Matching, il les numérote 1,2,3 et je perds la correspondance avec "race_black","race_hispan","race_white". Comment puis-je faire pour retrouver la correspondance entre les deux dans le tableau X de mdata.
Par ailleurs j'aimerais récupérer les valeurs des autres variables non inclues dans la procédure de matching, comme une variable nom des individus s'il y en avait une. Comment puis-je faire ?
Merci d'avance de votre aide
Voici mon code. L'obtention des données appariées est à la fin.
Code : Tout sélectionner
---
title: "R_matching"
output: html_document
editor_options:
chunk_output_type: console
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
Je suis ce qui est dans le cours "Estimation avec le score de propension sous R" de Simon Quantin pour l'INSEE en 2018 : https://www.insee.fr/fr/statistiques/3546202
https://www.insee.fr/fr/statistiques/fichier/3546202/M2018-01.pdf
```{r}
#install.packages("cobalt", dep=TRUE)
#install.packages("Matching", dep=TRUE)
```
```{r}
library(cobalt)
## Chargement des données : Lalonde Data from cobalt
data(lalonde,package="cobalt")
mybase <- lalonde
## Construction d'indicatrices à partir de la variable
## qualitative RACE
class(mybase$race)
mybase <- splitfactor(data = mybase, var.name = "race",
replace = FALSE,
drop.level=NULL,
drop.first=FALSE)
class(mybase$race)
## Changement d!unité des variables de revenus (en milliers de dollars)
mybase$re74 <- mybase$re74/1000
mybase$re75 <- mybase$re75/1000
## Créations d!indicatrices d!absence de revenus
mybase$u74 <- ifelse(test=(mybase$re74==0),1,0)
mybase$u75 <- ifelse(test=(mybase$re75==0),1,0)
## Création variable factorielle de groupe
mybase$group <- factor(x=mybase$treat,levels=c(0,1),
labels=c("Control","Treated"))
```
Comme le montre les résultats ci-dessous, les moyennes (ou proportions) des différentes cova-
riables diffèrent entre le groupe traité et non traité. Les deux groupes ne sont donc pas en moyenne
directement comparables, ce que l’appariement sur score de propension cherchera à dépasser.
```{r}
var <- names(mybase)[-c(1,4,8,15)]
round(t(sapply(mybase[var], function(x) tapply(x, mybase$treat, mean))), digits = 3)
```
Par ailleurs, nous utiliserons aussi une dataframe composée de deux vecteurs, old et new, per-
mettant d’affecter un libellé aux variables pour les représentations graphiques du package cobalt.
```{r}
## Définition de nouveaux noms de variables
## pour les représentations graphiques
new.names <- data.frame(
old = c("age", "educ", "race_black",
"race_hispan","race_white", "married",
"nodegree", "re74", "re75","u74",
"u75", "pscore"),
new = c("Age", "Years of Education", "Black",
"Hispanic", "White", "Married", "No Degree Earned",
"Earnings 1974", "Earnings 1975", "Employed in 1974",
"Employed in 1975", "Propensity Score"))
```
# 1. Estimation de la propension
L’estimation (ici paramétrique) du score du propension est un préalable commun aux troismé-
thodes présentées. Dans ce document, elle s’appuie sur l’estimation d’un modèle de régression logis-
tique, réalisée avec la fonction glm.
```{r}
ps.logit <- glm(treat ~ re74 + u74 + re75 + u75 + educ + race +
married + nodegree + age + I(age^2) +
nodegree:educ + re74:nodegree + u75:educ,
family = binomial(link="logit"),
data = mybase)
mybase$pscore <- ps.logit$fitted # pour garder le pscore dans la base de données
mybase$pscore
```
# Vérification de la proprité équilibrante
# 2 - Appariement sur le score de propension
```{r}
## Appariement k :1 sur score de propension et matching exact
library(Matching)
Y <- mybase$re78 # variable de revenu
Tr <- mybase$treat # variable de traitement
#Matching avec remise et appariement exact sur black et hisp
Match.out <- Match(Y = Y, Tr = Tr,
estimand="ATT",
M=2, # nombre de voisins
replace = TRUE, # avec remise
ties = TRUE, # prise en compte
# des voisins similaires
caliper = 0.15, # caliper à 0.15 ecart-type
# du score de propension
## déclaration des variables utilisées pour
## l'appariement exact ou non
X = cbind(mybase$pscore,mybase$race),
exact = c("FALSE","TRUE"))
Match.out
```
Sous R, quelle que soit la méthode utilisée (appariement, stratification, ou ajustement par les
pondérations), la fonction bal.tab du package cobalt peut être utilisée pour calculer les différences
de moyennes standardisées et les ratios de variances des covariables afin de vérifier la propriété
équilibrante du score de propension. Sa syntaxe est rappelée ci-dessous :
```{r}
#Vérifier la propriété équilibrante après appariement
##Différences standardisées et ratios de variance avant/après appariement
library(cobalt)
balMatch <- bal.tab(x=Match.out, #et non M= comme ds le pdf
formula = treat ~ re74 + u74 + re75 +
u75 + educ + race + married +
nodegree + age + pscore,
data = mybase,
un = TRUE,
disp.v.ratio = TRUE)
balMatch
## Graphique des différences de moyennes (standardisées) et de proportions
love.plot(x = balMatch,
stat = "mean.diffs",
abs = TRUE,
## options graphiques
var.order = "unadjusted",
var.names = new.names,
threshold = 0.2,
shape = 18)
#pour éviter le message d'erreur mais il faudrait utiliser les label pr que ce soit mieux. en fait ça chnage juste le fait que ce n'est aps la même mesure pour les variables quali et quanti qui est fait et que sinon c'est marqué pareil. là il y a une étoile devant les quali
love.plot(x = balMatch,
stat = "mean.diffs",
abs = TRUE,
## options graphiques
var.order = "unadjusted",
var.names = new.names,
stars = "raw",
threshold = 0.2,
shape = 23)
```
```{r}
#Densité du score de propension avant et après appariement
## Densité du score de propension avant et après appariement
## Les deux variables utilisées pour l!appariement sont déclarées dans formula
bal.plot(x = Match.out,
formula = treat ~ pscore + race,
data = mybase,
var.name = "pscore",
which = "both")
## Densité de l!âge avant et après appariement
## La variable age doit aussi être intégrée dans l!option formula
bal.plot(x = Match.out,
formula = treat ~ age + pscore + race,
data = mybase,
var.name = "age",
which = "both")
bal.plot(x = Match.out,
formula = treat ~ age + pscore + race,
data = mybase,
var.name = "race",
which = "both")
bal.plot(x = Match.out,
formula = treat ~ age + pscore + race,
data = mybase,
var.name = "pscore",
which = "both")
```
#3. Estimation de l'effet du traitement
```{r}
#Estimateur simple après appariement
## 2:1 matching avec remise et appariement exact sur race
Match.out <- Match(Y = Y, Tr = Tr,
X = cbind(mybase$pscore,mybase$race),
estimand = "ATT",
exact = c("FALSE","TRUE"),
M=2,
caliper = 0.15,
sample=TRUE)
summary(Match.out, full=TRUE) #L’option full=TRUE permet
#d’obtenir l’estimateur classique de la variance et celui proposé par Abadie et Imbens (2006) pour tenir
#compte (asymptotiquement) de la variabilité induite par la procédure d’appariement elle-même
```
```{r}
##Estimateur apparié corrigé du biais
# matching avec remise et appariement exact sur race
Match.out <- Match(Y = Y, Tr = Tr,
X = cbind(mybase$pscore,mybase$race),
exact = c("FALSE","TRUE"),
M=2,
caliper = 0.15,
## Correction du biais résiduel
BiasAdjust = TRUE,
## Covariables utilisées pour corriger du biais
Z = cbind(mybase$re74,mybase$re75,mybase$age),
sample=TRUE)
summary(Match.out, full=TRUE)
```
# récupération du fichier de données après appariement
```{r}
mdata <- Match.out[["mdata"]]
X <- mdata[["X"]]
donnees_appariees <- data.frame(mdata[["Y"]],mdata[["Tr"]], X[,1],X[,2])
colnames(donnees_appariees) = c("Y","Tr","pscore","race")
donnees_appariees
```