knitr::opts_chunk$set(echo = TRUE, warning=FALSE)

R Markdown

This is the code to reproduce my moderated mediation analysis on the relationship between violence exposure and externalizing behavior symptoms among male and female youth who have involvement with the child welfare system and who have manifested symptoms of Post-traumatic stress. http://rmarkdown.rstudio.com.

To read similar papers, click here.

library(tidyverse)
library(knitr)
library(lavaan)
library(psych)
library(MBESS)
library(xtable)
library(semPlot)
require(car)
library(dplyr)
library(lattice)
library(lavaanPlot)

First open the text file and chooise the variables to include. In this part of the analysis, I chose to explore the unadjusted relationships between symptoms of externalizing behavior and violence exposure as expressed through post-traumatic stress.

pts_vex <- read.csv("/media/barbozag/easystore/Research/
          Mediation-Moderation Models/Data/test7.csv")


analysis_vars <- c("tra1", "bcext1",  "GENDER",  "EV_W1")
pts_vex <- pts_vex[analysis_vars]
pts_vex <- pts_vex[complete.cases(pts_vex),]

names(pts_vex)[1]<-"M"
names(pts_vex)[2]<-"Y"
names(pts_vex)[3]<-"W"
names(pts_vex)[4]<-"X"

pts_vex  %>%
  headTail() %>%
  kable()

pts_vex %>%
  dplyr::select(X, M, Y) %>%
  pairs.panels(scale = FALSE, pch = ".", ci = TRUE, alpha=.1)

Including Plots

You can also embed plots, for example:

pts_vex$W_f  <- as.factor(pts_vex$W  )
pts_vex$W_f  <- factor(pts_vex$W , levels=c(0,1), labels=c("Male", "Female"))

scatterplot(Y~M|W_f  , data=pts_vex, pch=c(3,21), regLine=F, smooth=F,
            boxplots="xy", main="Scatterplot of Externalizing v PTS by Gender")

panel.smoother <- function(x, y) {
  panel.grid(h=-1, v=-1)
  panel.xyplot(x, y)
  panel.loess(x, y)
  panel.abline(h=65, lwd=2, lty=2, col="green")
}
xyplot(Y~M|W_f , data=pts_vex,
       scales=list(cex=.8, col="red"),
       panel=panel.smoother,
       xlab="PTS", ylab="EXT",
       main="EXT vs PTS by Gender",
       sub = "Dotted Lines Represent Clinical Cutoff", aspect=1)

pts_vex$X_cut<-cut(pts_vex$X, c(0,1,2,7), include.lowest = TRUE, right = TRUE)

xyplot(Y~M|W_f*factor(X_cut), data=pts_vex,
       scales=list(cex=.8, col="red"),
       panel=panel.smoother,
       xlab="PTS", ylab="EXT",
       main="EXT vs PTS by Gender",
       sub = "Dotted Lines Represent Clinical Cutoff", aspect=1)

model4 <- "# a path
           M ~ 1 + a*X

           # b path
           Y ~ 1 + b*M

           # c_prime path
           Y ~ cp*X

           #indirect and total effects
           ind_eff:= a*b
           tot_eff:=cp+ind_eff
"

set.seed(1234)
sem.fit1 <- sem(model4, data = pts_vex, se = "bootstrap", bootstrap = 10000)

summary(sem.fit1, standardized = TRUE)

parameterestimates(sem.fit1, boot.ci.type = "bca.simple", standardized = TRUE) %>%
  kable()

semPaths(sem.fit1, style="lisrel",
         whatLabels = "std", edge.label.cex = .8,
         label.prop=0.8, edge.label.color = "black", rotation = 2,
         equalizeManifests = FALSE,  node.width = 1.5,
         edge.width = 0.5, shapeMan = "rectangle",
         shapeInt = "triangle", sizeMan = 6, sizeInt = 2, sizeLat = 4,
         curve=1, unCol = "#070b8c")

labels <- list(X = "VEX",
               M = "PTS", Y = "EXT")

lavaanPlot(model = sem.fit1, labels=labels, node_options = 
          list(shape = "box", fontname = "Helvetica"),
          edge_options = list(color = "grey"), coefs = T, sig = 1, stars = TRUE)

estVals <- parameterEstimates(sem.fit1)

xtab2 =data.frame(e = estVals[estVals$op == '~', 'lhs'],
                  vals = paste0(format(estVals[estVals$op == '~', 'est'], 
                  digits = 3), ifelse(estVals[estVals$op == '~', 'pvalue'] 
                  < .05, '*', '')))

xtab2 <-
  estVals %>%
  filter(estVals$op=='~' | estVals$op==':=') %>%
  select(lhs, est, ci.lower, ci.upper) %>%
  mutate_if(is.numeric, round, 3)

print(xtab2, include.rownames = FALSE)
xtable(xtab2, caption="Mediation Effect of Post-Traumatic Stress Symptoms")

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.