day14-MatchedEstimationTesting.Rmd

---
title: Estimation and Testing for Stratified Matched Designs
date: '`r format(Sys.Date(), "%B %d, %Y")`'
author: ICPSR 2023 Session 1
bibliography:
 - BIB/abbrev-long.bib
 - BIB/refs.bib
 - BIB/master.bib
 - BIB/misc.bib
 - BIB/causalinference.bib
fontsize: 10pt
geometry: margin=1in
graphics: yes
biblio-style: authoryear-comp
biblatexoptions:
  - natbib=true
output:
  beamer_presentation:
    slide_level: 2
    keep_tex: true
    latex_engine: xelatex
    citation_package: biblatex
    template: icpsr.beamer
    includes:
        in_header:
           - defs-all.sty
---

<!-- Make this document using library(rmarkdown); render("day12.Rmd") -->


```{r include=FALSE, cache=FALSE}
# Some customization.  You can alter or delete as desired (if you know what you are doing).
# knitr settings to control how R chunks work.
rm(list = ls())

require(knitr)

## This plus size="\\scriptsize" from https://stackoverflow.com/questions/26372138/beamer-presentation-rstudio-change-font-size-for-chunk

knitr::knit_hooks$set(mysize = function(before, options, envir) {
  if (before) {
    return(options$size)
  } else {
    return("\\normalsize")
  }
})

knit_hooks$set(plotdefault = function(before, options, envir) {
  if (before) par(mar = c(3, 3, .1, .1), oma = rep(0, 4), mgp = c(1.5, .5, 0))
})

opts_chunk$set(
  tidy = "styler", # display code as typed
  echo = TRUE,
  results = "markup",
  strip.white = TRUE,
  fig.path = "figs/fig",
  cache = FALSE,
  highlight = TRUE,
  width.cutoff = 132,
  size = "\\scriptsize",
  out.width = ".7\\textwidth",
  fig.retina = FALSE,
  message = FALSE,
  comment = NA,
  mysize = TRUE,
  plotdefault = TRUE
)

if (!file.exists("figs")) dir.create("figs")

options(
  digits = 4,
  scipen = 8,
  width = 132,
  show.signif.stars = FALSE
)
```

```{r eval=FALSE, include=FALSE, echo=FALSE}
## Run this only once and then not again until we want a new version from github
library("devtools")
library("withr")
with_libpaths("./lib", install_github("markmfredrickson/RItools"), "pre")
```

```{r echo=FALSE}
library(dplyr)
library(ggplot2)
library(RItools, lib.loc = "./lib")
library(optmatch)
library(sandwich)
library(lmtest)
library(estimatr)
library(coin)
library(lme4)
```

## Today


  1. Agenda:  Estimation and Testing given a Matched Stratified Design.
  2. Reading: Review Gerber and Green on Block-Randomized Experiments.
  3. Questions arising from the reading or assignments or life?

# But first, review:


## What have we done so far?

 1. "Interpretable comparison" \citep{kinder1993experimental} versus "No causation
without manipulation." \citep{holland1986}
 2. The problem of overfitting in logistic regression; what this means
for propensity scores; what to do about it.
 3. What does "Effective sample size" mean in the optmatch output? Why should
    we care?
 4. Strategies for choosing a matched design to maximize information and
balance and substantive grounds for argument. (Using your computer to hunt for
a defensible matched design.)
 5. What does it mean to "compare a design to an experiment"? Why do this? How did we do this?


```{r echo=FALSE, cache=TRUE}
load(url("http://jakebowers.org/Data/meddat.rda"))
meddat <- mutate(meddat,
  HomRate03 = (HomCount2003 / Pop2003) * 1000,
  HomRate08 = (HomCount2008 / Pop2008) * 1000
)
row.names(meddat) <- meddat$nh
```

# Estimation and Testing

## Example design for the day

Imagine, that we had this matched design (including matching on missing data):

```{r}
## Make one of the covariates have missing data to
## demonstrate how to match on it
set.seed(12345)
whichmissing <- sample(1:45, 5)
meddat$nhPopD[whichmissing] <- NA
datNoNA <- fill.NAs(nhTrt ~ nhPopD + nhAboveHS + HomRate03,
  data = meddat
)
stopifnot(all.equal(row.names(datNoNA), row.names(meddat)))
datNoNA$id <- meddat$nh
datNoNA$HomRate08 <- meddat$HomRate08
## covs <- c("nhPopD","nhPopD.NA","nhAboveHS","HomRate03")
covs <- unique(c(names(meddat)[c(5:7, 9:24)], "HomRate03"))
balfmla <- reformulate(covs, response = "nhTrt")
mhdist <- match_on(balfmla, data = meddat, method = "rank_mahalanobis")
meddat$nhPopDmd <- datNoNA$nhPopD
meddat$nhPopD.NA <- datNoNA$nhPopD.NA
balfmla <- update(balfmla, . ~ . - nhPopD + nhPopDmd + nhPopD.NA)
mhdist <- match_on(balfmla, data = meddat, method = "rank_mahalanobis")
psmod <- arm::bayesglm(balfmla, data = meddat, family = binomial(link = "logit"))
stopifnot(any(abs(coef(psmod)) < 10))
psdist <- match_on(psmod, data = meddat)
## Make a scalar distance
tmp <- meddat$HomRate03
names(tmp) <- rownames(meddat)
absdist <- match_on(tmp, z = meddat$nhTrt, data = meddat)
```

## Example design for the day

```{r}
## Inspect the distance matrices to choose calipers if desired
quantile(as.vector(psdist), seq(0, 1, .1))
quantile(as.vector(mhdist), seq(0, 1, .1))
quantile(as.vector(absdist), seq(0, 1, .1))
## Match and require no more than 3 treated per control, and no more than 5 control per treated
fmMh <- fullmatch(psdist + caliper(psdist, 5) + caliper(absdist, 2)
  + caliper(mhdist, 52),
min.controls = 0, ## 1/3,
max.controls = Inf,
data = meddat, tol = .00001
)
summary(fmMh, min.controls = 0, max.controls = Inf, propensity.model = psmod)
meddat$fmMh <- factor(fmMh)
meddat$nhTrtF <- factor(meddat$nhTrt)
```

# Testing Hypotheses by Randomization Inference in a Block-Randomized Trial

## Testing Approach: By Hand

```{r}

newexp <- function(trt, b) {
  newtrt <- unsplit(lapply(split(trt, b), sample), b)
  return(newtrt)
}

mdwt1 <- function(y, trt, b) {
  datB <- data.frame(y, trt, b) %>%
    group_by(b) %>%
    summarise(ateb = mean(y[trt == 1]) - mean(y[trt == 0]), nb = n())
  ate_nbwt <- with(datB, sum(ateb * nb / sum(nb)))
  return(ate_nbwt)
}

mdwt2 <- function(y, trt, b) {
  datB <- data.frame(y, trt, b) %>%
    group_by(b) %>%
    summarise(
      ateb = mean(y[trt == 1]) - mean(y[trt == 0]),
      nb = n(),
      nTb = sum(trt),
      nCb = sum(1 - trt),
      pb = mean(trt),
      pbwt = pb * (1 - pb),
      hbwt1 = pbwt * nb,
      hbwt3 = (2 * (nCb * nTb) / (nTb + nCb))
    )
  ate_hbwt <- with(datB, sum(ateb * hbwt1 / sum(hbwt1)))
  return(ate_hbwt)
}
```

## Testing by hand

```{r}
wdat <- meddat %>% filter(!is.na(meddat$fmMh))

obsmd1 <- with(wdat, mdwt1(y = HomRate08, trt = nhTrt, b = fmMh))
obsmd2 <- with(wdat, mdwt2(y = HomRate08, trt = nhTrt, b = fmMh))

origtab <- with(wdat, table(trt = nhTrt, b = fmMh))
testtab <- with(wdat, table(trt = newexp(trt = nhTrt, b = fmMh), b = fmMh))
all.equal(origtab, testtab)
```

## Testing by hand

```{r}

set.seed(12345)
nulldist1 <- replicate(1000, with(wdat, mdwt1(y = HomRate08, trt = newexp(trt = nhTrt, b = fmMh), b = fmMh)))
set.seed(12345)
nulldist2 <- replicate(1000, with(wdat, mdwt2(y = HomRate08, trt = newexp(trt = nhTrt, b = fmMh), b = fmMh)))

p1 <- mean(nulldist1 <= obsmd1)
p2 <- mean(nulldist2 <= obsmd2)

var(nulldist1)
var(nulldist2)
```

```{r}

plot(density(nulldist1))
lines(density(nulldist2), lty = 2)
```



## Testing Approach: Faster

```{r}

xbTest1 <- xBalance(nhTrt ~ HomRate08, strata = list(fmMh = ~fmMh), data = wdat, report = "all")
xbTest1$results[, , "fmMh"]

xbTest2 <- balanceTest(nhTrt ~ HomRate08 + strata(fmMh), data = wdat, report = "all")
xbTest2$results[, , "fmMh"]
```

## Testing Approach: Faster

```{r}
wdat$nhTrtF <- factor(wdat$nhTrt)
meanTestAsym <- oneway_test(HomRate08 ~ nhTrtF | fmMh, data = wdat, distribution = "asymptotic")
set.seed(12345)
meanTestPerm <- oneway_test(HomRate08 ~ nhTrtF | fmMh, data = wdat, distribution = approximate(nresample = 1000))

pvalue(meanTestAsym)
pvalue(meanTestPerm)

rankTestAsym <- wilcox_test(HomRate08 ~ nhTrtF | fmMh, data = wdat, distribution = "asymptotic")
set.seed(12345)
rankTestPerm <- wilcox_test(HomRate08 ~ nhTrtF | fmMh, data = wdat, distribution = approximate(nresample = 1000))

pvalue(rankTestAsym)
pvalue(rankTestPerm)
```

See also `RItest` (in development) and the `ri2` and `ri` packages for R.

# Estimation

## Overview: Estimate and Test "as if block-randomized"

What are we estimating? Most people would say ACE=$\bar{\tau}=\bar{y}_1 - \bar{y}_0$. What estimator estimates this without bias?

```{r echo=TRUE}
meddat[names(fmMh), "fmMh"] <- fmMh
datB <- meddat %>%
  filter(!is.na(fmMh)) %>%
  group_by(fmMh) %>%
  summarise(
    Y = mean(HomRate08[nhTrt == 1]) - mean(HomRate08[nhTrt == 0]),
    nb = n(),
    nbwt = unique(nb / nrow(meddat)),
    nTb = sum(nhTrt),
    nCb = sum(1 - nhTrt),
    pb = mean(nhTrt),
    pbwt = pb * (1 - pb),
    hbwt1 = pbwt * nb,
    hbwt2 = pbwt * nbwt,
    hbwt3 = (2 * (nCb * nTb) / (nTb + nCb))
  )
datB


## Notice that all of these different ways to express the harmonic mean weight are the same.
datB$hbwt101 <- datB$hbwt1 / sum(datB$hbwt1)
datB$hbwt201 <- datB$hbwt2 / sum(datB$hbwt2)
datB$hbwt301 <- datB$hbwt3 / sum(datB$hbwt3)
stopifnot(all.equal(datB$hbwt101, datB$hbwt201))
stopifnot(all.equal(datB$hbwt101, datB$hbwt301))
```

## Using the weights: Set size weights

First, we could estimate the set-size weighted ATE. Our estimator uses the
size of the sets to estimate this quantity.

```{r}
## The set-size weighted version
atewnb <- with(datB, sum(Y * nb / sum(nb)))
atewnb
```

## Using the weights: Set size weights

Sometimes it is convenient to use `lm` because there are R functions for design-based standard errors and confidence intervals.

```{r warnings=FALSE}
meddat$id <- row.names(meddat)
## See Gerber and Green section 4.5 and also Chapter 3 on block randomized experiments. Also Hansen and Bowers 2008.
wdat <- meddat %>%
  filter(!is.na(fmMh)) %>%
  group_by(fmMh) %>%
  mutate(
    pb = mean(nhTrt),
    nbwt = nhTrt / pb + (1 - nhTrt) / (1 - pb),
    gghbwt = 2 * (n() / nrow(meddat)) * (pb * (1 - pb)), ## GG version,
    gghbwt2 = 2 * (nbwt) * (pb * (1 - pb)), ## GG version,
    nb = n(),
    nTb = sum(nhTrt),
    nCb = nb - nTb,
    hbwt1 = (2 * (nCb * nTb) / (nTb + nCb)),
    hbwt2 = nbwt * (pb * (1 - pb))
  )

row.names(wdat) <- wdat$id ## dplyr strips row.names
lm0b <- lm(HomRate08 ~ nhTrt, data = wdat, weight = nbwt)
coef(lm0b)["nhTrt"]
coeftest(lm0b, vcov = vcovHC(lm0b, type = "HC2"))[1:2, ]
```

## Using the weights: Set size weights

Also the linear-model-as-mean-difference-calculate provides convenient, asymptotic approximate confidence intervals.

```{r}
theci0 <- coefci(lm0b, parm = "nhTrt", vcov. = vcovHC(lm0b, type = "HC2"))
theci0
lmE0 <- lm_robust(HomRate08 ~ nhTrt, data = wdat, weight = nbwt)
lmE0
```

## Using the weights: precision weights

Set-size weighting is easy to explain but may differ in terms of precision:

```{r}
atewhb <- with(datB, sum(Y * hbwt1 / sum(hbwt1)))
atewhb
lm1 <- lm_robust(HomRate08 ~ nhTrt + fmMh, data = wdat)
summary(lm1)$coef[2, ]
summary(lmE0)$coef[2, ]
lm1a <- lm_robust(HomRate08 ~ nhTrt, fixed_effects = ~fmMh, data = wdat)
summary(lm1a)$coef[1, ]
```

## Precision weighting

Block-mean centering is another approach.

```{r}
wdat$HomRate08Cent <- with(wdat, HomRate08 - ave(HomRate08, fmMh))
wdat$nhTrtCent <- with(wdat, nhTrt - ave(nhTrt, fmMh))

lm2 <- lm_robust(HomRate08Cent ~ nhTrtCent, data = wdat)
summary(lm2)$coef[2, ]
```


## What about random effects?

```{r}
lmer1 <- lmer(HomRate08 ~ nhTrtF + (1 | fmMh),
  data = wdat,
  verbose = 2, start = 0,
  control = lmerControl(restart_edge = TRUE, optCtrl = list(maxfun = 1000))
)
summary(lmer1)$coef
```


## Which estimator to choose?

The  block-sized weighted approach is unbiased. But unbiased is not the only
indicator quality in an estimator.

##

```{r}
library(DeclareDesign)

thepop <- declare_population(wdat)
theassign <- declare_assignment(blocks = fmMh, block_m_each = table(fmMh, nhTrt))
po_fun <- function(data) {
  data$Y_Z_1 <- data$HomRate08
  data$Y_Z_0 <- data$HomRate08
  data
}
thepo <- declare_potential_outcomes(handler = po_fun)
thereveal <- declare_reveal(Y, Z) ## how does assignment reveal potential outcomes
thedesign <- thepop + theassign + thepo + thereveal

oneexp <- draw_data(thedesign)
## Test
all.equal(origtab, with(oneexp, table(trt = Z, b = fmMh)))
```

```{r}

estimand1 <- declare_estimand(ACE = mean(Y_Z_1 - Y_Z_0))

est1 <- declare_estimator(Y ~ Z,
  estimand = estimand1,
  model = difference_in_means,
  label = "E1: Ignoring Blocks"
)
est2 <- declare_estimator(Y ~ Z,
  fixed_effects = ~fmMh,
  estimand = estimand1, model = lm_robust,
  label = "E2: precision weights fe1"
)
est3 <- declare_estimator(Y ~ Z + fmMh,
  estimand = estimand1, model = lm_robust,
  label = "E3: precision weights fe2"
)

nbwt_est_fun <- function(data) {
  data$newnbwt <- with(data, (Z / pb) + ((1 - Z) / (1 - pb)))
  obj <- lm_robust(Y ~ Z, data = data, weights = newnbwt)
  res <- tidy(obj) %>% filter(term == "Z")
  return(res)
}

hbwt_est_fun <- function(data) {
  data$newnbwt <- with(data, (Z / pb) + ((1 - Z) / (1 - pb)))
  data$newhbwt <- with(data, newnbwt * (pb * (1 - pb)))
  obj <- lm_robust(Y ~ Z, data = data, weights = newhbwt)
  res <- tidy(obj) %>% filter(term == "Z")
  return(res)
}

est4 <- declare_estimator(handler = tidy_estimator(nbwt_est_fun), estimand = estimand1, label = "E4: direct block size weights")
est5 <- declare_estimator(handler = tidy_estimator(hbwt_est_fun), estimand = estimand1, label = "E5: direct precision weights")



direct_demean_fun <- function(data) {
  data$Y <- with(data, Y - ave(Y, fmMh))
  data$Z <- with(data, Z - ave(Z, fmMh))
  obj <- lm_robust(Y ~ Z, data = data)
  data.frame(
    term = "Z",
    estimate = obj$coefficients[[2]],
    std.error = obj$std.error[[2]],
    statistic = obj$statistic[[2]],
    p.value = obj$p.value[[2]],
    conf.low = obj$conf.low[[2]],
    conf.high = obj$conf.high[[2]],
    df = obj$df[[2]],
    outcome = "Y"
  )
}

est6 <- declare_estimator(handler = tidy_estimator(direct_demean_fun), estimand = estimand1, label = "E6: Direct Demeaning")

lmer_est_fun <- function(data) {
  thelmer <- lmer(Y ~ Z + (1 | fmMh),
    data = data,
    control = lmerControl(restart_edge = TRUE, optCtrl = list(maxfun = 1000))
  )
  obj <- summary(thelmer)
  cis <- confint(thelmer, parm = "Z")
  data.frame(
    term = "Z",
    estimate = obj$coefficients[2, 1],
    std.error = obj$coefficients[2, 2],
    statistic = obj$coefficients[2, 3],
    p.value = NA,
    conf.low = cis[1, 1],
    conf.high = cis[1, 2],
    df = NA,
    outcome = "Y"
  )
}


est7 <- declare_estimator(handler = tidy_estimator(lmer_est_fun), estimand = estimand1, label = "E7: Random Effects")
# est1(oneexp)
# est2(oneexp)
est7(oneexp)

thedesign_plus_est <- thedesign + estimand1 + est1 + est2 + est3 + est4 + est5 + est6 + est7
```

## Diagnosands and diagnosis

```{r diagnose, cache=TRUE}

set.seed(12345)

thediagnosands <- declare_diagnosands(
  bias = mean(estimate - estimand),
  rmse = sqrt(mean((estimate - estimand)^2)),
  power = mean(p.value < .05),
  coverage = mean(estimand <= conf.high & estimand >= conf.low),
  mean_estimate = mean(estimate),
  sd_estimate = sd(estimate),
  mean_se = mean(std.error),
  mean_estimand = mean(estimand)
)

diagnosis <- diagnose_design(thedesign_plus_est,
  sims = 10000, bootstrap_sims = 0,
  diagnosands = thediagnosands
)
```

## Results of the Simulation

```{r}
reshape_diagnosis(diagnosis, digits = 4)[, -c(1:2, 4)]

kable(reshape_diagnosis(diagnosis, digits = 4)[, -c(1:2, 4)])
```






## The direct permutation approach

```{r}
library(randomizr)
library(permute)

newExp <- function(z, b) {
  n <- length(z)
  h1 <- how(blocks = b)
  z[shuffle(n, control = h1)]
}

newExp2 <- function(b, nT) {
  block_ra(blocks = b, block_m = nT)
}

newz2 <- newExp2(b = wdat$fmMh, nT = datB$nTb)
testnewExp2 <- sapply(split(newz2, wdat$fmMh), function(x) {
  c(nb = length(x), nTb = sum(x))
})
stopifnot(all.equal(testnewExp2["nb", ], datB$nb, check.attributes = FALSE))
stopifnot(all.equal(testnewExp2["nTb", ], datB$nTb, check.attributes = FALSE))

newz1 <- newExp(wdat$nhTrt, wdat$fmMh)
testnewExp1 <- sapply(split(newz1, wdat$fmMh), function(x) {
  c(nb = length(x), nTb = sum(x))
})
stopifnot(all.equal(testnewExp1["nb", ], datB$nb, check.attributes = FALSE))
stopifnot(all.equal(testnewExp1["nTb", ], datB$nTb, check.attributes = FALSE))


hwtfn <- function(data) {
  tapply(data$Tx.grp, data$stratum.code, function(x) {
    2 * sum((x - mean(x))^2)
  })
}

setsizewtfn <- function(data) {
  tapply(data$Tx.grp, data$stratum.code, function(x) {
    length(x)
  })
}

trtsizewtfn <- function(data) {
  ## Assumes Tx.grp \in \{0,1\} and 1=assigned to treatment
  tapply(data$Tx.grp, data$stratum.code, function(x) {
    sum(x)
  })
}

wtMeanDiffTZ <- function(y, z, b, wtfn) {
  tzb <- mapply(function(yb, zb) {
    mean(yb[zb == 1]) - mean(yb[zb == 0])
  },
  yb = split(y, b),
  zb = split(z, b)
  )
  wts <- wtfn(data.frame(Tx.grp = z, stratum.code = b))
  sum(tzb * wts / sum(wts))
}

## Testing (compare to above)
obsHTZ <- wtMeanDiffTZ(wdat$HomRate08, wdat$nhTrt, wdat$fmMh, hwtfn)
obsNTZ <- wtMeanDiffTZ(wdat$HomRate08, wdat$nhTrt, wdat$fmMh, setsizewtfn)
obsTTZ <- wtMeanDiffTZ(wdat$HomRate08, wdat$nhTrt, wdat$fmMh, trtsizewtfn)
```

##  The direct permutation approach

Test the null hypothesis of no effects:

```{r cache=TRUE}

set.seed(12345)
nulldistHwt <- replicate(10000, wtMeanDiffTZ(wdat$HomRate08, newExp(wdat$nhTrt, wdat$fmMh), wdat$fmMh, hwtfn))
set.seed(12345)
nulldistNwt <- replicate(10000, wtMeanDiffTZ(wdat$HomRate08, newExp(wdat$nhTrt, wdat$fmMh), wdat$fmMh, setsizewtfn))

## Notice more precision with the Harmonic weight in thes p-values.
2 * min(mean(nulldistHwt >= obsHTZ), mean(nulldistHwt <= obsHTZ))
2 * min(mean(nulldistNwt >= obsNTZ), mean(nulldistNwt <= obsNTZ))
```

Comparing the reference distributions to each other and to their Normal approximations.

```{r}
plot(density(nulldistHwt), ylim = c(0, 3))
lines(density(nulldistNwt), lty = 2)
curve(dnorm(x, sd = sd(nulldistHwt)), from = -1, to = 1, col = "gray", add = TRUE)
curve(dnorm(x, sd = sd(nulldistNwt)), from = -1, to = 1, col = "gray", lty = 2, add = TRUE)
```


# Difference in Differences for Matched Designs

## Difference in Differences

Although we have adjusted for contemporaneous differences between neighborhoods and also
adjusted somewhat for time-varying differences within neighborhoods by
matching on baseline outcome, we *might* increase precision and diminish bias by
further adjusting after matching.

```{r}
bal1 <- balanceTest(update(balfmla, . ~ . + strata(fmMh)), data = meddat, report = "all")
t(bal1$results["HomRate03", , ])
```

```{r}
wdat$HDiff <- wdat$HomRate08 - wdat$HomRate03
ddnbwt <- lm(HDiff ~ nhTrt, data = wdat, weights = nbwt)
coef(ddnbwt)[2]
## compare to non-differenced version
atewnb
ddnbwt2 <- lm(HomRate08 ~ nhTrt + HomRate03, data = wdat, weights = nbwt)
coef(ddnbwt2)[2]
```

## Difference in Differences

```{r echo=FALSE}
newdat <- bind_rows(list(
  yr08 = wdat[, c("nhTrt", "fmMh", "HomRate08", "nbwt")],
  yr03 = wdat[, c("nhTrt", "fmMh", "HomRate03", "nbwt")]
),
.id = "year"
)
newdat$post <- as.numeric(newdat$year == "yr08")
newdat$Y <- ifelse(newdat$year == "yr08", newdat$HomRate08, newdat$HomRate03)
newdat$Z <- factor(newdat$nhTrt)
newdat$fm <- factor(newdat$fmMh)

g <- ggplot(data = newdat, aes(x = post, y = Y, color = Z, shape = fm)) +
  scale_shape_manual(values = 1:14) +
  # geom_point() +
  geom_jitter(width = .1) +
  geom_smooth(method = "lm", aes(x = post, y = Y, group = Z, weight = nbwt), se = FALSE)

print(g)
```


## Difference in Differences

```{r}
ddhbwt1 <- lm(HDiff ~ nhTrt + fmMh, data = wdat)
coef(ddhbwt1)[2]
atewhb

## Another method of Harmonic Mean Weighting
wdat <- wdat %>%
  group_by(fmMh) %>%
  mutate(
    HDiffMD = HDiff - mean(HDiff),
    nhTrtMD = nhTrt - mean(nhTrt)
  )
ddhbwt2 <- lm(HDiffMD ~ nhTrtMD, data = wdat)
coef(ddhbwt2)[2]
```
## Difference in Differences

In this case, we don't see big precision improvements (recall that we matched
quite closely on baseline outcome).


```{r eval=FALSE}
cinb <- coefci(ddnbwt, parm = "nhTrt", vcov. = vcov(ddnbwt, type = "HC2"))
cihb1 <- coefci(ddhbwt1, parm = "nhTrt", vcov. = vcov(ddhbwt1, type = "HC2"))
cihb2 <- coefci(ddhbwt2, parm = "nhTrtMD", vcov. = vcov(ddhbwt2, type = "HC2"))
```

Difference in Diffences as Covariance Adjustment

```{r}

dind1 <- lm_lin(HomRate08 ~ nhTrt, covariates = ~HomRate03, weight = nbwt, data = wdat)

dind1

## Versus

lmbasic <- lm_robust(HomRate08 ~ nhTrt, weight = nbwt, data = wdat)
lmbasic

lmfe <- lm_robust(HomRate08 ~ nhTrt, fixed_effects = ~fmMh, data = wdat)
lmfe
```




# In-class time for your work

## Time for your own work

Now that you've seen **a lot**, you can practice.



## Anything Else?

 - Recall the utility of `fill.NAs()`: You can and should match on
   missingness. No reason to throw away observations only because of covariate
   missingness.

 - EXTRA: How would we assess the claim that the sequential intersection union
   principle controls the family-wise error rate for balance tests?


## Summary

Questions? Comments?

## Using the weights: Set size weights

We illustrate an two more elaborate versions of this below (from Winston Lin via the Green Lab SOP).

```{r}
wdat <- wdat %>%
  group_by(fmMh) %>%
  mutate(fmwt = n() / nrow(wdat), nb = n())
X <- model.matrix(~ fmMh - 1, data = wdat)
XminusXbar <- apply(X, 2, function(x) {
  x - mean(x)
})
wrkdat <- cbind.data.frame(wdat, data.frame(XminusXbar))
tmpfmla <- reformulate(grep("fmMh1", names(wrkdat), value = TRUE)[-1], response = "HomRate08")
lmfmla <- update(tmpfmla, . ~ nhTrt * (.))
lmLin <- lm(lmfmla, data = wrkdat)
coef(lmLin)["nhTrt"]
## But, in this case, the HC2 Standard Error is undefined because some of our blocks have too few observations
## Can't calculate var(\bar{y}_T) with only 1 observation.
coeftest(lmLin, vcov = vcovHC(lmLin, type = "HC2"))[1:2, ]

## See Gerber and Green 4.5
wdat$Zf <- factor(wdat$pb)
Z <- model.matrix(~ Zf - 1, data = wdat)
ZminusZbar <- apply(Z, 2, function(z) {
  z - mean(z)
})
wrkdat <- cbind.data.frame(wrkdat, ZminusZbar)
tmpfmla <- reformulate(grep("Zf0", names(wrkdat), value = TRUE)[-1], response = "HomRate08")
lmfmlaZ <- update(tmpfmla, . ~ nhTrt * (.))
lmLinA <- lm(lmfmlaZ, data = wrkdat)
coef(lmLinA)["nhTrt"]
coeftest(lmLinA, vcov = vcovHC(lmLinA, type = "HC2"))[1:2, ]
```

```{r}
lmE2 <- lm_lin(HomRate08 ~ nhTrt, covariates = ~fmMh, data = wdat)
```

## Covariance Adjustment after Matching

@rubin:thom:2000 suggest covariance adjustment after matching.

## References