---
title: "Using PanelMatch"
author: "In Song Kim, Adam Rauh, Erik Wang, Kosuke Imai"
date: "`r Sys.Date()`"
output:
  html_document:
    df_print: paged
  pdf_document:
    number_sections: yes
  fig_width: 6 
  fig_height: 4 
vignette: |
  %\VignetteIndexEntry{Using PanelMatch} 
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
---

The goal of this vignette is to provide a quick overview of using the `PanelMatch` package, highlight important features and functions, and help users get the most out of their experience with the package. It assumes that you have been able to install the package successfully and are already familiar with the basics of R.

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

We will be working with the `dem` data set, which comes included with the package. Before doing any matching, estimation or further analysis, it is helpful to understand the distribution of the treatment variable within your data set. The package hopes to facilitate this with the `DisplayTreatment function. 


```{r}
library(PanelMatch)
DisplayTreatment(unit.id = "wbcode2",
                 time.id = "year", legend.position = "none",
                 xlab = "year", ylab = "Country Code",
                 treatment = "dem", data = dem)
```
In the plot, the x axis represents time, and the y axis displays the different units in your data set. Red tiles indicate periods where "treatment" is applied to a given unit and blue tiles indicate "control" periods. White spaces indicate missing data. In the above plot, we have a large number of unique units and time periods, so the axes become hard to read. The `DisplayTreatment` function uses `ggplot2` to create this plot, so any custom styling can be applied easily using `ggplot2` conventions. However, the function has some built in options to make cleaning up the plot a little easier. When the data set is particularly large, it can also help to use the `dense.plot` option. There are many ways to customize these plots. We will return to some of these later, but consult the documentation for the full list and descriptions of the arguments. 


```{r}
DisplayTreatment(unit.id = "wbcode2",
                 time.id = "year", legend.position = "none",
                 xlab = "year", ylab = "Country Code",
                 treatment = "dem", data = dem, 
                 hide.x.axis.label = TRUE, hide.y.axis.label = TRUE) # axis label options
```

```{r}
DisplayTreatment(unit.id = "wbcode2",
                 time.id = "year", legend.position = "none",
                 xlab = "year", ylab = "Country Code",
                 treatment = "dem", data = dem, 
                 hide.x.axis.label = TRUE, hide.y.axis.label = TRUE, 
                 dense.plot = TRUE) # setting dense.plot to TRUE
```

Next, we will move to the `PanelMatch` function. The primary purposes of this function are to 1) create sets matching treated units to control units and 2) determine weights for each control unit in a given matched set. These weights are then used in the estimation stage in an intuitive way: units with higher weights factor more heavily into the estimations.

(1) is achieved by matching units that receive treatment after previously being untreated (ie. units that move from `control` to `treatment` at a certain time) to control units that have matching treatment histories in a specified lag window, while also remaining untreated during the same period that the treated unit receives treatment. For example, if unit 4 is a control unit until 1992, when it receives treatment, then, for a specified lag = 4, it will be matched with control units that share an identical treatment history with unit 4 from 1988-1991, while also *remaining control units* in 1992. 

(2) is achieved by defining which variables should be used for measuring similarity/distance between units, determining the most comparable control units to be included in the matched set, and assigning weights to control units as appropriate. There are many parameters that can be tuned for this step. Users must choose a refinement method ("mahalanobis", "ps.match", "CBPS.match", "ps.weight", "CBPS.weight", "ps.msm.weight", "CBPS.msm.weight", or "none"). The "matching" or mahalanobis refinement methods will assign equal weights to the `size.match` most similar control units in a matched set. The "weighting" methods will generate weights in such a way that control units more similar to treated units will be assigned higher weights. 

Users must also define which covariates should be used in this process for defining similarity between units. This is set using the `covs.formula` argument, which takes the form of a one side formula object. The variables defined on the right hand side of the formula are the variables used in these calculations. Users can included "lagged" versions of variables using `I(lag(name.of.var, 0:n))`.


The first example sets `refinement.method` to `none`, meaning all control units will receive equal weights and no refinement is performed. This will also be helpful to refer back to when we are evaluating the impact of refinement.

The `PanelMatch` function returns a `PanelMatch` object, which contains a `matched.set` object. Please consult the wiki page about [matched set objects for more about them](https://github.com/insongkim/PanelMatch/wiki/Matched-Set-Objects).

```{r}
PM.results.none <- PanelMatch(lag = 4, time.id = "year", unit.id = "wbcode2", 
                         treatment = "dem", refinement.method = "none", 
                         data = dem, match.missing = TRUE, 
                         size.match = 5, qoi = "att", outcome.var = "y",
                         lead = 0:4, forbid.treatment.reversal = FALSE, 
                         use.diagonal.variance.matrix = TRUE)
```


Below, we will use the `mahalanobis` option for `refinement.method` and will use only contemporaneous values of the `tradewb` to define similarity.

```{r}
PM.results <- PanelMatch(lag = 4, time.id = "year", unit.id = "wbcode2", 
                         treatment = "dem", refinement.method = "mahalanobis", # use Mahalanobis distance 
                         data = dem, match.missing = TRUE, 
                         covs.formula = ~ tradewb, 
                         size.match = 5, qoi = "att" , outcome.var = "y",
                         lead = 0:4, forbid.treatment.reversal = FALSE, 
                         use.diagonal.variance.matrix = TRUE)
```

Next, we will include 4 lags of the `tradewb` variable and the outcome variable, excluding any contemporaneous values. 

While there are no hard rules for how to set the many different `PanelMatch` parameters, in general, you want to find a balance between having a good number of matched sets and having matched sets that are large enough in size. Having many small matched sets will lead to larger standard errors. You also want to create sets with good covariate balance, which will be discussed later. 

```{r}
PM.results <- PanelMatch(lag = 4, time.id = "year", unit.id = "wbcode2", 
                         treatment = "dem", refinement.method = "mahalanobis", 
                         data = dem, match.missing = TRUE, 
                         covs.formula = ~ I(lag(tradewb, 1:4)) + I(lag(y, 1:4)), # lags
                         size.match = 5, qoi = "att", outcome.var = "y",
                         lead = 0:4, forbid.treatment.reversal = FALSE, 
                         use.diagonal.variance.matrix = TRUE)
```

We can also apply listwise deletion of units for missing data. 
```{r}
PM.results1 <- PanelMatch(lag = 4, time.id = "year", unit.id = "wbcode2", 
                         treatment = "dem", refinement.method = "mahalanobis", 
                         data = dem, match.missing = FALSE, listwise.delete = TRUE, # listwise deletion used 
                         covs.formula = ~ I(lag(tradewb, 1:4)) + I(lag(y, 1:4)), 
                         size.match = 5, qoi = "att", outcome.var = "y",
                         lead = 0:4, forbid.treatment.reversal = FALSE, 
                         use.diagonal.variance.matrix = TRUE)
```

Let's try out a weighting method using propensity scores and then compare performance.
```{r}
PM.results2 <- PanelMatch(lag = 4, time.id = "year", unit.id = "wbcode2", 
                         treatment = "dem", refinement.method = "ps.weight", 
                         data = dem, match.missing = FALSE, listwise.delete = TRUE, 
                         covs.formula = ~ I(lag(tradewb, 1:4)) + I(lag(y, 1:4)), 
                         size.match = 5, qoi = "att", outcome.var = "y",
                         lead = 0:4, forbid.treatment.reversal = FALSE, 
                         use.diagonal.variance.matrix = TRUE)
```


Now that we've created matched sets, you might be interested in visualizing and evaluating a few aspects of these results. First, let's make sure the results are sensible, by using `DisplayTreatment` to see if our treated and control units meet the parameters outlined previously for one of our matched sets. 

```{r}
# extract the first matched set
mset <- PM.results.none$att[1]

DisplayTreatment(unit.id = "wbcode2",
                 time.id = "year", legend.position = "none",
                 xlab = "year", ylab = "Country Code",
                 treatment = "dem", data = dem,
                 matched.set = mset, # this way we highlight the particular set
                 show.set.only = TRUE)

```

We can also get an idea of the number and size of our matched sets by using the `summary` and `plot` methods for `matched.set` objects. By default, the `plot` method excludes empty matched sets (treated units that could not be matched to any control units) by default, showing the number of empty sets in a separate vertical bar at x = 0. This can be turned off by setting `include.empty.sets` to `TRUE`.

```{r}
summary(PM.results.none$att)

plot(PM.results.none$att)

plot(PM.results.none$att, include.empty.sets = TRUE) # The red tiny bar that would otherwise indicate empty sets is now part of the grey bar
```

Now, let's see the impact of performing refinement by comparing the balance of covariates before and after. See the documentation for the `get_covariate_balance` for more about how the package calculates covariate balance. We hope to see low values (.2 is a helpful threshold value as a rule of thumb) in the results of our covariate balance calculations. We also hope to see an improvement in covariate balance after applying some method of refinement to our results. 


We can use the `get_covariate_balance` function to assess our results. To view the raw results, set `plot = FALSE`. 

```{r}
# get covariate balance for sets that are unrefined

get_covariate_balance(PM.results.none$att,
                      data = dem,
                      covariates = c("tradewb", "y"),
                      plot = FALSE)

# compare with sets that have had various refinements applied

get_covariate_balance(PM.results$att,
                      data = dem,
                      covariates = c("tradewb", "y"),
                      plot = FALSE)

get_covariate_balance(PM.results1$att,
                      data = dem,
                      covariates = c("tradewb", "y"), 
                      plot = FALSE)

get_covariate_balance(PM.results2$att,
                      data = dem,
                      covariates = c("tradewb", "y"), 
                      plot = FALSE)


```

You can also check the balance of unrefined matched sets by setting the `use.equal.weights` argument to `TRUE`. From these examples we can see that different refinement methods have various levels of success in improving covariate balance. Refinement using mahalanobis distance and the `use.diagonal.variance.matrix` argument set to `TRUE` performed well. Refinement using propensity score weighting also helped improve covariate balance, relative to that of the unrefined matched sets. 

The parameters and refinement methods that work best will depend on your data and application/context.

We can also create plots showing covariate balance throughout the lag window period by setting the `plot` argument to `TRUE`, as shown below.
```{r}
get_covariate_balance(PM.results1$att,
                      data = dem,
                      covariates = c("tradewb", "y"), 
                      plot = TRUE, # visualize by setting plot to TRUE
                      ylim = c(-.2, .2))
```

We can also evaluate our results using the `balance_scatter` function:
```{r}
balance_scatter(non_refined_set = PM.results.none$att,
               refined_list = list(PM.results$att, PM.results2$att),
               data = dem,
               covariates = c("y", "tradewb"))
```



We now move to the next major part of the package: obtaining point estimates and standard errors using `PanelEstimate`. The package uses bootstrapping to obtain standard errors. By default, 1000 bootstrap iterations are used, and .95 is the default confidence level. 

The function returns a `PanelEstimate` object, which behaves like a list. As such, you can access the various elements just as you would in a list.

```{r}
PE.results <- PanelEstimate(sets = PM.results1, data = dem)

names(PE.results)

# View the point estimates 

PE.results[["estimates"]]

```

`PanelEsimate` objects have custom `summary` and `plot` methods defined. The `plot` method can be customized with all of the same arguments/operations as the regular `plot` function in base R. This method shows the point estimates for the specified periods, along with the standard errors.


```{r}
summary(PE.results)

plot(PE.results)
```

This last plot makes it clear that the effect of treatment on treated units (att) in this configuration is statistically insignificant.




