rm(list=ls(all=TRUE))  # clear all variables

# Programmed by Greg Francis (gfrancis@purdue.edu) November 2019

# These functions take various statistics and computes post hoc power (assuming the same sample size as the original test)
# Always assumes a two-tailed test

# This function takes a Wald statistic, which is really a Chi-square 
# This works only for 1 degree of freedom tests
WaldPower=function(W=NULL, alpha=0.05) {
		z = sqrt(W)
		cv = qnorm(1-alpha/2)
		power = pnorm(-cv, mean=z, lower.tail=TRUE) + pnorm(cv, mean=z, lower.tail=FALSE) 
		cat("Power for Wald test=\t", power, "\n")
}

# Various tests based on a normal distribution
zPower=function(Z=NULL, alpha=0.05) {
		cv = qnorm(1-alpha/2)
		power = pnorm(-cv, mean=Z, lower.tail=TRUE) + pnorm(cv, mean=Z, lower.tail=FALSE) 
		cat("Power for Z test=\t", power, "\n")
}

# This function takes a X^2 value (one degree of freedom)
Chi2Power=function(X2=NULL, alpha=0.05) {
		z = sqrt(X2)
		cv = qnorm(1-alpha/2)
		power = pnorm(-cv, mean=z, lower.tail=TRUE) + pnorm(cv, mean=z, lower.tail=FALSE) 
		cat("Power for Chi square test=\t", power, "\n")
}


# This function takes a F value and degrees of freedom
FPower=function(Fvalue=NULL, dfNum=NULL, dfDenom=NULL, alpha=0.05){
	
	# Compute critical value of F
	Fcv = qf(p=alpha, df1=dfNum, df2=dfDenom, lower.tail=FALSE)
	
	# Compute estimate of non-centrality parameter
	lambda = dfNum * Fvalue
	
	# Compute unbiased estimate of non-centrality parameter
	lambdaUnbiased = (dfNum-2)*lambda/dfDenom - dfNum
	
	if(lambdaUnbiased>=0){lambda=lambdaUnbiased}
	
	power = pf(q=Fcv, df1=dfNum, df2=dfDenom, ncp=lambda, lower.tail=FALSE)
	
	cat("Power for F test=\t", power, "\n")
}

# This function takes a t value and degrees of freedom
tPower=function(t=NULL, df=NULL, alpha=0.05) {
	
		# Convert to F test
		Fvalue = t*t
		dfNum=1
		dfDenom = df
		# Compute critical value of F
		Fcv = qf(p=alpha, df1=dfNum, df2=dfDenom, lower.tail=FALSE)
		
		# Compute estimate of non-centrality parameter
		lambda = dfNum * Fvalue
		
		# Compute unbiased estimate of non-centrality parameter
		lambdaUnbiased = (dfNum-2)*lambda/dfDenom - dfNum
		
		if(lambdaUnbiased>=0){lambda=lambdaUnbiased}
		
		power = pf(q=Fcv, df1=dfNum, df2=dfDenom, ncp=lambda, lower.tail=FALSE)
	
		cat("Power for t test=\t", power, "\n")
}