day13-MatchingInformation.Rmd

---
title: 'Matching Tools\: Propensity scores, Balance.'
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
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
    incremental: true
    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(chemometrics) ## for drawMahal
library(mvtnorm)
library(splines)
library(Rsolnp)
library(parallel)
library(arm) ## See also https://github.com/stan-dev/rstanarm
```

## Today

  1.  Agenda: How to avoid overfit propensity score functions; How to characterize the information in a matched design (use
  what we know about the variance of estimators from block-randomized experiments).
  2. Reading for tomorrow and next week: DOS 8--9, 13 and \cite[\S~9.5]{gelman2006dau}, and \cite{hans:04} (on full matching versus fixed ratio matching versus pair matching), \cite{gelman2008weakly} on "separation" problems in logistic regression other references on the problem of "overfitting" in general.
  3. Updated day-plan.pdf
  4. Questions arising from the reading or assignments or life?

# But first, review:

## What have we done so far?

  - Science is a social process: we persuade ourselves and each other about
    theoretical explanations using logic and observation.
  - Common questions:  "Why should I believe this?", "Why did I do this?",
    "What about another explanation for the observed comparison?"
  - Our answers: (1)  Anticipate alternative explanations; (2)  Develop
    standards (tests should have controlled false positive rate and powerful,
    estimators should be unbiased, precise, if not also consistent, observed
    comparisons should compare favorably  to a randomized comparisons, etc.)

## Decision Points in using algorithmic matching to  create research designs

 - Which covariates and their scaling and coding. (For example, exclude covariates with no variation!)
 - Which distance matrices (scalar distances for one or two important variables, Mahalanobis distances (rank  transformed or not), Propensity distances (using linear predictors)).
 - (Possibly) which calipers (and how many, if any, observations to drop. Note about ATT as a random quantity and ATE/ACE as fixed.)
 - (Possibly) which exact matching or strata
 - (Possibly) which structure of sets (how many treated per control, how many controls per treated)
 - Which remaining differences are  tolerable from a substantive perspective?
 - How well does the resulting research design compare to an equivalent block-randomized study?
 - (Possibly) How much statistical power does this design provide for the quantity of interest?
 - Other questions to ask about a research design aiming to help clarify comparisons.



```{r include=FALSE,   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
```


# More topics: Separation in logistic regression

## The separation problem

```{r echo=TRUE}
thecovs <- unique(c(names(meddat)[c(5:7,9:24)],"HomRate03"))
balfmla<-reformulate(thecovs,response="nhTrt")
psfmla <- update(balfmla,.~.+ns(HomRate03,2)+ns(nhPopD,2)+ns(nhHS,2))
glm0 <- glm(balfmla,data=meddat,family=binomial(link="logit"))
glm1 <- glm(psfmla,data=meddat,family=binomial(link="logit"))
bayesglm0 <- bayesglm(balfmla,data=meddat,family=binomial(link="logit"))
bayesglm1 <- bayesglm(psfmla,data=meddat,family=binomial(link="logit"))
psg1 <- predict(glm1,type="response")
psg0 <- predict(glm0,type="response")
psb1 <- predict(bayesglm1,type="response")
psb0 <- predict(bayesglm0,type="response")
```
## The separation problem

Logistic regression is excellent at discriminating between groups \ldots often **too excellent** for us \autocite{gelman2008weakly}. First evidence of this is big and/or missing coefficients in the propensity score model. See the coefficients below:

```{r echo=FALSE}
thecoefs <- rbind(glm0=coef(glm0)[1:20],
      glm1=coef(glm1)[1:20],
      bayesglm0=coef(bayesglm0)[1:20],
      bayesglm1=coef(bayesglm1)[1:20]
      )
thecoefs[,1:5]
```


## The separation problem


```{r, echo=FALSE, out.width=".9\\textwidth"}
par(mfrow=c(1,2))
matplot(t(thecoefs),axes=FALSE)
axis(2)
axis(1,at=0:19,labels=colnames(thecoefs),las=2)

matplot(t(thecoefs),axes=FALSE,ylim=c(-15,10))
axis(2)
axis(1,at=0:19,labels=colnames(thecoefs),las=2)

legend("topright",col=1:4,lty=1:4,legend=c("glm0","glm1","bayesglm0","bayesglm1"))
```


## The separation problem in logistic regression

```{r out.width=".9\\textwidth", echo=FALSE}
par(mfrow=c(2,2),mar=c(3,3,2,.1))
boxplot(psg0~meddat$nhTrt,main=paste("Logit",length(coef(glm0))," parms",sep=" "))
stripchart(psg0~meddat$nhTrt,vertical=TRUE,add=TRUE)
boxplot(psg1~meddat$nhTrt,main=paste("Logit",length(coef(glm1))," parms",sep=" "))
stripchart(psg1~meddat$nhTrt,vertical=TRUE,add=TRUE)
boxplot(psb0~meddat$nhTrt,main=paste("Shrinkage Logit",length(coef(bayesglm0))," parms",sep=" "))
stripchart(psb0~meddat$nhTrt,vertical=TRUE,add=TRUE)
boxplot(psb1~meddat$nhTrt,main=paste("Shrinkage Logit",length(coef(bayesglm1))," parms",sep=" "))
stripchart(psb1~meddat$nhTrt,vertical=TRUE,add=TRUE)
```

## Matching on the propensity score

```{r}
psdist <- match_on(bayesglm1,data=meddat)
psdist[1:4,1:4]
fmPS <- fullmatch(psdist,data=meddat)
summary(fmPS,min.controls=0,max.controls=Inf,propensity.model=bayesglm1)
```

## Matching on the propensity score

```{r}
stopifnot(all.equal(names(fmPS),row.names(meddat)))
meddat$fmPS <- fmPS ## because same order
meddat$pscore <- predict(bayesglm1) ## because no missing data
setdiffs <- meddat %>% group_by(fmPS) %>% summarize(minHR03=min(HomRate03),
						    maxHR03=max(HomRate03),
						    meanHR03=mean(HomRate03),
						    nset=n(),
						    nT=sum(nhTrt==1),
						    nC=nset - nT)
setdiffs

xbFMPS <- balanceTest(update(balfmla,.~.+pscore + strata(fmPS)),
		      data=meddat,report="all")

xbFMPS
xbFMPS$overall[,]
xbFMPS$results["HomRate03",,]
```

# Balance

## Balance Hunt using the SIUP 

 Hansen and Sales (2015) suggest one way to stop iterating between
 \texttt{fullmatch} and \texttt{xBalance} when you have one caliper. 
 
<!-- The idea
 is that if you would reject the null of balance with one caliper, you would
 also certainly reject it with a wider caliper. That is, the idea is that
 hypothesis tests about balance using calipers can be understood as nested, or
 ordered. Rosenbaum (2008) talks about this in his paper ``Testing Hypotheses
 in Order'' and Hansen and Sales (2008) how these ideas can help us choose a
 matched design: -->

>  ``The SIUP[sequential intersection union principle] states that if a
>  researcher pre-specifies a sequence of hypotheses and corresponding
>  level-$\alpha$ tests, tests those hypotheses in order, and stops testing
>  after the first non-rejected hypothesis, then the probability of incorrectly
>  rejecting at least one correct hypothesis is at most $\alpha$.'' (@hansensales2015, page 2)



## Balance Hunt using the SIUP

Let us try this out and also try to assess it. Say, we start by saying that we will reject the null of balance at $\alpha=.50$.


Imagine, for example we had this matched design:

```{r}
balfmla <- nhTrt ~ nhPopD + nhAboveHS + HomRate03
mhdist <- match_on(balfmla,data=meddat)
psmod <- arm::bayesglm(balfmla,data=meddat,family=binomial(link="logit"))
psdist <- match_on(psmod,data=meddat)
tmp <- meddat$HomRate03
names(tmp) <- rownames(meddat)
absdist <- match_on(tmp, z = meddat$nhTrt,data=meddat)

summary(psdist)
summary(mhdist)
summary(absdist)
```

## Principled Balance Search using SIUP


```{r}
fmMh <- fullmatch(mhdist,data=meddat,tol=.00001)
summary(fmMh,min.controls=0,max.controls=Inf)
```

## Principled Balance Search using SIUP

If we just want to find the set of calipers and optmatch options which maximize
balance, why not do a search?

> The SIUP[sequential intersection union principle] states that if a researcher pre-specifies a sequence of hypotheses and corresponding level-$\alpha$ tests, tests those hypotheses in order, and stops testing after the first non-rejected hypothesis, then the probability of incorrectly rejecting at least one correct hypothesis is at most $\alpha$.'' \autocite{saleshansen2014}(page 2)


```{r}
matchAndBalance<-function(x,balfmla,distmat,thedata){
	#x is a caliper width
	thefm<-fullmatch(distmat+caliper(distmat,x),data=thedata,tol=.00001)
  ## This next is very annoying but there are scope problems with balanceTest and balfmla
  ## And I don't want to add thefm to meddat at each iteration, which would really slow things down.
	thexb<-balanceTest(update(balfmla,.~.+strata(thefm)),
			data=cbind(thedata,thefm),
			report=c("chisquare.test"))
	return(c(x=x,d2p=thexb$overall["thefm","p.value"]))
}
```


## Balance Search using SIUP


```{r cache=FALSE}
## Start with the the largest distance between a treated and control unit.
maxpsdist<-max(as.vector(psdist))
minpsdist<-min(as.vector(psdist))
psdistsum <- summary(psdist)
quantile(as.vector(psdist),seq(0,1,.1))
```

```{r cache=TRUE}
results1<-sapply(seq(3,minpsdist,length=100),function(thecal){
                   matchAndBalance(thecal,balfmla,distmat=psdist,thedata=meddat)})

apply(results1,1,summary)
apply(results1[,results1["d2p",]>.8],1,summary)
```


```{r echo=FALSE}
## Reorder the data from low to high to make cummax work better
results1 <- data.frame(t(results1))
results1o <- results1[order(results1$x,decreasing=TRUE),]
results1o$maxp <- cummax(results1o$d2p)
```

## Balance Search using SIUP

Keeping the maximum produces a set of nested tests:  rejecting balance at some caliper implies that any caliper tighter
than the chosen one would have less balance (a smaller $p$, more information
against the null that our design is like a well randomized block randomized
study).


```{r echo=FALSE, out.width=".8\\textwidth"}
with(results1o,{
       plot(x,d2p,xlab="PS Caliper",ylab="d2 p",cex=.6)
       points(x,maxp,col="blue")
})
```

\note{
  Sometimes we want our matched designs to relate well not only to an
	equivalent block-randomized experiment, but also to help us make the
	argument that our comparisons are comparing specific kinds of like
	with like and/or that our comparisons are statistically powerful.
	That is, among matched designs that we might call "balanced", we might
	one which drops the fewest observations, and perhaps one that has
	specially good balance on certain special covariates (like baseline
	outcomes). So, here is one example, of doing such a search.

	In this case, we are not doing strictly nested hypothesis testing, but are
	using the $p$ values to tell us about information against the null of
	balance rather than using them strictly speaking to reject this null, or
	not-reject it.
  }

## Design Search for both precision and balance

Here I demonstrate searching for two calipers.

```{r gridsearch, cache=FALSE}

findbalance<-function(x){
	##message(paste(x,collapse=" "))
	thefm<-try(fullmatch(psdist+caliper(mhdist,x[2])+caliper(psdist,x[1]),data=meddat,tol=.00001))

	if(inherits(thefm,"try-error")){
		return(c(x=x,d2p=NA,maxHR03diff=NA,n=NA,effn=NA))
	}

	thexb<-try(balanceTest(update(balfmla,.~.+strata(thefm)),
			    data=cbind(meddat,thefm),
			    report=c("chisquare.test","p.values")),silent=TRUE)

	if(inherits(thexb,"try-error")){
		return(c(x=x,d2p=NA,maxHR03diff=NA,n=NA,effn=NA))
	}

	maxHomRate03diff<-max(unlist(matched.distances(thefm,distance=absdist)))

	return(c(x=x,d2p=thexb$overall["thefm","p.value"],
           maxHR03diff=maxHomRate03diff,
           n=sum(!is.na(thefm)),
           effn=summary(thefm)$effective.sample.size))

}

```

## Design Search for both precision and balance

```{r echo=FALSE, cache=TRUE, warning=FALSE}
## Test the function
## findbalance(c(3,3))
## Don't worry about errors for certain combinations of parameters
maxmhdist<-max(as.vector(mhdist))
minmhdist<-min(as.vector(mhdist))
set.seed(123455)
system.time({
	results<-replicate(1000,findbalance(c(runif(1,minpsdist,maxpsdist),
					      runif(1,minmhdist,maxmhdist))))
}
)

```

```{r eval=FALSE, echo=FALSE}
## If you have a mac or linux machine you can speed this up:
system.time({
	resultsList<-mclapply(1:5000,function(i){
				      findbalance(c(runif(1,minpsdist,maxpsdist),
						    runif(1,minmhdist,maxmhdist)))},
			      mc.cores=detectCores())
	resultsListNA<-sapply(resultsList,function(x){ any(is.na(x)) })
	resultsArr<-simplify2array(resultsList[!resultsListNA])
}
)

```


## Which matched design might we prefer?

Now, how might we interpret the results of this search for matched designs?
Here are a few ideas.

```{r }
if(class(results)=="list"){
	resAnyNA<-sapply(results,function(x){ any(is.na(x)) })
	resNoNA<-simplify2array(results[!resAnyNA])
} else {
	resAnyNA<-apply(results,2,function(x){ any(is.na(x)) })
	resNoNA<-simplify2array(results[,!resAnyNA])
}
apply(resNoNA,1,summary)
highbalres<-resNoNA[,resNoNA["d2p",]>.5]
apply(highbalres,1,summary)
```

## Which matched design might we prefer?

```{r eval=TRUE, echo=FALSE}
# color points more dark for smaller differences
plot(resNoNA["d2p",],resNoNA["n",],
     xlab='d2p',ylab='n',
     col=gray(1- ( resNoNA["maxHR03diff",]/max(resNoNA["maxHR03diff",]))),
     pch=19)

## identify(resNoNA["d2p",],resNoNA["n",],labels=round(resNoNA["maxHR03diff",],3),cex=.7)
```

## Which matched design might we prefer?

```{r eval=TRUE,echo=TRUE}
interestingDesigns<- (resNoNA["d2p",]>.5 & resNoNA["n",]>=10 &
		      resNoNA["maxHR03diff",]<=1 & resNoNA["effn",] > 6)
candDesigns <- resNoNA[,interestingDesigns,drop=FALSE]
str(candDesigns)
apply(candDesigns,1,summary)
candDesigns<-candDesigns[,order(candDesigns["d2p",],decreasing=TRUE)]
```

## How would we use this information in `fullmatch`?

```{r bigmatch}
stopifnot(nrow(candDesigns)==1)
fm4<-fullmatch(psdist+caliper(psdist,candDesigns["x1"])+caliper(mhdist,candDesigns["x2"]),data=meddat,tol=.00001)

summary(fm4,min.controls=0,max.controls=Inf)

meddat$fm4<-NULL ## this line exists to prevent confusion with new fm4 objects
meddat[names(fm4),"fm4"]<-fm4

xb3<-balanceTest(update(balfmla,.~.+strata(fm4)),
	      data=meddat, report=c("all"))
xb3$overall[,1:3]
zapsmall(xb3$results["HomRate03",,])
```

## Another approach: more fine tuned optimization

```{r eval=TRUE,cache=FALSE}
matchAndBalance2<-function(x,distmat,alpha){
	#x is a caliper widths
	if(x>max(as.vector(distmat)) | x<min(as.vector(distmat))){ return(99999) }
	thefm<-fullmatch(distmat+caliper(distmat,x),data=meddat,tol=.00001)
	thexb<-xBalance(balfmla,
			strata=data.frame(thefm=thefm),
			data=meddat,
			report=c("chisquare.test"))
	return(thexb$overall[,"p.value"])
}

maxpfn<-function(x,distmat,alpha){
	## here x is the targeted caliper width and x2 is the next wider
	## caliper width
	p1<-matchAndBalance2(x=x[1],distmat,alpha)
	p2<-matchAndBalance2(x=x[2],distmat,alpha)
	return(abs( max(p1,p2) - alpha) )
}

maxpfn(c(minpsdist,minpsdist+1),distmat=psdist,alpha=.25)
#quantile(as.vector(psdist),seq(0,1,.1))
#sort(as.vector(psdist))[1:10]
```

## Another approach: more fine tuned optimization

```{r solnp, warning=FALSE, message=FALSE, cache=TRUE}
### This takes a long time
results3<-gosolnp(fun=maxpfn,
		  ineqfun=function(x,distmat,alpha){ x[2] - x[1] },
		  ineqLB = 0,
		  ineqUB = maxpsdist,
		  LB=c(minpsdist,minpsdist+.01),
		  UB=c(maxpsdist-.01,maxpsdist),
		  n.restarts=2,
		  alpha=.25,distmat=psdist,
		  n.sim=500,
		  rseed=12345,
		  control=list(trace=1)
)
```

## Another approach: more fine tuned optimization

```{r}
maxpfn(results3$pars,distmat=psdist,alpha=.25)
matchAndBalance2(results3$pars[1],distmat=psdist,alpha=.25)
matchAndBalance(results3$par[1],balfmla=balfmla,distmat=psdist,thedata=meddat)
```

## Next:

Estimating causal effects, testing hypotheses about causal effects.



## References