R/unequal_unordered_cov.R

In powerLATE: Generalized Power Analysis for LATE

#' @title Subsidiary powerLATE Function
#' @description Subsidiary function to perform power calculation with covariates under unequal assignment probability without ordered mean assumption.
#' @param pZ            probability of being assigned to treatment.
#' @param pi            compliance rate. Equivalently, average causal effect of Z on D.
#' @param N             total number of observations
#' @param kappa         effect size
#' @param sig.level     significance level (Type I error probability).
#' @param power         power of test (1 minus Type II error probability)
#' @param r2dw			proportion of variation in D left unexplained by Z that is explained by W.
#' @param r2yw 			proportion of variation in Y left unexplained by Z that is explained by W.
#' @return A vector of values for one in \{kappa, N, power\} that is not supplied by the user.
#' @note This function is called internally and thus should not be used directly.
#' @author Kirk Bansak and Eddie Yang
#' @seealso \code{\link{equal.unordered.cov}}, \code{\link{equal.ordered.cov}}, \code{\link{unequal.ordered.cov}}.
#' @importFrom stats pnorm qnorm
#' @references Bansak, K. (2020). A Generalized Approach to Power Analysis for Local Average Treatment Effects. Statistical Science, 35(2), 254-271.

unequal.unordered.cov <- function(
	power = NULL,
	sig.level = NULL,
	pi = NULL,
	kappa = NULL,
	N = NULL,
	pZ = NULL,
	r2dw = NULL,
	r2yw = NULL){

	S <- 1-r2yw
	T <- 1-r2dw
	J <- pZ*(1-pZ)

	# mdes
	if (!is.null(power) && !is.null(N)){
		beta <- 1 - power
		M <- qnorm(1-(sig.level/2)) + qnorm(1-beta)
		numer <- 2*M^2*sqrt(S*T) + 4*M*pi*sqrt(N*S*J)
		denom <- 4*N*pi^2*J - M^2*T
		upper.effect <- numer/denom
		return(upper.effect)
	}

	# sample size
	if (!is.null(power) && !is.null(kappa)){
		beta <- 1 - power
		M <- qnorm(1-(sig.level/2)) + qnorm(1-beta)
		numer <- M^2*(kappa^2*T + 4*kappa*sqrt(S*T) + 4*S)
		denom <- 4*kappa^2*pi^2*J
		upper.N <- numer/denom
		return(upper.N)
	}

	# power
	if (!is.null(N) && !is.null(kappa)){
		c.val <- qnorm(1-(sig.level/2)) # 1-alpha/2
		effect.bound <- (kappa*pi*sqrt(J*N))/sqrt(S + 0.25*kappa^2*T + kappa*sqrt(S*T))
		power <- pnorm(-c.val + effect.bound) + pnorm(-c.val - effect.bound)
		return(power)
	}
}