-
Notifications
You must be signed in to change notification settings - Fork 14
/
optmatch.r
70 lines (49 loc) · 1.92 KB
/
optmatch.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
library(optmatch)
pause <- function(str) {
invisible(readline(paste("\nNext:", str, "<press ret>")))
}
data(nuclearplants)
pause("A propensity score distance")
(dist.glm <- mdist(glm(pr~.-(pr+cost), family=binomial(), data=nuclearplants)))
pause("A Mahalanobis distance")
(dist.mahal <- mdist(pr ~ t1 + t2, data = nuclearplants))
pause("Absolute difference on a scalar-distance")
sdiffs <- function(treatments, controls) {
abs(outer(treatments$t1, controls$t1, `-`))
}
(dist.abs <- mdist(sdiffs, structure.fmla = pr ~ 1, data = nuclearplants))
pause("Pair match")
plantspm <- pairmatch(dist.glm)
stratumStructure(plantspm)
(plantspm.d <- matched.distances(plantspm,dist.glm))
mean(unlist(plantspm.d))
pause("A 1:2 match")
plantstm <- pairmatch(dist.glm,controls=2)
stratumStructure(plantstm)
mean(unlist(matched.distances(plantstm,dist.glm)))
unlist(lapply(matched.distances(plantstm,dist.glm), max))
pause("A full match")
plantsfm <- fullmatch(dist.mahal)
stratumStructure(plantsfm)
mean(unlist(matched.distances(plantsfm,dist.glm)))
pause("Full matching with restrictions")
plantsfm1 <- fullmatch(dist.mahal,
min.controls=2, max.controls=3)
stratumStructure(plantsfm1)
mean(unlist(matched.distances(plantsfm1,dist.mahal)))
all(plantsfm1[sort(names(plantsfm1))]
== fullmatch(
t(dist.mahal$m),
min.controls=(1/3),
max.controls=(1/2)) )
pause("Mahalanobis matching with propensity calipers")
plantsfm2 <- fullmatch(caliper(dist.glm, width = .25) + dist.mahal)
matched.distances(plantsfm2,dist.glm)
matched.distances(plantsfm2,dist.mahal)
pause("Combining matches with original data")
### The usual scenario is that rownames will be in the same
### order as the original data, though it safe to check:
all.equal(names(plantsfm2), row.names(nuclearplants))
cbind(nuclearplants, matches=plantsfm2)
### If the match is not in original order:
cbind(nuclearplants, matches = plantsfm2[row.names(nuclearplants)])