analysis.Rmd

---
title: "Examining Associations Between Psychopathic Traits and Executive Functions in Incarcerated Violent Offenders"
author: "Carl Delfin"
date: "30 May 2018"
output:
  pdf_document: default
  html_document: default
  word_document: default
---

# Introduction

This RMarkdown document contains the R code used for all analysis in the study [Examining Associations Between Psychopathic Traits and Executive Functions in Incarcerated Violent Offenders](https://www.frontiersin.org/articles/10.3389/fpsyt.2018.00310/full).

# Clear environment

Clear any existing items from the global environment.

```{r}
rm(list = ls())
```

# Initial set-up

## Start timer

Not really necessary *per se*, but the analysis does take some time to run, owing to the MCMC procedure.

```{r}
# load tictoc package
library("tictoc")

# clear timer and timer log
tic.clear()
tic.clearlog()

# start timer
tic("Full analysis")
```

## Paths and keys

Paths are stored in an external script for security reasons. 

```{r}
source("hidden/path.R")
```

## Required packages 

```{r}
requiredPackages <- c("foreign",
                      "psych",
                      "ReporteRs",
                      "BayesMed",
                      "reshape2",
                      "data.table",
                      "cowplot",
                      "ggplot2")
```

The *ipak.R* script goes through the `requiredPackages` vector and checks if the required packages are installed, installs them if not, then loads them. Note that other packages are loaded as dependencies. For a complete overivew, see *sessioninfo.txt*.

```{r message=FALSE, warning=FALSE}
source("scripts/ipak.R")
ipak(requiredPackages)
rm(requiredPackages, ipak)
```

## Data and seed

The data is read from an SPSS file and the seed is set for reproducibility (I always use the year the analysis was conducted, in this case 2018, to avoid "seed hacking"). The `seed` function uses [Mersenne Twister](https://en.wikipedia.org/wiki/Mersenne_Twister) by default.

```{r}
data <- read.spss(DATALOCATION, to.data.frame = TRUE)
rm(DATALOCATION)
seed <- 2018
```

# Data preparation

## Select variables

Selec which variables to include in the analysis.

```{r}
# outcome variables
psychopathicTraits <- c("facet1", "facet2", "facet3", "facet4")

# EF predictors
executiveFunctions <- c("IEDstages", "IEDtotaler",
                        "SWMtotaler", "SWMstrategy",
                        "SSTSSRT", "SSTMRT",
                        "SOCMITT5", "SOCPS")

# PCL-R items
pclItems <- c("pcl1", "pcl2", "pcl3", "pcl4", "pcl5",
              "pcl6", "pcl7", "pcl8", "pcl9", "pcl10",
              "pcl11", "pcl12", "pcl13", "pcl14", "pcl15",
              "pcl16", "pcl17", "pcl18", "pcl19", "pcl20")

# other variables of interest
otherVariables <- c("totalscore", "age")

# subset data to selected variables
data <- data[, c(pclItems, otherVariables, psychopathicTraits, executiveFunctions)]

# drop missing data
data <- data[complete.cases(data[, c(executiveFunctions, psychopathicTraits)]), ]

# convert ms to s
data$SOCMITT5 <- data$SOCMITT5 / 1000
data$SSTSSRT <- data$SSTSSRT / 1000
data$SSTMRT <- data$SSTMRT / 1000
```

# Descriptive information

## Table 1

Table 1 in the manuscript consists of an overview of numerical variables. See *table1.R* for details.

```{r}
source("scripts/table1.R")
```

## Internal consistency

[Internal consistency](https://en.wikipedia.org/wiki/Internal_consistency) of the PCL-R scores is assessed using `psych::alpha`.

```{r}
source("scripts/internalconsistency.R")
```

# Zero-order correlations

## Classical analysis

Classical in this case is just plain old [null hypothesis statistical testing](https://en.wikipedia.org/wiki/Statistical_hypothesis_testing).

```{r}
source("scripts/classicalanalysis.R")
```

## Bayesian analysis

**Note!** This step uses [Markov chain Monte Carlo](https://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo) to estimate the posterior distribution, which takes a significant amount of time.

```{r message=FALSE, warning=FALSE, include=FALSE}
# select number of MCMC samples and number of burn-in samples
niter <- 20000
nburnin <- 2000

# start timer and grab some popcorn
tic("Bayesian analysis using MCMC")
source("scripts/bayesiananalysis.R")
toc()
```

# Figure 1

Figure 1 is a four-panel overview Pearson's *r* values, *p*-values, posterior probability and Bayes factors. See *figure1.R* for details.

```{r}
source("scripts/figure1.R")
ggsave(filename = "hidden/figures/figure1.eps",
       plot = bayesplot,
       height = 7.5,
       width = 13,
       dpi = 900,
       units = "in")
rm(bayesplot)
```

# Session info

Save a text file with full session info for future reference.

```{r}
sessInfo <- capture.output(sessionInfo())
write(sessInfo, file = "log/sessioninfo.txt")
```

# Stop timer

Analysis is done! Let's stop the timer.

```{r}
toc()
```