"prodclin" <-
function(a,sea,b,seb,rho,alpha)
{
#Prodclin.r - MacKinnon, Fritz, Williams, and Lockwood, 2005       
#Dept. of Psychology, Arizona State University
#Necessary files: prodclin.r, prodclin.exe                       

#The prodclin.r program calculates asymmetric confidence intervals for 
#the mediated effect using the distribution of the product of two    
#random variables and the Fortran program Prodclin.exe.  To run, type
#'source("F:/prodclin.r")'  in the R GUI, where F:/ is replaced with the
#location of the file prodclin.r  Then, the location of the Fortran program
#must be entered in four places (both files should be in the same location).
#Next, type 'prodclin(a,sea,b,seb,rho,alpha)' in the R GUI, a is the value
#of the a path, sea is the st. error of a, b is the value of the b path, 
#seb is the st. error of b, rho is the correlation between a and b
#(usually set to zero), and alpha is the type 1 error rate.  The output
#consists of the user defined values of a, sea, b, seb, rho, alpha, ab,
#the first order standard error (sobel), the lower and upper standardized 
#critical values from Prodclin [low,high], the asymmetric confidence limits 
#[prodlow,produp] and the corresponding normal theory confidence limits 
#[normlow,normup].

#All places that need user input have been labeled with comments.
     
#Change file address to match the location of the file prodclin.exe
write(c(a,sea,b,seb,rho,alpha),file="G:/raw.txt",ncolumns=6)

#Change file address to match the location of the file prodclin.exe
setwd("G:/")

#Change file address to match the location of the file prodclin.exe
shell("G:/prodclin.exe")

#Change file address to match the location of the file prodclin.exe
critval <- read.table("G:/critval.txt",header=FALSE)
lcrit <- critval[1,1]
ucrit <- critval[1,2]
ab <- a*b
sobel <- sqrt(a*a*seb*seb+b*b*sea*sea)
da<-a/sea
db<-b/seb
dadb<-da*db
sedadb<-sqrt(da*da+db*db+1)
low <- (lcrit-dadb)/sedadb
high <- (ucrit-dadb)/sedadb

prodlow <- ab+low*sobel
produp <- ab+high*sobel

nl<-qnorm(alpha/2)
normlow <- ab+nl*sobel
normup <- ab-nl*sobel

if(lcrit==99999 || ucrit==99999) {
   out<- "Did not converge.  Please check input values or try again with fewer decimal places" }
else {
   out<-list(a=a, sea=sea, b=b, seb=seb, ab=ab, "Sobel se"=sobel, rho=rho,"type1 error"=alpha, 
	"Normal lower limit"=normlow,"Normal upper limit"=normup,"Prodclin lower critical value"=low,
	"Prodclin upper critical value"=high,"Prodclin lower limit"=prodlow, "Prodclin upper limit"=produp)}
out
}