## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(psychTools)
library(psych)
library(OpenMx)
library(semtree)

## -----------------------------------------------------------------------------
data(spi)
data(spi.keys)

# score all available scales and extract two broad dimensions
spi_scored <- scoreFast(keys = spi.keys, items = spi)


selected_traits <- c("Agree-A", "Extra-A")

trait_scores <- spi_scored[, selected_traits, drop = FALSE]

# potential predictors available in the dataset
candidate_covariates <- intersect(c("age", "education", "gender"), names(spi))
spi_analysis <- cbind(trait_scores, spi[, candidate_covariates, drop = FALSE])
spi_analysis <- na.omit(spi_analysis)

selected_traits <- c("Agreeableness", "Extraversion")
names(spi_analysis)[1:2] <- selected_traits

# add some noise
set.seed(104358)
spi_analysis$noise1 <- ordered(sample(c(0,1), nrow(spi_analysis), TRUE))
spi_analysis$noise2 <- ordered(sample(c(1:10), nrow(spi_analysis), TRUE))
spi_analysis$noise3 <- ordered(sample(c(1:10), nrow(spi_analysis), TRUE))
spi_analysis$noise4 <- ordered(sample(c(1:20), nrow(spi_analysis), TRUE))

spi_analysis$education <- ordered(spi_analysis$education)

## -----------------------------------------------------------------------------

two_dim_model <- mxModel(
  "TwoDimensions",
  type = "RAM",
  manifestVars = selected_traits,
  mxData(spi_analysis[, selected_traits, drop = FALSE], type = "raw"),
  # freely estimate the variances and their covariance
  mxPath(
    from = selected_traits,
    arrows = 2,
    connect = "unique.pairs",
    free = TRUE,
    values = c(1, 0.2, 1),
    labels = c("var_trait1", "cov_traits", "var_trait2")
  ),
  # estimate manifest means
  mxPath(
    from = "one",
    to = selected_traits,
    arrows = 1,
    free = TRUE,
    values = 0,
    labels = paste0("mean_", selected_traits)
  )
)

two_dim_fit <- mxRun(two_dim_model)

## ----message=FALSE, warning=FALSE, results="hide"-----------------------------
spi_tree <- semtree(
  model = two_dim_fit,
  data = spi_analysis,
  control = semtree_control(method="score")
)


## ----message=FALSE, warning=FALSE, results="hide"-----------------------------
bonf_ctrl <- semtree.control(bonferroni = TRUE, 
                             alpha = 0.05,
                             method = "score")
spi_tree_bonf <- semtree(
  model = two_dim_fit,
  data = spi_analysis,
  control = bonf_ctrl
)

## -----------------------------------------------------------------------------
plot(spi_tree)
plot(spi_tree_bonf)