# To plot a simple interaction: Y = X + Z + XZ

# Function will ask for y, x, z and data and whether x and z are already centered
# Plot will be an simple x,y scatter plot
# with ablines added (Y=X for +/- 1 sd(z))

# Plot limits are set as +/- 2 SD(x,y) respectively



# REF: Aiken & West (1991)
# REF: CCWA (2003)

## Steps (or f(x)):
## 1. get data, mean center (x,z), create XZ and then run lm(): moderate.lm
## 2. from moderate.lm - calculate simple slopes and return ss and their se
## 3. graph simple slopes on an x,y plot see above

#### Additional notes for dealing with categorical x or z
#### must not center, treating as factor with appropriate contrasts will do this
#### need to add if/else statements to test for ?is.factor
#### probably will have to write new f(x)s for graphing

########################################################################
# Functions follow #
# INDEX of FUNCTIONS

# 	mean.center
#	moderate.lm
#	sim.slopes
#	graph.mod
#	

mean.center <- function (x,data)
{
	mcx <- 	x - mean(x, na.rm=TRUE)
	return(mcx)
	}
	
# initial moderation model 

moderate.lm <- function (x,z,y,data,mc=FALSE)
{
	
	if(!mc)
	{	mcx <- 	x - mean(x, na.rm=TRUE)
		mcz <- 	z - mean(z, na.rm=TRUE)}
	else {
		mcx <- x
		mcz <- z 	}
		
	lm1 <- lm(y ~ mcx*mcz, na.action=na.omit)	
		}
	

### simple slopes - based on object (moderate.lm or other lm with 3 predictors (x,z,xz))
### SS are the slopes (and SE) for the slopes of lines y~x at certain values of z (+/- 1 sd)

sim.slopes <- function (mod,z)
{ 	
		
	zhi <- mean(z) + sd(z, na.rm=TRUE)
	zlo <- mean(z) - sd(z, na.rm=TRUE)
	zme <- mean(z)
		b0 <- summary(mod)$coef[1,1]
		b1 <- summary(mod)$coef[2,1]
		b2 <- summary(mod)$coef[3,1]
		b3 <- summary(mod)$coef[4,1]
	
	x.zhi <- (b1 + b3*zhi)    # simple slope
	x.zlo <- (b1 + b3*zlo)    # simple slope
	x.zme <- (b1 + b3*zme)
	int.zhi <- (b0 + b2*zhi)  # simple intercept
	int.zlo <- (b0 + b2*zlo)  # simple intercept
	int.zme <- (b0 + b2*zme)
	

		seb.zhi <- sqrt(vcov(mod)[2,2] + 2*zhi*vcov(mod)[2,4] + zhi^2*vcov(mod)[4,4])
		seb.zlo <- sqrt(vcov(mod)[2,2] + 2*zlo*vcov(mod)[2,4] + zlo^2*vcov(mod)[4,4])
		seb.zme <- sqrt(vcov(mod)[2,2] + 2*zme*vcov(mod)[2,4] + zme^2*vcov(mod)[4,4])
		td <- qt(.975, df = summary(mod)$df[2])
		
	zhi.u <- x.zhi + td*seb.zhi  # upper limit of simple slope
	zhi.l <- x.zhi - td*seb.zhi  # lower limit of simple slope
	zlo.u <- x.zlo + td*seb.zlo  # upper limit of simple slope
	zlo.l <- x.zlo - td*seb.zlo  # lower limit of simple slope
	zme.u <- x.zme + td*seb.zme
	zme.l <- x.zme - td*seb.zme
	
	# create a matrix with int & slope of y~x at 2 levels of z; se.z; 95% CI 
	mat <- matrix(NA,3,5)
	colnames(mat) <- c("INT", "Slope", "SE", "LCL", "UCL")
	rownames(mat) <- c("at zHigh", "at zMean", "at zLow")
	mat[1,] <- c(int.zhi,x.zhi,seb.zhi,zhi.l,zhi.u)
	mat[2,] <- c(int.zme,x.zme,seb.zme,zme.l,zme.u)
	mat[3,] <- c(int.zlo,x.zlo,seb.zlo,zlo.l,zlo.u)
	return(data.frame(mat))
	
	}
	

##### GRAPH using Simple slopes and plot(x,y)

### requires mean.center of X 
### X & Y should come from a data.frame

graph.mod <- function(ssmod, x, y, data, title="",xlab="Centered X",ylab="Raw Y",xll=-2,yuu=2, ...)
{
	
	mcx <- mean.center(x,data)
	plot(mcx,y,main=title, xlab=xlab, ylab=ylab)
	abline(ssmod[1,1],ssmod[1,2],col=4)
	abline(ssmod[3,1],ssmod[3,2],col=2)
	abline(ssmod[2,1],ssmod[2,2],lty=2, col=1)
	legend (xll,yuu, c("zHigh", "zLow"), lty=1, col=c(4,2))
	}

	
#######

#### Examples #####

# create some data depicting an interaction

x <- rnorm(100)
z <- .28*x + .9*rnorm(100)
xz <- x*z
y <- .3*x + .3*z - .4*xz + .2*rnorm(100)

d <- data.frame(x,z,y)

tlm <- moderate.lm(x,z,y,d)
tss <- sim.slopes(tlm,mean.center(d$z))
tss
graph.mod(tss,x,y,d,"Interaction")