/*
cls
cd "C:/Users/sttp/Documents/nhb/STATA/StataHacks/docs/psmatch2"
do2md using returned_variables_caliper_match, css("../../css/nhb.css") //show
*/
/*
cls
cd "C:\Users\sttp\documents\nhb\stata\StataHacks\docs\Coding Stata\psmatch2"
use http://www.stata-press.com/data/r13/cattaneo2, clear
save cattaneo2, replace
use cattaneo2, clear
*/
/***
---
title: Returned variables from -psmatch2-
...


# Returned variables from -psmatch2-

## Aim

Here, the convenience variables generated by the Stata command **psmatch2** 
are explained using the example dataset *cattaneo2*.

At [stephenporter.org](https://stephenporter.org/understanding-weight-calculations-in-statas-psmatch2/)
an unknown international expert on matching is quoted for:

	I can never figure out what psmatch2 is doing, and that’s actually one of 
	the reasons I don’t use it much and don’t tend to recommend it. 
	The documentation isn’t very clear and there are a number of things where 
	I think it should be doing something but then it’s not clear that it is.

This leads Stephen Porter to conclude:

	Not exactly a compelling endorsement. 
	This is why I now recommend to my students that they avoid psmatch2: 
	how can you justify using a software procedure if you have no idea what it is doing? 
	Especially considering so many other procedures available for Stata, as well as matchit in R.

The main problem for Stephen Porter is the deriving of the weights, *_weight*.

It turns out that -psmatch2- quite elegant transforms trimming and caliper 
matching into simple weights.

From Stata version 13, there has been -teffects- in [TE] and  
[the standard errors for ATT are conservative](https://www.statalist.org/forums/forum/general-stata-discussion/general/2841-standard-errors-for-ate-with-psmatch2) 
Note as written in the link above, that there now is an option **ai(#)** or 
**altvariance** which makes -psmatch2- return the  Abadie-Imbens standard errors.

The beauty of -psmatch2- is that the returned variables makes the user 
a lot freer in the following analysis.
Hence, it is still interesting to use -psmatch2-.


### The -psmatch2- help file

From **psmatch2** version 4.0.12 30jan2016 E. Leuven, B. Sianesi, the help file 
says:

	Convenience variables in Stata’s -psmatch2- (Caliper matching)

	-psmatch2- creates a number of variables for the convenience of the user:

		*_treated* is a variable that equals 0 for control observations and 1 for 
		treatment observations.

		*_support* is an indicator variable with equals 1 if the observation is 
		on the common support and 0 if the observation is off the support.

		*_pscore* is the estimated propensity score or a copy of the one 
		provided by pscore().

		*_outcome_variable* for every treatment observation stores the value of 
		the matched outcome.

		*_weight*. For nearest neighbor matching, it holds the frequency with 
		which the observation is used as a match; 
		with option ties and k-nearest neighbors matching it holds the normalized
		weight; 
		for kernel matching, and llr matching with a weight other than 
		stata's tricube, it stores the overall weight given to the matched 
		observation. 
		When estimating att only *_weight* = 1 for the treated.

		*_id* In the case of one-to-one and nearest-neighbors matching, a new 
		identifier created for all observations.

		*_nk* In the case of one-to-one and nearest-neighbors matching, for every 
		treatment observation, it stores the observation number of the k-th 
		matched control observation. 
		Do not forget to sort by _id if you want to use the observation 
		number (id) of for example the 1st nearest neighbor as in

			sort _id
			g x_of_match = x[_n1]

		*_nn* In the case of nearest-neighbors matching, for every treatment 
		observation, it stores the number of matched control observations.

This quoted version of the -psmatch2- help comes after the Porter post.

It might be that the help for -psmatch2- has been modified after the post.

With respect to the points mentioned in the Porter post, I think that the 
-psmatch2- help is quite clear.

	
## What is returned from -psmatch2-

We now look into what is returned from -psmatch2- using the same data as 
Stephen Porter does.
Propensities are added.

First we look into what is returned in the one-to-one case, then the one-to-many. 
In both cases, estimates and comparison tests are reported.


### The Stephen Porter example data

```
use http://www.stata-press.com/data/r13/cattaneo2, clear
```
***/
//OFF
use cattaneo2, clear
//ON
/***
Creating the propensities
***/
/**/logit mbsmoke prenatal1 fbaby mmarried medu fedu mage fage mrace frace 
/**/predict p
/***
### The returned variables from -psmatch2-, the one to one case
***/
/***
Using the propensities in the **psmatch2** command
***/
/**/summarize p
display `=0.25*r(sd)'
psmatch2 mbsmoke, outcome(bweight) pscore(p) neighbor(1) common ///
	caliper(`=0.25*r(sd)') noreplacement logit
/**/label variable _treated "Treatment"
/***
For the description of the convenience variables only the outcome, the exposure 
and the propensity score is needed.
***/
/**/keep bweight mbsmoke p _*
/**/format p _pscore %6.3f
/***
Five non-treated rows
***/
preserve
keep if !_treated
list in 1/5, noobs
restore
/***
Five treated rows
***/
preserve
keep if _treated
list in 1/5, noobs
restore
/***
#### *_id*

The variable *_id* a new identifier created for all observations.
The values are row numbers.


#### *_nn*

The variable *_nn* counts the number of neighbors found.
Note variable *_nn* can be zero when no neighbors are found.


#### *_n#*

The variable *_n1* is the *_id* for the first neighbor etc.
Note variable *_n#* can be missing when no neighbors are found.


#### *_treated*

This is a copy of the treatment variable
***/
compare _treated mbsmoke
/***
#### *_support*

The non-treated are always having *_support* equal to one.

Depending on choices in eg options **caliper** and **trim** the treated can be 
in common support (*_support* equal one) or not (*_support* equal zero).

Note that *_support* also is zero if no match is found or already taken 
(*_n1* is missing and *_nn* is zero).
***/
list _* if _treated & !_nn, noobs
/***
#### *_pscore*

This is the same as the variable containing the propensity score.
***/
compare _pscore p
/***
#### *_bweight*

When there is only one matched value it is given in the variable _outcome, 
here *_bweight*.
To see how *_bweight* is the matched value of *bweight*, the rows must be sorted 
*_id*. Then the matched value can be found using the variable *_n1*.
***/
/**/sort _id
/**/generate matched_bweight = bweight[_n1]
compare matched_bweight _bweight
/**/drop matched_bweight
/***
#### *_weight*

The treated always get value one if the variable *_support* is one, ie that the 
row is in the common support.
When there is no replacement, the matched non-treated row also get the value 1.
If the matching is with replacement, the value of *_weight* is the number of 
times the non-treated value is matched to a treated value.

The weights sums up to the same number (the number of treated on support) 
for non-treated and treated.
***/
sumat _weight, statistics(sum) rowby(_treated) decimals(0)
/***
#### *_pdif*

The variable *_pdif* is the absolute difference between the propensity scores 
for the treated and non-treated.
***/
/**/sort _id
/**/generate double tmpdif = abs(_pscore - _pscore[_n1])
compare tmpdif _pdif
/**/drop tmpdif
/***
### Reproducing the estimates and test

The ATT can be reproduced using the command -ttest- (compare with ATT estimates 
from -psmatch2- above)
***/
ttest bweight == _bweight if _nn, unpaired
/***
or regression with frequency weights (compare with ATT estimates from -psmatch2- 
above).
***/
regress bweight i._treated [aw=_weight], noheader
/***
### The variable *_weight* from -psmatch2-, the one to many case


The data is (again)
```
use http://www.stata-press.com/data/r13/cattaneo2, clear
```

The scores are found the same way.
***/
//OFF
use cattaneo2, clear
//ON
/**/logit mbsmoke prenatal1 fbaby mmarried medu fedu mage fage mrace frace 
/**/predict p
/***
And -psmatch2- is called similar as above now requring 5 neighbors and trimming 
5 percent in each end.
***/
/**/summarize p
display `=0.25*r(sd)'
psmatch2 mbsmoke, outcome(bweight) pscore(p) neighbor(5) caliper(`=0.25*r(sd)') logit
/**/label variable _treated "Treatment"
/**/keep bweight mbsmoke p _*
/**/format p _pscore %6.3f
/***
The sum of *_weight* is (again) the same for non-treated and treated, the number
of treated in common support.
***/
sumat _weight, statistics(sum) rowby(_treated) decimals(0)
/***
To understand how the weights are calculated we look at the last of the treated.
***/
/**/sort _treated _id
/**/local treated = _id[_N]
/***
Then all the matched id's for the last treated is found and saved in the local 
*untreat*. A variable *slct* is generated to mark the rows where the matched 
id's appear.

```
local untreat ""
generate slct = .
forvalues val = 1/5 {
	local untreated`val' = _n`val'[_N]
	if `untreated`val'' < . {
		local untreat "`untreat', `untreated`val''"
		replace slct = inlist(`untreated`val'', _n1, _n2, _n3, _n4, _n5)
	}
}
```
***/
//OFF
local untreat ""
generate slct = .
forvalues val = 1/5 {
	local untreated`val' = _n`val'[_N]
	if `untreated`val'' < . {
		local untreat "`untreat', `untreated`val''"
		replace slct = inlist(`untreated`val'', _n1, _n2, _n3, _n4, _n5)
	}
}
//ON
/****/display `"We see that the last treated with id `treated' is matched with the non-treated id's (`=usubstr("`untreat'", 3,.)')."'
/****/display `"`list bweight _* if inlist(_id, `treated' `untreat'), noobs`
/***/list bweight _* if inlist(_id, `treated' `untreat'), noobs
/***
The value of *_bweight* for the treated id is the mean of the non-treated values.
***/
/****/display `"`summarize bweight if inlist(_id `untreat')`
/***/summarize bweight if inlist(_id `untreat')
/***
The *_weight* values for the non-treated are:
***/
/**/mkmat _weight if inlist(_id `untreat'), rownames(_id )
/**/matrix roweq _weight = "non-treated id's" 
/****/matprint _weight, style(md) decimals(4)
/***
To understand the weights for the non-treated one must look where they are matched.
***/
list bweight _* if slct, noobs
/***
Eg, the non-treated with id 3775 appear in the means for the treated with 
id's 4640 (n = 5), 4641 (n = 5) and 4642 (n = 3).

So the weight for the non-treated with id 3775 is:

$$\frac{1}{5} + \frac{1}{5} + \frac{1}{3} = \frac{11}{15} = 0.7333$$

The weight for the non-treated with id 3776 is:

$$\frac{1}{5} + \frac{1}{3} = \frac{8}{15} = 0.5333$$

Finally, the weight for the non-treated with id 3777 is:

$$\frac{1}{3} = 0.3333$$

Quite simple.


### Reproducing the estimates and test

Again, the ATT can be reproduced using the command -ttest- 
(compare with ATT estimates from -psmatch2- above)
***/
ttest bweight == _bweight if _nn, unpaired
/***
or regression with frequency weights (compare with ATT estimates from -psmatch2- 
above).
***/
regress bweight i._treated [aw=_weight], noheader