21-2way_simple_slopes.R

#########################################################################
##################### 2-way Simple Slopes Testing #######################
#########################################################################

#####
# Date: TODAY'S DATE
# By: *INSERT NAME HERE*
# Description: PROJECT DISCRIPTION
# Version of R used: CURRENT VERSION OF R


#####################################
### Import data & load libraries ####
#####################################

## Install the required script packages if not yet installed

# Install pacman, jmv, & reghelper package if necessary
if(!"pacman" %in% rownames(installed.packages())) install.packages("pacman")
if(!"jmv" %in% rownames(installed.packages())) install.packages("jmv")
if(!"reghelper" %in% rownames(installed.packages())) install.packages("reghelper")
if(!"see" %in% rownames(installed.packages())) install.packages("see")

pacman::p_load(rio, pequod, jtools, interactions, 
               apaTables, stargazer, psych, tidyverse, 
               parameters, performance, effectsize,
               janitor, ggstance, patchwork) # this row is needed for plotting

## load data

# RData files work the best in R. 
# Try to only open RData files to avoid any issues.
# CSV works the next best in R. 
# Try to only save and open CSV files to avoid any issue if you cannot load RData files.
# SPSS files can be buggy to import, especially factors and labels

# the following command will open a dialog box and allow you to select the file you wish to laod
dat <- import(file.choose()) 


# check to see that you loaded the correct dataset
View(dat)

# list variables in dataset
glimpse(dat)

# NB: Alt + Shift + K --> will bring up keyboard shortcuts

################################
###### Summary Statistics ######
################################

## descriptive summary table 
# note: na.omit() removes any NAs contained within each of the IVs

dat %>%
      dplyr::select(iv1, iv2, dv) %>%
      describe()

####### center and scale IVs
dat$c_iv1 <- effectsize::standardize(dat$iv1, two_sd = FALSE, force = TRUE)
dat$c_iv2 <- effectsize::standardize(dat$iv2, two_sd = FALSE, force = TRUE)

# verify centering
dat %>%
  dplyr::select(starts_with("c_")) %>%
  describe()


##############################
###### 2-way Regression ######
##############################

####### test 2-way regression interaction

## regression function 
# simplified code to run regressions on all DVs at once

reg_models <- dat %>% 
                     dplyr::select(starts_with("avg_")) %>%   # this line tells the map() only use your DVs (all start "avg_" in my datasets)
                      map(~summ(lm(. ~ c_iv1 * c_iv2, data = dat))) 

reg_models 

### linear regression to dive into any significant 2-ways
step1.1 <- lm(dv ~ c_iv1 + c_iv2, data=dat)
step2.1 <- lm(dv ~ c_iv1 * c_iv2, data=dat)

##  check GLM assumptions for:
#heteroskedastic (error variance)
check_heteroscedasticity(step2.1)

# autocorrelation (independence of errors)
check_autocorrelation(step2.1)

# normality (normality of residuals)
check_normality(step2.1)

# multicollinearity (predictor independence)
check_collinearity(step2.1)

# tests for outliers in model then iteratively removes outliers and re-runs the model
# consider add "mcd" method to detect outliers (Mahalanobis et al., 2018)
check_outliers(step2.1, method = c("cook", "mahalanobis"))

# can also plot GLM assumptions
performance::check_model(step1.1)
performance::check_model(step2.1)


## SPSS-like regression summary
jmv::linReg(data = dat,
            dep = dv,
            covs = vars(c_iv1, c_iv2),
            blocks = list(
              list("c_iv1", "c_iv2"),
              list(c("c_iv1", "c_iv2"))),
            r2 = FALSE, r2Adj = TRUE, ci = TRUE, stdEst = TRUE,
            ciEmm = FALSE, emmPlots = FALSE, emmWeights = FALSE)


##############################################################
######## could also achive this differently by doing: ########
##############################################################

# hierarchical linear regression
reghelper::build_model(dv, c(c_iv1 + c_iv2), 
                       c(c_iv1 * c_iv2), 
                       data=dat, model='lm') %>% summary()

# step 1 betas
model_parameters(step1.1, standardize = "basic")

# step 2 betas
model_parameters(step2.1, standardize = "basic")


######################################################################
############### Simple Slope Testing Automatically ###################
######################################################################

### test simple slopes 

# Johnson-Neyman intervals with plots
sim_slopes(step2.1, pred = c_iv1, modx = c_iv2, jnplot = TRUE)
sim_slopes(step2.1, pred = c_iv2, modx = c_iv1, jnplot = TRUE)

# simple slopes plots with Johnson-Neyman intervals in output
probe_interaction(step2.1, 
                  pred = c_iv1, modx = c_iv2,
                  modx.values = "plus-minus")


probe_interaction(step2.1, 
                  pred = c_iv2, modx = c_iv1,
                  modx.values = "plus-minus")


###############################################
######### Plotting 2-way interaction ##########
###############################################

# the follow commands will create APA-style, MS ready figures

### create x, z, and y columns by renaming IVs
dat$z <- dat$c_iv1        # x-axis variable here
dat$x <- dat$c_iv2        # moderator_1 variable here
dat$y <- dat$dv           # outcome variable here


## create labels for figure
# must use quotes for labels
# change labels in quotes to be what you want them to be

y_label        <- "dv_name"
y_axis_high    <- 9.0                    # high descrete numeric value displayed on y-axis
y_axis_low     <- 1.0                    # low descrete numeric value displayed on y-axis
y_increment    <- 1.0                    # increments for y-axis
z_label        <- "z_non-moderator_name" # X-axis variable (non-moderator)
z_values       <- c("low", "high")       # non-moderator values
z_range        <- c(-1.0, 1.0)           # non-moderator numerical values
modx_label     <- "x_moderator_name"     # figure legend (moderator)
modx_values    <- c("low", "high")       # moderator values

# run but only modify these if needed
# line type and color are set for APA 7ed
line_types    <- c("longdash", "solid")   # can be “solid”, “dashed”, “dotted”, “dotdash”, “longdash”, “twodash”
line_colors   <- c("#9E9E9E", "#000000")  # these are colorblind friendly color schemes, don't change unless needed

# change legend location if needed
# in the format of (x,y), can be any number between 0 and 1
# (0.8, 0.8) is upper-right corner, (0.2, 0.8) is the top left of the figure
legend_loc    <- c(0.2, 0.8)

# run next line and figure will be automatically created
source("https://raw.githubusercontent.com/rastlab/R_templates/master/99991-ggplot2_2way_SS.R")


#####################################
####### Here are the figures ########
#####################################

# verify figures look similar to the probe_interaction() commands at line 146 or 151
figure_1

# if you are happy with figures, save them
# can change dimensions, file type, and dpi as per journal requirement specifications
ggsave('./figures/figure_1.png', panel_1, width = 8, height = 8, unit = 'in', dpi = 320)


#######################################
###### Create descriptives table ######
#######################################

# the below script is set up so that once YOU update the variable names in your analysis, 
# then the below commands will automatically create a descriptives table that can be placed
# directly into your MS, presentation, or poster
# we'll also remove NA values to make this simpler

dat3 = na.omit(dat %>% 
                 dplyr::select(iv1, iv2, dv) %>% 
                 rename(NEW_NAME_IV1 = iv1, # relabel whatever you want your variables to be named in the manuscript, cannot contain spaces though
                        NEW_NAME_IV2 = iv2, 
                        NEW_NAME_DV  = dv))

# correlation table
apa.cor.table(dat3, filename = "./tables/correlation_table.doc", table.number = 1,
              show.conf.interval = FALSE, landscape = TRUE)


#######################################
###### Saving Data and Workspace ######
#######################################

# Save current workspace: 
export(dat, "./data/21_two_way_regression.RData")

# save R data file as CSV
export(dat, "./data/21_two_way_regression.csv")

# save R data file as SAV SPSS file
export(dat, "./data/21_two_way_regression.sav")