Read in data and relevant packages

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(psych)
library(tidyLPA)
## Note that an update to tidyLPA is forthcoming; see vignette('introduction-to-major-changes') for details!
#stats2 <- read.csv("constrained_parameters.csv")
stats <- read.csv("stats_rev.csv")
colnames(stats)[1] <- "ID"
options(scipen=999)
options(digits=2)
library(stargazer)
## 
## Please cite as:
##  Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.2. https://CRAN.R-project.org/package=stargazer
hU <- read.csv("hUdata.csv")
hUhddm <- inner_join(hU[c(1,14)],stats, by="ID")
hUfull <- read.csv("RNT_wide_final_V7.csv")
#which(colnames(hUfull)=="PHQ9_tot_1")
#which(colnames(hUfull)=="EmotionalImpulsivity_1")
hUhddm <- inner_join(hUhddm,hUfull[c(1,31,49)], by="ID")
hUhddm[c(2:16)] <- scale(hUhddm[c(2:16)])

Regression models using HDDM parameters

RNT regressed on a, t, and the four v parameters from HDDM

lmod1 <- lm(RNT_1 ~ a_mean + t_mean + v_IncSwitch_mean + v_ConSwitch_mean+ v_ConStay_mean + v_IncStay_mean, data=hUhddm)
stargazer(lmod1, type="html")
Dependent variable:
RNT_1
a_mean -0.110
(0.120)
t_mean -0.080
(0.120)
v_IncSwitch_mean 0.360
(0.220)
v_ConSwitch_mean -0.640***
(0.220)
v_ConStay_mean 0.130
(0.230)
v_IncStay_mean 0.077
(0.210)
Constant 0.000
(0.100)
Observations 90
R2 0.140
Adjusted R2 0.082
Residual Std. Error 0.960 (df = 83)
F Statistic 2.300** (df = 6; 83)
Note: p<0.1; p<0.05; p<0.01

RNT is linked with lower drift rate for congruent shift trials, which indicates poorer set shifting.

Cluster models (Latent Profile Analysis) using HDDM parameters

For all cluster models, the optimal number of clusters was determined by comparing BIC values. For parsimony, only the best fitting models are shown here.

Cluster model using means

lpaMod <- estimate_profiles(hUhddm[c(3,5,7,9,11,13)], 
                  a_mean,t_mean, v_ConStay_mean, v_ConSwitch_mean, v_IncStay_mean, v_IncSwitch_mean, 
                  n_profiles=4, return_orig_df = T)
## Fit Equal variances and covariances fixed to 0 (model 1) model with 4 profiles.
## LogLik is 602.575
## BIC is 1353.643
## Entropy is 0.95
plot_profiles(lpaMod, to_center=T, plot_error_bars = F)

Cluster model using SDs

lpaMod2 <- estimate_profiles(hUhddm[c(4,6,8,10,12,14)], 
                            a_std, t_std, v_ConStay_std, v_ConSwitch_std, v_IncStay_std, v_IncSwitch_std, 
                            n_profiles=3, return_orig_df = T)
## Fit Equal variances and covariances fixed to 0 (model 1) model with 3 profiles.
## LogLik is 652.713
## BIC is 1422.421
## Entropy is 0.885
plot_profiles(lpaMod2, to_center=T, plot_error_bars = F)

Cluster model using both means and SDs

lpaMod3 <- estimate_profiles(hUhddm[c(3:14)], 
                             a_mean,a_std,t_mean, t_std, v_ConStay_mean,v_ConStay_std, v_ConSwitch_mean,v_ConSwitch_std, v_IncStay_mean,v_IncStay_std, v_IncSwitch_mean,v_IncSwitch_std, 
                             n_profiles=4, return_orig_df = T)
## Fit Equal variances and covariances fixed to 0 (model 1) model with 4 profiles.
## LogLik is 1223.929
## BIC is 2731.346
## Entropy is 0.973
plot_profiles(lpaMod3, to_center=T, plot_error_bars = F)

huLPA <- cbind(hUhddm, lpaMod[c(7)], lpaMod2[c(7)], lpaMod3[c(13)])

colnames(huLPA)[c(18,19)] <- c("profileSD","profileComb")

Test whether these clusters (profiles) predict differences in RNT using regression

RNT regressed on means clusters

mod <- lm(RNT_1 ~ profile, data=huLPA)
stargazer(mod, type="html")
Dependent variable:
RNT_1
profile2 -0.210
(0.280)
profile3 -0.670**
(0.280)
profile4 -0.730
(0.730)
Constant 0.340
(0.220)
Observations 90
R2 0.076
Adjusted R2 0.043
Residual Std. Error 0.980 (df = 86)
F Statistic 2.300* (df = 3; 86)
Note: p<0.1; p<0.05; p<0.01

Profile 3 is associated with lower RNT.

RNT regressed on SD clusters

mod2 <- lm(RNT_1 ~ profileSD, data=huLPA)
stargazer(mod2, type="html")
Dependent variable:
RNT_1
profileSD2 0.064
(0.390)
profileSD3 -0.100
(0.220)
Constant 0.042
(0.160)
Observations 90
R2 0.003
Adjusted R2 -0.020
Residual Std. Error 1.000 (df = 87)
F Statistic 0.150 (df = 2; 87)
Note: p<0.1; p<0.05; p<0.01

RNT regressed on means and SD combined clusters

mod3 <- lm(RNT_1 ~ profileComb, data=huLPA)
stargazer(mod3, type="html")
Dependent variable:
RNT_1
profileComb2 -0.250
(0.300)
profileComb3 -0.620**
(0.280)
profileComb4 -0.510
(0.440)
Constant 0.380
(0.230)
Observations 90
R2 0.060
Adjusted R2 0.027
Residual Std. Error 0.990 (df = 86)
F Statistic 1.800 (df = 3; 86)
Note: p<0.1; p<0.05; p<0.01

Profile 3 is associated with lower RNT.

Read in behavioral data

hUBehavioral <- hUfull[c(1,6:10,33)]

hUBehavioral[c(2:7)] <- scale(hUBehavioral[c(2:7)])
hUBehavioral <- hUBehavioral[complete.cases(hUBehavioral),]

Note that RNT is called “PTQ” in this dataset but it is the same variable.

Regression model with behavioral data

modlmB <- lm(PTQ_tot_1 ~ sw_costRT + sw_avg_ac + sw_inc_ac + sw_acc + sw_iacc, data=hUBehavioral)
stargazer(modlmB, type="html")
Dependent variable:
PTQ_tot_1
sw_costRT 0.037
(0.130)
sw_avg_ac 0.130
(0.790)
sw_inc_ac 0.280
(0.600)
sw_acc -0.710
(0.750)
sw_iacc 0.330
(0.540)
Constant -0.052
(0.100)
Observations 86
R2 0.071
Adjusted R2 0.013
Residual Std. Error 0.970 (df = 80)
F Statistic 1.200 (df = 5; 80)
Note: p<0.1; p<0.05; p<0.01

There are no significant associations between the behavioral data and RNT.

Clustering and regression with behavioral data

lpaBehavioral <- estimate_profiles(hUBehavioral[c(2:6)], 
                             sw_costRT, sw_avg_ac, sw_inc_ac, sw_acc, sw_iacc, 
                             n_profiles=4, return_orig_df = T)
## Fit Equal variances and covariances fixed to 0 (model 1) model with 4 profiles.
## LogLik is 342.351
## BIC is 809.423
## Entropy is 0.963
plot_profiles(lpaBehavioral, to_center=T, plot_error_bars = F)
## Warning: attributes are not identical across measure variables;
## they will be dropped

hULPABehavioral <- cbind(hUBehavioral, lpaBehavioral[c(6)])


clustModB <- lm(PTQ_tot_1 ~ profile, data=hULPABehavioral)
stargazer(clustModB, type="html")
Dependent variable:
PTQ_tot_1
profile2 0.100
(0.340)
profile3 0.230
(0.340)
profile4 0.190
(0.380)
Constant -0.200
(0.280)
Observations 86
R2 0.007
Adjusted R2 -0.030
Residual Std. Error 0.990 (df = 82)
F Statistic 0.180 (df = 3; 82)
Note: p<0.1; p<0.05; p<0.01

There are no significant associations.