We still use the lalonde data from the MatchIt package and use the propensity score model that we found out in R example 5
library("MatchIt")
data("lalonde")
model <- glm(treat ~ re74 + race + married + I(re74^2) + re74:race, data = lalonde, family = "binomial")
eps <- predict(model, type = "response")
lalonde
Then we perform a first check of the weights by drawing the histogram or boxplot of the weights. At this moment, the figures are not very informative as there are units with extremely large weights.
## Calculate the raw weights (before normalization)
## Weights to estimate ATE
n.treated <- sum(lalonde$treat == 1)
n.control <- sum(lalonde$treat == 0)
weights <- ifelse(lalonde$treat == 1, 1/eps, 1/(1 - eps))
## Check the weights histogram
library(ggplot2)
temp.data <- data.frame(weights = weights, treated = as.factor(lalonde$treat))
ggplot(temp.data, aes(x = weights, fill = treated, color = treated)) +
geom_histogram(alpha = 0.5, position = "identity") +
xlab("Weights")
## We can also draw boxplots
ggplot(temp.data, aes(x = treated, y = weights, color = treated)) +
geom_boxplot()
- Check covariate balancing. If we do not perform any trimming and just look at the full data, we can check covariate balancing from the love plot or summary data of the standardized mean difference. Here, we need to write our own function to draw the love plot.
## Need to change race (categorical) into indicators (numerical)
lalonde$black <- lalonde$race == "black"
lalonde$hispan <- lalonde$race == "hispan"
lalonde$white <- lalonde$race == "white"
## Draw love plot
love.plot = function(cov, treat, ## cov is the matrix of covariates and treat is a vector of treatment assignment
weights = rep(1, length(treat)),
plot = F)
{
## mean with normalized weights \sum w_i x_i / (\sum w_i)
treat.means <- colSums(cov[treat == 1,] * weights[treat == 1])/sum(weights[treat == 1])
treat.var <- colSums(t(t(cov[treat == 1,]) - treat.means)^2 *
weights[treat == 1])/sum(weights[treat == 1])
control.means <- colSums(cov[treat == 0,] * weights[treat == 0])/sum(weights[treat == 0])
control.var <- colSums(t(t(cov[treat == 0,]) - control.means)^2 *
weights[treat == 0])/sum(weights[treat == 0])
## the standardized mean differences for every covariate
smd <- (treat.means - control.means)/sqrt((treat.var + control.var)/2)
names(smd) <- colnames(cov)
if (plot == T) {
plot.data <- data.frame(smd = smd, covariates = names(smd))
range <- max(abs(smd))
ggplot(plot.data) + geom_point(aes(x = as.numeric(smd), y = covariates)) +
geom_vline(xintercept = 0) + xlim(-range, range) +
labs(x = 'Standardized Difference in Means')
}
return(smd)
}
colnames(lalonde)
[1] "treat" "age" "educ" "race" "married" "nodegree" "re74" "re75"
[9] "re78" "black" "hispan" "white"
raw.smd <- love.plot(lalonde[, c(2:3, 5:9, 10:12)], lalonde$treat)
weighted.smd <- love.plot(lalonde[, c(2:3, 5:9, 10:12)], lalonde$treat, weights = weights)
plot.data <- data.frame(smd = c(raw.smd, weighted.smd),
covariates = c(names(raw.smd), names(weighted.smd)),
category = c(rep("Original", length(raw.smd)), rep("IPW", length(weighted.smd))))
range <- max(abs(plot.data$smd))
ggplot(plot.data) + geom_point(aes(x = as.numeric(smd), y = covariates, color = category)) +
geom_vline(xintercept = c(-0.1, -0.05, 0, 0.05, 0.1),
linetype = c("solid", "dashed", "solid", "dashed", "solid")) +
xlim(-range, range) +
labs(x = 'Standardized Difference in Means')
An alternative approach (no love plot) is to use the survey package to create a summary table
library(survey)
## Another simplier way to check balance (no love plot)
weighteddata <- svydesign(ids = ~ 1, data = lalonde, weights = ~weights)
weightedtable <- svyCreateTableOne(vars = colnames(lalonde)[2:8], strata = "treat", data = weighteddata, test = F)
print(weightedtable, smd = T) ## The first row is the summation of weights
Stratified by treat
0 1 SMD
n 615.43 586.79
age (mean (SD)) 27.01 (10.81) 25.46 (6.20) 0.177
educ (mean (SD)) 10.15 (2.84) 10.95 (1.93) 0.328
race (%) 0.043
black 244.6 (39.7) 243.3 (41.5)
hispan 72.0 (11.7) 71.2 (12.1)
white 298.8 (48.6) 272.3 (46.4)
married (mean (SD)) 0.41 (0.49) 0.35 (0.48) 0.110
nodegree (mean (SD)) 0.62 (0.49) 0.52 (0.50) 0.199
re74 (mean (SD)) 4508.71 (6374.58) 3678.23 (4795.67) 0.147
re75 (mean (SD)) 2090.95 (3097.31) 1919.51 (3273.26) 0.054
Now we perform trimming and check covariate balancing again
## Histogram of estimated propensity score
temp.data <- data.frame(eps = eps, treated = as.factor(lalonde$treat))
ggplot(temp.data, aes(x = eps, fill = treated, color = treated)) +
geom_histogram(alpha = 0.5, position = "identity") + xlim(c(0, 1)) +
ggtitle("Histogram of eps before trimming")
rm.idx <- which(eps < 0.1 | eps > 0.9)
## remove control units
length(rm.idx)
[1] 237
## Check the histogram of eps again
ggplot(temp.data[-rm.idx, ], aes(x = eps, fill = treated, color = treated)) +
geom_histogram(alpha = 0.5, position = "identity") + xlim(c(0, 1)) +
ggtitle("Histogram of eps after trimming")
temp.data <- data.frame(weights = weights, treated = as.factor(lalonde$treat))
## We can also draw boxplots
ggplot(temp.data[-rm.idx, ], aes(x = treated, y = weights, color = treated)) +
geom_boxplot()
Check covariate balancing after trimming. Covariate balancing is much better.
raw.smd <- love.plot(lalonde[-rm.idx, c(2:3, 5:9, 10:12)], lalonde$treat[-rm.idx])
weighted.smd <- love.plot(lalonde[-rm.idx, c(2:3, 5:9, 10:12)], lalonde$treat[-rm.idx],
weights = weights[-rm.idx])
plot.data <- data.frame(smd = c(raw.smd, weighted.smd),
covariates = c(names(raw.smd), names(weighted.smd)),
category = c(rep("Original", length(raw.smd)), rep("IPW", length(weighted.smd))))
range <- max(abs(plot.data$smd))
ggplot(plot.data) + geom_point(aes(x = as.numeric(smd), y = covariates, color = category)) +
geom_vline(xintercept = c(-0.1, -0.05, 0, 0.05, 0.1),
linetype = c("solid", "dashed", "solid", "dashed", "solid")) +
xlim(-range, range) +
labs(x = 'Standardized Difference in Means')
- Estimate causal effect by weighted least square regression and estimate the variance by Sandwich estimator
lm.result <- lm(re78 ~ treat, weights = weights[-rm.idx], data = lalonde[-rm.idx, ])
summary(lm.result)
Call:
lm(formula = re78 ~ treat, data = lalonde[-rm.idx, ], weights = weights[-rm.idx])
Weighted Residuals:
Min 1Q Median 3Q Max
-17233 -6644 -2553 4239 64304
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5135.9 480.0 10.700 <2e-16 ***
treat 1203.2 676.4 1.779 0.0761 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 9335 on 375 degrees of freedom
Multiple R-squared: 0.008368, Adjusted R-squared: 0.005723
F-statistic: 3.164 on 1 and 375 DF, p-value: 0.07607
library(sandwich)
SE <- sqrt(diag(vcovHC(lm.result, type = "HC2")))[2]
## get the 95% CI
result <- c(lm.result$coefficients[2], SE, c(tau_hat- 1.96 * SE, tau_hat + 1.96 * SE))
names(result) <- c("est", "sd", "CI_lower", "CI_upper")
result
est sd CI_lower CI_upper
1203.2270 768.3267 -393.7199 2618.1209
LS0tCnRpdGxlOiAiUiBFeGFtcGxlIDY6IEludmVyc2UgcHJvcGVuc2l0eSBzY29yZSB3ZWlnaHRpbmciCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCldlIHN0aWxsIHVzZSB0aGUgbGFsb25kZSBkYXRhIGZyb20gdGhlIE1hdGNoSXQgcGFja2FnZSBhbmQgdXNlIHRoZSBwcm9wZW5zaXR5IHNjb3JlIG1vZGVsIHRoYXQgd2UgZm91bmQgb3V0IGluIFIgZXhhbXBsZSA1CgpgYGB7cn0KbGlicmFyeSgiTWF0Y2hJdCIpCmRhdGEoImxhbG9uZGUiKQoKbW9kZWwgPC0gZ2xtKHRyZWF0IH4gcmU3NCArIHJhY2UgKyBtYXJyaWVkICsgSShyZTc0XjIpICsgcmU3NDpyYWNlLCBkYXRhID0gbGFsb25kZSwgZmFtaWx5ID0gImJpbm9taWFsIikKZXBzIDwtIHByZWRpY3QobW9kZWwsIHR5cGUgPSAicmVzcG9uc2UiKQoKbGFsb25kZQpgYGAKClRoZW4gd2UgcGVyZm9ybSBhIGZpcnN0IGNoZWNrIG9mIHRoZSB3ZWlnaHRzIGJ5IGRyYXdpbmcgdGhlIGhpc3RvZ3JhbSBvciBib3hwbG90IG9mIHRoZSB3ZWlnaHRzLiBBdCB0aGlzIG1vbWVudCwgdGhlIGZpZ3VyZXMgYXJlIG5vdCB2ZXJ5IGluZm9ybWF0aXZlIGFzIHRoZXJlIGFyZSB1bml0cyB3aXRoIGV4dHJlbWVseSBsYXJnZSB3ZWlnaHRzLgpgYGB7cn0KIyMgQ2FsY3VsYXRlIHRoZSByYXcgd2VpZ2h0cyAoYmVmb3JlIG5vcm1hbGl6YXRpb24pCiMjIFdlaWdodHMgdG8gZXN0aW1hdGUgQVRFCm4udHJlYXRlZCA8LSBzdW0obGFsb25kZSR0cmVhdCA9PSAxKQpuLmNvbnRyb2wgPC0gc3VtKGxhbG9uZGUkdHJlYXQgPT0gMCkKd2VpZ2h0cyA8LSBpZmVsc2UobGFsb25kZSR0cmVhdCA9PSAxLCAxL2VwcywgMS8oMSAtIGVwcykpCgojIyBDaGVjayB0aGUgd2VpZ2h0cyBoaXN0b2dyYW0KbGlicmFyeShnZ3Bsb3QyKQp0ZW1wLmRhdGEgPC0gZGF0YS5mcmFtZSh3ZWlnaHRzID0gd2VpZ2h0cywgdHJlYXRlZCA9IGFzLmZhY3RvcihsYWxvbmRlJHRyZWF0KSkKZ2dwbG90KHRlbXAuZGF0YSwgYWVzKHggPSB3ZWlnaHRzLCBmaWxsID0gdHJlYXRlZCwgY29sb3IgPSB0cmVhdGVkKSkgKyAKICBnZW9tX2hpc3RvZ3JhbShhbHBoYSA9IDAuNSwgcG9zaXRpb24gPSAiaWRlbnRpdHkiKSArIAogIHhsYWIoIldlaWdodHMiKSAKCiMjIFdlIGNhbiBhbHNvIGRyYXcgYm94cGxvdHMKZ2dwbG90KHRlbXAuZGF0YSwgYWVzKHggPSB0cmVhdGVkLCB5ID0gd2VpZ2h0cywgY29sb3IgPSB0cmVhdGVkKSkgKyAKICBnZW9tX2JveHBsb3QoKSAKYGBgCgotIENoZWNrIGNvdmFyaWF0ZSBiYWxhbmNpbmcuCklmIHdlIGRvIG5vdCBwZXJmb3JtIGFueSB0cmltbWluZyBhbmQganVzdCBsb29rIGF0IHRoZSBmdWxsIGRhdGEsIHdlIGNhbiBjaGVjayBjb3ZhcmlhdGUgYmFsYW5jaW5nIGZyb20gdGhlIGxvdmUgcGxvdCBvciBzdW1tYXJ5IGRhdGEgb2YgdGhlIHN0YW5kYXJkaXplZCBtZWFuIGRpZmZlcmVuY2UuIEhlcmUsIHdlIG5lZWQgdG8gd3JpdGUgb3VyIG93biBmdW5jdGlvbiB0byBkcmF3IHRoZSBsb3ZlIHBsb3QuIApgYGB7cn0KIyMgTmVlZCB0byBjaGFuZ2UgcmFjZSAoY2F0ZWdvcmljYWwpIGludG8gaW5kaWNhdG9ycyAobnVtZXJpY2FsKQpsYWxvbmRlJGJsYWNrIDwtIGxhbG9uZGUkcmFjZSA9PSAiYmxhY2siCmxhbG9uZGUkaGlzcGFuIDwtIGxhbG9uZGUkcmFjZSA9PSAiaGlzcGFuIgpsYWxvbmRlJHdoaXRlIDwtIGxhbG9uZGUkcmFjZSA9PSAid2hpdGUiCgojIyBEcmF3IGxvdmUgcGxvdApsb3ZlLnBsb3QgPSBmdW5jdGlvbihjb3YsIHRyZWF0LCAgIyMgY292IGlzIHRoZSBtYXRyaXggb2YgY292YXJpYXRlcyBhbmQgdHJlYXQgaXMgYSB2ZWN0b3Igb2YgdHJlYXRtZW50IGFzc2lnbm1lbnQKICAgICAgICAgICAgICAgICAgICAgd2VpZ2h0cyA9IHJlcCgxLCBsZW5ndGgodHJlYXQpKSwKICAgICAgICAgICAgICAgICAgICAgcGxvdCA9IEYpIAp7CiAgICAKICAgICMjIG1lYW4gd2l0aCBub3JtYWxpemVkIHdlaWdodHMgXHN1bSB3X2kgeF9pIC8gKFxzdW0gd19pKQogIHRyZWF0Lm1lYW5zIDwtIGNvbFN1bXMoY292W3RyZWF0ID09IDEsXSAqIHdlaWdodHNbdHJlYXQgPT0gMV0pL3N1bSh3ZWlnaHRzW3RyZWF0ID09IDFdKQogIHRyZWF0LnZhciA8LSBjb2xTdW1zKHQodChjb3ZbdHJlYXQgPT0gMSxdKSAtIHRyZWF0Lm1lYW5zKV4yICoKICAgICAgICAgICAgICAgICAgICAgICAgICB3ZWlnaHRzW3RyZWF0ID09IDFdKS9zdW0od2VpZ2h0c1t0cmVhdCA9PSAxXSkKICAKICBjb250cm9sLm1lYW5zIDwtIGNvbFN1bXMoY292W3RyZWF0ID09IDAsXSAqIHdlaWdodHNbdHJlYXQgPT0gMF0pL3N1bSh3ZWlnaHRzW3RyZWF0ID09IDBdKQogIGNvbnRyb2wudmFyIDwtIGNvbFN1bXModCh0KGNvdlt0cmVhdCA9PSAwLF0pIC0gY29udHJvbC5tZWFucyleMiAqCiAgICAgICAgICAgICAgICAgICAgICAgICAgd2VpZ2h0c1t0cmVhdCA9PSAwXSkvc3VtKHdlaWdodHNbdHJlYXQgPT0gMF0pCiAgCiAgIyMgdGhlIHN0YW5kYXJkaXplZCBtZWFuIGRpZmZlcmVuY2VzIGZvciBldmVyeSBjb3ZhcmlhdGUKICBzbWQgPC0gKHRyZWF0Lm1lYW5zIC0gY29udHJvbC5tZWFucykvc3FydCgodHJlYXQudmFyICsgY29udHJvbC52YXIpLzIpCiAgbmFtZXMoc21kKSA8LSBjb2xuYW1lcyhjb3YpCiAgCiAgaWYgKHBsb3QgPT0gVCkgewogICAgcGxvdC5kYXRhIDwtIGRhdGEuZnJhbWUoc21kID0gc21kLCBjb3ZhcmlhdGVzID0gbmFtZXMoc21kKSkKICAgIHJhbmdlIDwtIG1heChhYnMoc21kKSkKICAgIGdncGxvdChwbG90LmRhdGEpICsgZ2VvbV9wb2ludChhZXMoeCA9IGFzLm51bWVyaWMoc21kKSwgeSA9IGNvdmFyaWF0ZXMpKSArCiAgICAgIGdlb21fdmxpbmUoeGludGVyY2VwdCA9IDApICsgeGxpbSgtcmFuZ2UsIHJhbmdlKSArCiAgICAgIGxhYnMoeCA9ICdTdGFuZGFyZGl6ZWQgRGlmZmVyZW5jZSBpbiBNZWFucycpCiAgfQogIHJldHVybihzbWQpCn0KCmNvbG5hbWVzKGxhbG9uZGUpCnJhdy5zbWQgPC0gbG92ZS5wbG90KGxhbG9uZGVbLCBjKDI6MywgNTo5LCAxMDoxMildLCBsYWxvbmRlJHRyZWF0KQp3ZWlnaHRlZC5zbWQgPC0gbG92ZS5wbG90KGxhbG9uZGVbLCBjKDI6MywgNTo5LCAxMDoxMildLCBsYWxvbmRlJHRyZWF0LCB3ZWlnaHRzID0gd2VpZ2h0cykKCgpwbG90LmRhdGEgPC0gZGF0YS5mcmFtZShzbWQgPSBjKHJhdy5zbWQsIHdlaWdodGVkLnNtZCksIAogICAgICAgICAgICAgICAgICAgICAgICBjb3ZhcmlhdGVzID0gYyhuYW1lcyhyYXcuc21kKSwgbmFtZXMod2VpZ2h0ZWQuc21kKSksCiAgICAgICAgICAgICAgICAgICAgICAgIGNhdGVnb3J5ID0gYyhyZXAoIk9yaWdpbmFsIiwgbGVuZ3RoKHJhdy5zbWQpKSwgcmVwKCJJUFciLCBsZW5ndGgod2VpZ2h0ZWQuc21kKSkpKQpyYW5nZSA8LSBtYXgoYWJzKHBsb3QuZGF0YSRzbWQpKQoKZ2dwbG90KHBsb3QuZGF0YSkgKyBnZW9tX3BvaW50KGFlcyh4ID0gYXMubnVtZXJpYyhzbWQpLCB5ID0gY292YXJpYXRlcywgY29sb3IgPSBjYXRlZ29yeSkpICsKICAgICAgZ2VvbV92bGluZSh4aW50ZXJjZXB0ID0gYygtMC4xLCAtMC4wNSwgMCwgMC4wNSwgMC4xKSwKICAgICAgICAgICAgICAgICBsaW5ldHlwZSA9IGMoInNvbGlkIiwgImRhc2hlZCIsICJzb2xpZCIsICJkYXNoZWQiLCAic29saWQiKSkgKyAKICAgICAgeGxpbSgtcmFuZ2UsIHJhbmdlKSArCiAgICAgIGxhYnMoeCA9ICdTdGFuZGFyZGl6ZWQgRGlmZmVyZW5jZSBpbiBNZWFucycpCmBgYAoKQW4gYWx0ZXJuYXRpdmUgYXBwcm9hY2ggKG5vIGxvdmUgcGxvdCkgaXMgdG8gdXNlIHRoZSBzdXJ2ZXkgcGFja2FnZSB0byBjcmVhdGUgYSBzdW1tYXJ5IHRhYmxlCmBgYHtyfQpsaWJyYXJ5KHN1cnZleSkKIyMgQW5vdGhlciBzaW1wbGllciB3YXkgdG8gY2hlY2sgYmFsYW5jZSAobm8gbG92ZSBwbG90KQp3ZWlnaHRlZGRhdGEgPC0gc3Z5ZGVzaWduKGlkcyA9IH4gMSwgZGF0YSA9IGxhbG9uZGUsIHdlaWdodHMgPSB+d2VpZ2h0cykKd2VpZ2h0ZWR0YWJsZSA8LSBzdnlDcmVhdGVUYWJsZU9uZSh2YXJzID0gY29sbmFtZXMobGFsb25kZSlbMjo4XSwgc3RyYXRhID0gInRyZWF0IiwgZGF0YSA9IHdlaWdodGVkZGF0YSwgdGVzdCA9IEYpCnByaW50KHdlaWdodGVkdGFibGUsIHNtZCA9IFQpICMjIFRoZSBmaXJzdCByb3cgaXMgdGhlIHN1bW1hdGlvbiBvZiB3ZWlnaHRzCmBgYAoKTm93IHdlIHBlcmZvcm0gdHJpbW1pbmcgYW5kIGNoZWNrIGNvdmFyaWF0ZSBiYWxhbmNpbmcgYWdhaW4KYGBge3J9CiMjIEhpc3RvZ3JhbSBvZiBlc3RpbWF0ZWQgcHJvcGVuc2l0eSBzY29yZQp0ZW1wLmRhdGEgPC0gZGF0YS5mcmFtZShlcHMgPSBlcHMsIHRyZWF0ZWQgPSBhcy5mYWN0b3IobGFsb25kZSR0cmVhdCkpCmdncGxvdCh0ZW1wLmRhdGEsIGFlcyh4ID0gZXBzLCBmaWxsID0gdHJlYXRlZCwgY29sb3IgPSB0cmVhdGVkKSkgKyAKICBnZW9tX2hpc3RvZ3JhbShhbHBoYSA9IDAuNSwgcG9zaXRpb24gPSAiaWRlbnRpdHkiKSArIHhsaW0oYygwLCAxKSkgKwogIGdndGl0bGUoIkhpc3RvZ3JhbSBvZiBlcHMgYmVmb3JlIHRyaW1taW5nIikKCnJtLmlkeCA8LSB3aGljaChlcHMgPCAwLjEgfCBlcHMgPiAwLjkpCiMjIHJlbW92ZSBjb250cm9sIHVuaXRzCmxlbmd0aChybS5pZHgpCgojIyBDaGVjayB0aGUgaGlzdG9ncmFtIG9mIGVwcyBhZ2FpbgpnZ3Bsb3QodGVtcC5kYXRhWy1ybS5pZHgsIF0sIGFlcyh4ID0gZXBzLCBmaWxsID0gdHJlYXRlZCwgY29sb3IgPSB0cmVhdGVkKSkgKyAKICBnZW9tX2hpc3RvZ3JhbShhbHBoYSA9IDAuNSwgcG9zaXRpb24gPSAiaWRlbnRpdHkiKSArIHhsaW0oYygwLCAxKSkgKwogIGdndGl0bGUoIkhpc3RvZ3JhbSBvZiBlcHMgYWZ0ZXIgdHJpbW1pbmciKQoKCnRlbXAuZGF0YSA8LSBkYXRhLmZyYW1lKHdlaWdodHMgPSB3ZWlnaHRzLCB0cmVhdGVkID0gYXMuZmFjdG9yKGxhbG9uZGUkdHJlYXQpKQojIyBXZSBjYW4gYWxzbyBkcmF3IGJveHBsb3RzCmdncGxvdCh0ZW1wLmRhdGFbLXJtLmlkeCwgXSwgYWVzKHggPSB0cmVhdGVkLCB5ID0gd2VpZ2h0cywgY29sb3IgPSB0cmVhdGVkKSkgKyAKICBnZW9tX2JveHBsb3QoKSAKYGBgCgpDaGVjayBjb3ZhcmlhdGUgYmFsYW5jaW5nIGFmdGVyIHRyaW1taW5nLiBDb3ZhcmlhdGUgYmFsYW5jaW5nIGlzIG11Y2ggYmV0dGVyLgpgYGB7cn0KcmF3LnNtZCA8LSBsb3ZlLnBsb3QobGFsb25kZVstcm0uaWR4LCBjKDI6MywgNTo5LCAxMDoxMildLCBsYWxvbmRlJHRyZWF0Wy1ybS5pZHhdKQp3ZWlnaHRlZC5zbWQgPC0gbG92ZS5wbG90KGxhbG9uZGVbLXJtLmlkeCwgYygyOjMsIDU6OSwgMTA6MTIpXSwgbGFsb25kZSR0cmVhdFstcm0uaWR4XSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgd2VpZ2h0cyA9IHdlaWdodHNbLXJtLmlkeF0pCnBsb3QuZGF0YSA8LSBkYXRhLmZyYW1lKHNtZCA9IGMocmF3LnNtZCwgd2VpZ2h0ZWQuc21kKSwgCiAgICAgICAgICAgICAgICAgICAgICAgIGNvdmFyaWF0ZXMgPSBjKG5hbWVzKHJhdy5zbWQpLCBuYW1lcyh3ZWlnaHRlZC5zbWQpKSwKICAgICAgICAgICAgICAgICAgICAgICAgY2F0ZWdvcnkgPSBjKHJlcCgiT3JpZ2luYWwiLCBsZW5ndGgocmF3LnNtZCkpLCByZXAoIklQVyIsIGxlbmd0aCh3ZWlnaHRlZC5zbWQpKSkpCnJhbmdlIDwtIG1heChhYnMocGxvdC5kYXRhJHNtZCkpCmdncGxvdChwbG90LmRhdGEpICsgZ2VvbV9wb2ludChhZXMoeCA9IGFzLm51bWVyaWMoc21kKSwgeSA9IGNvdmFyaWF0ZXMsIGNvbG9yID0gY2F0ZWdvcnkpKSArCiAgICAgIGdlb21fdmxpbmUoeGludGVyY2VwdCA9IGMoLTAuMSwgLTAuMDUsIDAsIDAuMDUsIDAuMSksCiAgICAgICAgICAgICAgICAgbGluZXR5cGUgPSBjKCJzb2xpZCIsICJkYXNoZWQiLCAic29saWQiLCAiZGFzaGVkIiwgInNvbGlkIikpICsgCiAgICAgIHhsaW0oLXJhbmdlLCByYW5nZSkgKwogICAgICBsYWJzKHggPSAnU3RhbmRhcmRpemVkIERpZmZlcmVuY2UgaW4gTWVhbnMnKQpgYGAKCi0gRXN0aW1hdGUgY2F1c2FsIGVmZmVjdCBieSB3ZWlnaHRlZCBsZWFzdCBzcXVhcmUgcmVncmVzc2lvbiBhbmQgZXN0aW1hdGUgdGhlIHZhcmlhbmNlIGJ5IFNhbmR3aWNoIGVzdGltYXRvcgoKYGBge3J9CmxtLnJlc3VsdCA8LSBsbShyZTc4IH4gdHJlYXQsIHdlaWdodHMgPSB3ZWlnaHRzWy1ybS5pZHhdLCBkYXRhID0gbGFsb25kZVstcm0uaWR4LCBdKQpzdW1tYXJ5KGxtLnJlc3VsdCkKCmxpYnJhcnkoc2FuZHdpY2gpClNFIDwtIHNxcnQoZGlhZyh2Y292SEMobG0ucmVzdWx0LCB0eXBlID0gIkhDMiIpKSlbMl0KCiMjIGdldCB0aGUgOTUlIENJCnJlc3VsdCA8LSBjKGxtLnJlc3VsdCRjb2VmZmljaWVudHNbMl0sIFNFLCBjKHRhdV9oYXQtIDEuOTYgKiBTRSwgdGF1X2hhdCArIDEuOTYgKiBTRSkpCm5hbWVzKHJlc3VsdCkgPC0gYygiZXN0IiwgInNkIiwgIkNJX2xvd2VyIiwgIkNJX3VwcGVyIikKcmVzdWx0CmBgYAo=