Fleeting

A temporary collection of notes

Quick and Dirty Matching

| Comments

So, you want to compare groups (intervention units) who are similar on baseline but who differ on the intervention. Here is a very very quick and dirty overview of going from data to an estimate of the ATE (taken in part from the documentation for optmatch see also Mark Fredrickson’s work through.)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

library(optmatch)
data(nuclearplants)

## Make a propensity score
ppty <- glm(pr ~ . - (pr + cost), family = binomial(), data = nuclearplants)
ppty.dist <- match_on(ppty) ## a propensity score distance matrix

### Make a mahalanobis (MH) distance matrix (similarity on X) could also rank first.
mhd <- match_on(pr ~ t1 + t2, data = nuclearplants)

## Prepare to match on MH distance but require observations be no more than 
## 2sd on propensity score apart
mhd.pptyc <- mhd + caliper(ppty.dist, width = 2)

pm <- pairmatch(mhd.pptyc,data=nuclearplants) ## for pairs
fm <- fullmatch(mhd.pptyc,data=nuclearplants) ## for a full match 

summary(pm)
summary(fm)

You can then add the matching-factors to the data for later analysis:

1

nuclearplants[names(fm),"fm"]<-as.factor(fm) ## to get rid of attributes and such

For example, you could estimate an ATE with lm:

1
2

lm1 <- lm(cost~pr+fm,data=nuclearplants)
coef(lm1)["pr"]

And you could approximate a randomization standard error for the ATE using an HC2 correction Winston Lin’s 2013 paper:

1
2
3
4

library(sandwich)
library(lmtest)
lm1.hc2<-coeftest(lm1,vcov=vcovHC(lm1,type="HC2"))
lm1.hc2[1:2,]

Here is my hack to get an HC2 CI:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

confint.HC<-function (object, parm, level = 0.95, thevcov, ...) {
  ## a copy of the confint.lm function adding "thevcov" to
    cf <- coef(object)
    pnames <- names(cf)
    if (missing(parm))
        parm <- pnames
    else if (is.numeric(parm))
        parm <- pnames[parm]
    a <- (1 - level)/2
    a <- c(a, 1 - a)
    fac <- qt(a, object$df.residual)
    pct <- stats:::format.perc(a, 3)
    ci <- array(NA, dim = c(length(parm), 2L), dimnames = list(parm,
        pct))
    ##ses <- sqrt(diag(vcov(object)))[parm]
    ses <- sqrt(diag(thevcov))[parm]
    ci[] <- cf[parm] + ses %o% fac
    ci
}


confint.HC(lm1,level=.95,parm="pr",thevcov=vcovHC(lm1,type="HC2"))

Comments