R/bal.ms.psa.R

In PSAgraphics: Propensity Score Analysis Graphics

#' Balance for Continuous Covariate: Random Strata as part of a PSA
#'
#' Function provides a measure (based on the trimmed mean) of the balance
#' achieved between control and treatment groups for a continuous covariate
#' from user defined strata. This statistic is compared to the same measure for
#' randomly permuted strata.
#'
#' This function measures the balance achieved across K strata for a continuous
#' covariate.  If \eqn{ \mu_{ik} } is the covariate trimmed (as specified by
#' user) mean of cases in stratum k, treatment i, then the statistic is the sum
#' over all K of \eqn{ |\mu_{0k} - \mu_{1k}| }.  A permutation distribution is
#' generated by randomly assigning cases to strata, thus generating B permuted
#' stratifications and the associated B permutation statistics.  The
#' permutation stratifications are generated under a fixed marginals model to
#' retain comparability with the original stratification.  A histogram of the
#' permutation statistics is produced with the original statistic referenced as
#' a red dot.
#'
#' @param continuous Quantitative covariate that is being balanced within
#' strata in a PSA.  If \code{continuous} has three columns, then the second
#' and third are assumed to be the treatment and strata respectively.  Missing
#' values are not allowed.
#' @param treatment Binary variable of same length as \code{continuous};
#' generally 0 for 'control,' 1 for 'treatment.'
#' @param strata Integer variable; a vector of same length as \code{continuous}
#' indicating the derived strata from estimated propensity scores.
#' @param trim Fraction (0 to 0.5) of observations to be trimmed from each end
#' of stratum-treatment level before the mean is computed. See
#' \code{\link{mean}}.
#' @param B Numeric; number of randomly generated iterations of the balance
#' measure are created for the comparison distribution.
#' @param main Title passed to \code{histogram}.
#' @return In addition to the histogram, a list with the following components
#' is returned: \item{balance.orig }{Balance measure of user defined strata.}
#' \item{rank.orig}{Rank of original balance measure in comparison with the B
#' randomly generated values.}
#' @author James E. Helmreich \email{ James.Helmreich@@Marist.edu}
#'
#' Robert M. Pruzek \email{RMPruzek@@yahoo.com}
#' @seealso \code{bal.ks.psa}, \code{bal.cws.psa}, \code{bal.cs.psa}
#' @keywords htest
#' @examples
#' #Balance stat should be close to zero
#' meas<-rnorm(500)
#' continuous<-c(meas,meas+rnorm(500,0,.1))
#' treatment<-c(rep(0,500),rep(1,500))
#' strata<-rep(c(rep(1,100),rep(2,100),rep(3,100),rep(4,100),rep(5,100)),2)
#' bal.ms.psa(continuous,treatment,strata)
#'
#'
#' #Balance stat should be close to .4
#' meas<-rnorm(500)
#' continuous<-c(meas, meas[1:250] + runif(250,0,.2),
#'    meas[251:500]-runif(250,0,.2))
#' treatment<-c(rep(0,500),rep(1,500))
#' strata<-rep(c(rep(1,100), rep(2,100), rep(3,100),
#'    rep(4,100),rep(5,100)),2)
#' bal.ms.psa(continuous, treatment, strata, B=200)
#'
#' @export bal.ms.psa
bal.ms.psa <-
	function(continuous,
			 treatment = NULL,
			 strata = NULL,
			 trim = 0,
			 B = 1000,
			 main = NULL) {
		#Compares means within randomly generated strata for a continuous covariate.
		#The analogue of bal.cs.psa


		n <- length(continuous)
		nstrat <- dim(table(strata))

		#If "continuous" has three columns, treat as c, t, s.
		if (dim(as.data.frame(continuous))[2] == 3) {
			treatment   <- continuous[, 2]
			strata      <-
				continuous[, 3]
			continuous <-
				continuous[, 1]
		}


		meas.mns <- tapply(continuous, list(treatment, strata), mean, trim = trim)
		sum.abs.diff.original <-
			sum((abs(meas.mns[1, ] - meas.mns[2, ])) * table(strata) / n)

		sum.abs.diff <- NULL
		for (i in 1:B) {
			rstrat <- sample(strata, n)
			meas.mns <- tapply(continuous, list(treatment, rstrat), mean, trim = trim)
			sum.abs.diff <-
				c(sum.abs.diff, sum((abs(
					meas.mns[1, ] - meas.mns[2, ]
				)) * table(strata) / n))
		}

		res <- c(sum.abs.diff.original, sum.abs.diff)
		rnk <- NULL
		rnk <- rank(c(sum.abs.diff.original[1], sum.abs.diff))[1]
		hist(
			c(sum.abs.diff.original, sum.abs.diff),
			xlab = paste(
				"Balance Statistics for",
				B,
				"Randomly Permuted Stratifications"
			),
			main = main
		)
		points(
			sum.abs.diff.original,
			0,
			col = "red",
			cex = 4,
			pch = 20
		)
		legend(
			x = "topright",
			legend = list(
				paste("Original Balance:",
					  round(sum.abs.diff.original, 2)),
				paste("Rank of Original:", round(rnk, 2), "of", B + 1)
			),
			pch = c(19, 19),
			col = c(2, 0)
		)

		out <- list(balance.orig = sum.abs.diff.original, rank.orig = rnk)
		return(out)
	}