################################################
# Functions for estimating a confidence interval on standardized effect sizes.
# Load the function by running the function first. 
# Several examples are presented at the end
# See Biesanz, J. C. (2009).  Constructing Confidence Intervals for Standardized Effect Sizes. Manuscript under review.
# for references and discussion.
#################################################

################################################
#CIR.r:  Function for estimating the confidence interval for an observed correlation or partial correlation
#Assumes that variables are randomly sampled as opposed to a fixed predictor
#Variable definitions:
#  r:  		observed correlation or partial correlation
#  df:  	degrees of freedom from overall analysis
#  conf:	desired confidence level (two-sided)
#  iter:	how many iterations of the stochastic approximation to run.  Over 1000 may be quite slow.
#Returns:
#  1.  Estimates of the upper and lower confidence interval endpoints.
################################################

CIR.r<- function(r, df, conf, iter){
	M <-200000
	alphaL <- (1-conf)/2
	alphaU <- 1- alphaL
	hr <- (rnorm(M,mean=0,sd=1) +r*sqrt(rchisq(n= M,df=df,ncp=0))/(sqrt(1-r*r)) )/sqrt(rchisq(n= M,df=df+1,ncp=0)) 
	CID <- hr/sqrt(hr*hr+1)
	
	initial <- quantile(CID,probs=c(alphaL, alphaU))
	initialqL <- initial[1]
	initialqU <- initial[2]
	lastqL <- initialqL
	lastqU <- initialqU

	lastdenL <- density(CID, from= lastqL, to = lastqL +1)$y[1]
	initialdenL <- lastdenL
	lastdenU <- density(CID, from= lastqU, to = lastqU +1)$y[1]
	initialdenU <- lastdenU
	
	for (n in 1: iter){	
		hr <- (rnorm(M,mean=0,sd=1) +r*sqrt(rchisq(n= M,df=df,ncp=0))/(sqrt(1-r*r)) )/sqrt(rchisq(n= M,df=df+1,ncp=0)) 
		CID <- hr/sqrt(hr*hr+1)

		fnL <- (1-1/n)*lastdenL + (1/n)*(sum(abs(CID-lastqL)<=(1/sqrt(n))))/(2*M*(1/sqrt(n)))
		SnL <- lastqL + (1/(n*max(fnL,initialdenL/sqrt(n))))*(alphaL - sum(CID <= lastqL)/M)
		
		fnU <- (1-1/n)*lastdenU + (1/n)*(sum(abs(CID-lastqU)<=(1/sqrt(n))))/(2*M*(1/sqrt(n)))
		SnU <- lastqU + (1/(n*max(fnU,initialdenU/sqrt(n))))*(alphaU - sum(CID <= lastqU)/M)
		lastqL <- SnL
		lastdenL <- fnL
		lastqU <- SnU
		lastdenU <- fnU
	}
	c(SnL,SnU)
}

################################################
#CIF.d:  Function for estimating the confidence interval for an observed standardized mean difference
#Assumes that the predictor is fixed.
#Variable definitions:
#  d:  		observed standardized mean difference
#  n1		sample size for group 1 
#  n2		sample size for group 2 
#  conf:	desired confidence level (two-sided)
#  iter:	how many iterations of the stochastic approximation to run.  Over 1000 may be quite slow.
#Returns:
#  1.  Estimates of the upper and lower confidence interval endpoints.
################################################

CIF.d<- function(d,n1,n2,conf, iter){
	M <-200000
	alphaL <- (1-conf)/2
	alphaU <- 1- alphaL
	df <- n1+n2-2
	CID <- rnorm(M,mean=0,sd=1)*sqrt(1/n1 + 1/n2) +d*sqrt(rchisq(n= M,df=df,ncp=0)/(df)) 
	initial <- quantile(CID,probs=c(alphaL, alphaU))
	initialqL <- initial[1]
	initialqU <- initial[2]
	lastqL <- initialqL
	lastqU <- initialqU

	lastdenL <- density(CID, from= lastqL, to = lastqL +1)$y[1]
	initialdenL <- lastdenL
	lastdenU <- density(CID, from= lastqU, to = lastqU +1)$y[1]
	initialdenU <- lastdenU
	
	for (n in 1: iter){	
		CID <- rnorm(M,mean=0,sd=1)*sqrt(1/n1 + 1/n2) +d*sqrt(rchisq(n= M,df=df,ncp=0)/(df)) 

		fnL <- (1-1/n)*lastdenL + (1/n)*(sum(abs(CID-lastqL)<=(1/sqrt(n))))/(2*M*(1/sqrt(n)))
		SnL <- lastqL + (1/(n*max(fnL,initialdenL/sqrt(n))))*(alphaL - sum(CID <= lastqL)/M)
		
		fnU <- (1-1/n)*lastdenU + (1/n)*(sum(abs(CID-lastqU)<=(1/sqrt(n))))/(2*M*(1/sqrt(n)))
		SnU <- lastqU + (1/(n*max(fnU,initialdenU/sqrt(n))))*(alphaU - sum(CID <= lastqU)/M)
		lastqL <- SnL
		lastdenL <- fnL
		lastqU <- SnU
		lastdenU <- fnU
	}
	c(SnL,SnU)
}

################################################
#CIR.d:  Function for estimating the confidence interval for an observed standardized mean difference
#Assumes that the predictor is randomly sampled.
#Variable definitions:
#  d:  		observed standardized mean difference
#  n1		sample size for group 1 
#  n2		sample size for group 2 
#  conf:	desired confidence level (two-sided)
#  iter:	how many iterations of the stochastic approximation to run.  Over 1000 may be quite slow.
#Returns:
#  1.  Estimates of the upper and lower confidence interval endpoints.
################################################

CIR.d<- function(d,n1,n2,conf){
	M <-200000
	alphaL <- (1-conf)/2
	alphaU <- 1- alphaL
	tobs<- d/sqrt(1/n1 + 1/n2)
	df <- n1+n2-2
	CID <- sqrt(1/n1 + 1/n2)*(rnorm(M,mean=0,sd=1)*sqrt(df) +tobs*sqrt(rchisq(n= M,df=df,ncp=0)))/sqrt(rchisq(n= M,df=(df+1),ncp=0)) 
	initial <- quantile(CID,probs=c(alphaL, alphaU))
	initialqL <- initial[1]
	initialqU <- initial[2]
	lastqL <- initialqL
	lastqU <- initialqU

	lastdenL <- density(CID, from= lastqL, to = lastqL +1)$y[1]
	initialdenL <- lastdenL
	lastdenU <- density(CID, from= lastqU, to = lastqU +1)$y[1]
	initialdenU <- lastdenU
	
	for (n in 1: iter){	
		CID <- sqrt(1/n1 + 1/n2)*(rnorm(M,mean=0,sd=1)*sqrt(df) +tobs*sqrt(rchisq(n= M,df=df,ncp=0)))/sqrt(rchisq(n= M,df=(df+1),ncp=0)) 

		fnL <- (1-1/n)*lastdenL + (1/n)*(sum(abs(CID-lastqL)<=(1/sqrt(n))))/(2*M*(1/sqrt(n)))
		SnL <- lastqL + (1/(n*max(fnL,initialdenL/sqrt(n))))*(alphaL - sum(CID <= lastqL)/M)
		
		fnU <- (1-1/n)*lastdenU + (1/n)*(sum(abs(CID-lastqU)<=(1/sqrt(n))))/(2*M*(1/sqrt(n)))
		SnU <- lastqU + (1/(n*max(fnU,initialdenU/sqrt(n))))*(alphaU - sum(CID <= lastqU)/M)
		lastqL <- SnL
		lastdenL <- fnL
		lastqU <- SnU
		lastdenU <- fnU
	}
	c(SnL,SnU)
}


################################################
#CIF.beta:  Function for estimating the confidence interval for an observed standardized regression coefficient
#Assumes that the predictor is fixed.
#Variable definitions:
#  beta:  	observed standardized regression coefficient
#  tobs		t-test for the regression coefficient
#  df		df for the t-test
#  R2full	R-squared from the overall analysis
#  R2X.Xk   R-squared when the predictor of interest is predicted by the other variables in the model 
#  conf:	desired confidence level (two-sided)
#  iter:	how many iterations of the stochastic approximation to run.  Over 1000 may be quite slow.
#Returns:
#  1.  Estimates of the upper and lower confidence interval endpoints.
################################################	
CIF.beta<- function(beta,tobs,df,R2full,R2X.Xk,conf, iter){
	M <-200000
	alphaL <- (1-conf)/2
	alphaU <- 1- alphaL
	R2change <- (tobs*tobs)*(1-R2full)/df
	R2reduced <- R2full - R2change
	D <- sqrt((1-R2reduced)/(1-R2X.Xk))
	Af <- rnorm(n=M,mean=0,sd=1)/sqrt(df) + tobs*sqrt(rchisq(n= M,df=(df),ncp=0))/df
	
	CID <- Af*D/sqrt(1+Af*Af)
	initial <- quantile(CID,probs=c(alphaL, alphaU))
	initialqL <- initial[1]
	initialqU <- initial[2]
	lastqL <- initialqL
	lastqU <- initialqU

	lastdenL <- density(CID, from= lastqL, to = lastqL +1)$y[1]
	initialdenL <- lastdenL
	lastdenU <- density(CID, from= lastqU, to = lastqU +1)$y[1]
	initialdenU <- lastdenU
	
	for (n in 1: iter){	
		Af <- rnorm(n=M,mean=0,sd=1)/sqrt(df) + tobs*sqrt(rchisq(n= M,df=(df),ncp=0))/df
		CID <- Af*D/sqrt(1+Af*Af)

		fnL <- (1-1/n)*lastdenL + (1/n)*(sum(abs(CID-lastqL)<=(1/sqrt(n))))/(2*M*(1/sqrt(n)))
		SnL <- lastqL + (1/(n*max(fnL,initialdenL/sqrt(n))))*(alphaL - sum(CID <= lastqL)/M)
		
		fnU <- (1-1/n)*lastdenU + (1/n)*(sum(abs(CID-lastqU)<=(1/sqrt(n))))/(2*M*(1/sqrt(n)))
		SnU <- lastqU + (1/(n*max(fnU,initialdenU/sqrt(n))))*(alphaU - sum(CID <= lastqU)/M)
		lastqL <- SnL
		lastdenL <- fnL
		lastqU <- SnU
		lastdenU <- fnU
	}
	c(SnL,SnU)
}


################################################
#CIR.beta:  Function for estimating the confidence interval for an observed standardized regression coefficient
#Assumes that the predictor is randomly sampled.
#Variable definitions:
#  beta:  	observed standardized regression coefficient
#  tobs		t-test for the regression coefficient
#  df		df for the t-test
#  R2full	R-squared from the overall analysis
#  R2X.Xk   R-squared when the predictor of interest is predicted by the other variables in the model 
#  conf:	desired confidence level (two-sided)
#  iter:	how many iterations of the stochastic approximation to run.  Over 1000 may be quite slow.
#Returns:
#  1.  Estimates of the upper and lower confidence interval endpoints.
################################################
CIR.beta<- function(beta,tobs,df,R2full,R2X.Xk,conf, iter){
	M <-200000
	alphaL <- (1-conf)/2
	alphaU <- 1- alphaL
	R2change <- (tobs*tobs)*(1-R2full)/df
	R2reduced <- R2full - R2change
	D <- sqrt((1-R2reduced)/(1-R2X.Xk))
	Ar <- (1/sqrt(rchisq(n= M,df=(df+1),ncp=0)))*(tobs*sqrt(rchisq(n= M,df=(df),ncp=0))/sqrt(df) + rnorm(n=M,mean=0,sd=1)) 
	CID <- Ar*D/sqrt(1+Ar*Ar)
	
	initial <- quantile(CID,probs=c(alphaL, alphaU))
	initialqL <- initial[1]
	initialqU <- initial[2]
	lastqL <- initialqL
	lastqU <- initialqU

	lastdenL <- density(CID, from= lastqL, to = lastqL +1)$y[1]
	initialdenL <- lastdenL
	lastdenU <- density(CID, from= lastqU, to = lastqU +1)$y[1]
	initialdenU <- lastdenU
	
	for (n in 1: iter){	
		Ar <- (1/sqrt(rchisq(n= M,df=(df+1),ncp=0)))*(tobs*sqrt(rchisq(n= M,df=(df),ncp=0))/sqrt(df) + rnorm(n=M,mean=0,sd=1)) 
		CID <- Ar*D/sqrt(1+Ar*Ar)

		fnL <- (1-1/n)*lastdenL + (1/n)*(sum(abs(CID-lastqL)<=(1/sqrt(n))))/(2*M*(1/sqrt(n)))
		SnL <- lastqL + (1/(n*max(fnL,initialdenL/sqrt(n))))*(alphaL - sum(CID <= lastqL)/M)
		
		fnU <- (1-1/n)*lastdenU + (1/n)*(sum(abs(CID-lastqU)<=(1/sqrt(n))))/(2*M*(1/sqrt(n)))
		SnU <- lastqU + (1/(n*max(fnU,initialdenU/sqrt(n))))*(alphaU - sum(CID <= lastqU)/M)
		lastqL <- SnL
		lastdenL <- fnL
		lastqU <- SnU
		lastdenU <- fnU
	}
	c(SnL,SnU)
}

################################################
#CI.b:  Function for estimating the confidence interval for an observed unstandardized regression coefficient
#Predictor may either be randomly sampled or fixed -- the result is the same
#Variable definitions:
#  b:  	observed standardized regression coefficient
#  tobs		t-test for the regression coefficient
#  df		df for the t-test
#  conf:	desired confidence level (two-sided)
#  iter:	how many iterations of the stochastic approximation to run.  Over 1000 may be quite slow.
#Returns:
#  1.  Estimates of the upper and lower confidence interval endpoints.
################################################

CI.b<- function(b,tobs,df,conf, iter){
	M <-200000
	alphaL <- (1-conf)/2
	alphaU <- 1- alphaL
	sb <- b/tobs
	CID <- rnorm(n=M,mean=0,sd=1)*sqrt(df)*sb/sqrt(rchisq(n= M,df=(df),ncp=0)) + b
	
	initial <- quantile(CID,probs=c(alphaL, alphaU))
	initialqL <- initial[1]
	initialqU <- initial[2]
	lastqL <- initialqL
	lastqU <- initialqU

	lastdenL <- density(CID, from= lastqL, to = lastqL +1)$y[1]
	initialdenL <- lastdenL
	lastdenU <- density(CID, from= lastqU, to = lastqU +1)$y[1]
	initialdenU <- lastdenU
	
	for (n in 1: iter){	
		CID <- rnorm(n=M,mean=0,sd=1)*sqrt(df)*sb/sqrt(rchisq(n= M,df=(df),ncp=0)) + b

		fnL <- (1-1/n)*lastdenL + (1/n)*(sum(abs(CID-lastqL)<=(1/sqrt(n))))/(2*M*(1/sqrt(n)))
		SnL <- lastqL + (1/(n*max(fnL,initialdenL/sqrt(n))))*(alphaL - sum(CID <= lastqL)/M)
		
		fnU <- (1-1/n)*lastdenU + (1/n)*(sum(abs(CID-lastqU)<=(1/sqrt(n))))/(2*M*(1/sqrt(n)))
		SnU <- lastqU + (1/(n*max(fnU,initialdenU/sqrt(n))))*(alphaU - sum(CID <= lastqU)/M)
		lastqL <- SnL
		lastdenL <- fnL
		lastqU <- SnU
		lastdenU <- fnU
	}
	c(SnL,SnU)
}

#Examples from Biesanz (2009)
CIR.r(r=.612, df=14, conf=.95, iter=200)

CIF.d(d=1.80, n1=17, n2=16, conf=.95, iter=200)