#
#
############ Moderated regression and simple slopes analysis (with graphs) ############
#
#
# This script assumes you have downloaded and installed all available packages.
# Packages needed for this script will be noted as "library(xxxxx)" where xxxxx
# refers to the required package/library to be loaded.
# Make sure you have the R-Commander library loaded.
library(foreign)
library(Rcmdr)
# Use the 'foreign' library to load the "Cars_Trimmed.sav" data file.
CarsTrimmed <- read.spss("http://researchsupport.unt.edu/class/Jon/R_SC/Module9/Car_Trimmed.sav",
use.value.labels=TRUE, max.value.labels=Inf, to.data.frame=TRUE)
# From here on, simply use this script; submitted or pasted into the
# R Console.
attach(CarsTrimmed)
names(CarsTrimmed)
nrow(CarsTrimmed)
cor(CarsTrimmed)
##########################################################################
# Analyze the moderating influence of accel on the relationship
# between mpg and weight.
# First, center the predictor (weight) and the moderator (accel); which is
# simply subtracting the mean of each variable from each variable.
mean(weight)
mean(accel)
c.weight <- weight - 2960.374
c.accel <- accel - 15.54296
cars1 <- cbind(CarsTrimmed, c.weight, c.accel)
rm(c.accel, c.weight)
names(cars1)
detach(CarsTrimmed)
attach(cars1)
# Create the interaction product term (c.weight * c.accel).
interact <- c.weight * c.accel
cars2 <- cbind(cars1, interact)
detach(cars1)
rm(cars1,interact)
attach(cars2)
head(cars2)
# Run two regressions and compare to determine if the interaction is meaningfully contributing
# to the model.
RegMod.1 <- lm(mpg~c.weight+c.accel)
summary(RegMod.1)
RegMod.2 <- lm(mpg~c.weight+c.accel+interact)
summary(RegMod.2)
anova(RegMod.1,RegMod.2)
# The Akaike Information Criterion (AIC).
AIC(RegMod.1)
AIC(RegMod.2)
# The Bayesian Information Criterion (BIC)
AIC(RegMod.1, k = log(nrow(CarsTrimmed)))
AIC(RegMod.2, k = log(nrow(CarsTrimmed)))
##########################################################################
########## Simple Slopes Analysis (t-test) ##########
# Construct the resulting linear regression equation using the "Coefficients"
# table from
RegMod.2$coefficients
# mpg = 23.09228 + -.00738(c.weight) + .38335(c.accel) + -.00045(interact).
# Or, we could write the same equation like this:
# mpg = 23.09228 + -.00738(c.weight) + .38335(c.accel) + -.00045(c.weight * c.accel).
# First, calculate the simple slope for the 'low group' where low corresponds to
# 1 standard deviation below the mean of c.accel (-1 SD = -2.77640).
# You may be confused because we centered the predictor and moderator; but centering
# only changes the location of the mean, it does not change the standard deviation.
sd(c.accel)
sd(accel)
# So, we use the regression equation from above and substitute -2.77640 (one standard
# deviation below the mean) in place of 'c.accel', which results in the following equation:
# y = 23.091 + (-.00738*c.weight) + (.38335*-2.77640) + (-.00045*c.weight*-2.77640)
# which can be reduced down to give us the simple slope we want:
# y = 23.091 + (.38335*-2.77640) + (-.00738*c.weight) + (-.00045*c.weight*-2.77640)
# y = 23.091 + -1.064333 + [-.00738 + (-.00045*-2.77640)]*c.weight
# y = 22.02667 + [-.00738 + .00124938]*c.weight
# y = 22.02667 + -.00613062*c.weight
# So, the simple slope for the 'low group' is -.00613062.
# Second, calculate the simple slope for the 'high group' where high corresponds
# to 1 standard deviation above the mean of c.accel (+1 SD = 2.77640).
# y = 23.091 + (-.00738*c.weight) + (.38335*2.77640) + (-.00045*c.weight*2.77640)
# which can be reduced down to give us the simple slope we wanted:
# y = 24.15533 + -0.00862938*c.weight
# So, the simple slope for the 'high group' is -.00862938.
# Next, we have to calculate the standard error for each simple slope/group.
# To do this; we need to the coefficient covariances.
summary(RegMod.2)$cov.unscaled
# We will need 3 things from this table.
# (1) the variance of the coefficient for c.weight = 4.630701e-09
# (2) the variance of the coefficient for interact = 4.732109e-10
# (3) the COVARIANCE of the coefficients for c.weight and interact = 4.708424e-10
# Now we can calculate the Standard Error of each group specific simple slope using
# the formula:
# std.err. = sqrt((1) + S*S*(2) + 2*S*(3))
# which can also be written as:
# std.err. = sqrt of (4.630701e-09 + S*S*4.732109e-10 + 2*S*4.708424e-10)
# where S corresponds to the group specific value for standard deviation;
# for the low group: S = -2.77640 and
# for the high group: S = 2.77640.
# For the low group: S = -2.77640.
std.err.low = sqrt(4.630701e-09 + -2.77640*-2.77640*4.732109e-10 + 2*-2.77640*4.708424e-10)
std.err.low
# Which equates to just 0.000075259
# For the high group: S = 2.77640.
std.err.high = sqrt(4.630701e-09 + 2.77640*2.77640*4.732109e-10 + 2*2.77640*4.708424e-10)
std.err.high
# Which equates to just 0.000104369
# Finally, we can conduct the t-tests to determine if each group specific simple slope is
# significantly different from zero.
# We calculate our t-value by simply dividing the slope by the standard error (for each group).
# BUT, remember we are doing multiple t-tests here, so we risk inflation of type 1 error rates.
# To overcome this, we must apply a simple Bonferroni correction to our alpha
# level (.05 / 2 comparisons = .025).
# For the low group:
low.slope.t <- -.00613062 / 0.000075259
low.slope.t
# For the high group:
high.slope.t <- -.00862938 / 0.000104369
high.slope.t
# To find our critical t-value, we use degrees of freedom (df) = N - k - 1 where N is the
# number of cases, and k is the number of predictors (including the predictor, moderator,
# & interaction term).
# df = 398 - 3 - 1 = 394
# This df yeilds a two-tailed critical t-value of 3.347468181818182 with alpha = .001.
# We can see that both simple slopes are significantly different from zero.
##########################################################################
# Splitting 'accel' into two equal groups (high & low based on ranking
# or sorting of accel).
a.group <- c(rep(1,199),rep(2,199))
head(cars2)
tail(cars2)
cars3 <- data.frame(cars2)
cars3$rank <- rank(cars3$accel, ties.method = "first")
rownames(cars3) <- cars3$rank
cars4 <- cars3[as.character(sort(cars3$rank)),]
cars5 <- cbind(cars4,a.group)
rm(a.group,cars2,cars3,cars4)
head(cars5)
tail(cars5)
############################################################################
# Graphing the interaction.
# Read what it below the two commands prior to running them.
library(latticist)
latticist(cars5,use.playwith=FALSE)
# SELECT 2: gWidgetscltk
Selection: 2
# Expect errors; but the errors do not prevent functioning.
# When the 'Latticist' GUI pops up, put 'weight' in the 'x=' box,
# 'mpg' in the 'y=' box, 'a.group' in the 'Groups/Color' box AND
# put 'a.group' in the 'Conditioning' box. This should produce three
# graphs, with each group of 'a.group' stacked and the total sample
# (both groups together) to the right. Each graph is also color
# coded so the points of our groups are differentiated.
# If you only want a single graph with both groups color coded
# similar to the third one produced above, then in the Latticist
# GUI:
# Click on the button for "superpose" which will remove
# 'a.group' from the 'Conditioning' box and produce the single
# desired graph with different colors for each group's points and
# their lines.
# There are several ways to display data in the Latticist GUI, so
# spend some time changing things around a bit.
# For instance, try putting 'weight' in the 'x=' box,
# 'mpg' in the 'y=' box, 'accel' in the 'Groups/Color' box AND
# put 'accel' in the 'Conditioning' box.
# As an example of the kinds of graphs one can produce using Latticist,
# copy and paste the following URL into your browser and hit the enter key:
# http://researchsupport.unt.edu/class/Jon/R_SC/Module9/LatticistGraph001.pdf