## ----echo = FALSE, message = FALSE, output = "hide"---------------------------
library(knitr)
library(kableExtra)
library(PROsetta)
library(dplyr)

## -----------------------------------------------------------------------------
d <- loadData(
  response  = response_dep,
  itemmap   = itemmap_dep,
  anchor    = anchor_dep
)

## -----------------------------------------------------------------------------
head(response_dep)

## -----------------------------------------------------------------------------
head(itemmap_dep)

## -----------------------------------------------------------------------------
head(anchor_dep)

## ----freq---------------------------------------------------------------------
freq_table <- runFrequency(d)
head(freq_table)

## ----fig.align = "center", fig.width = 7, fig.height = 7----------------------
plot(d, scale = "combined", title = "Combined scale")

## ----eval = FALSE-------------------------------------------------------------
# plot(d, scale = 1, title = "Scale 1") # not run
# plot(d, scale = 2, title = "Scale 2") # not run

## ----desc---------------------------------------------------------------------
desc_table <- runDescriptive(d)
head(desc_table)

## ----alpha--------------------------------------------------------------------
classical_table <- runClassical(d)
summary(classical_table$alpha$combined)

## ----alpha2-------------------------------------------------------------------
classical_table <- runClassical(d, scalewise = TRUE)
summary(classical_table$alpha$`2`)

## ----omega, eval = FALSE------------------------------------------------------
# classical_table <- runClassical(d, scalewise = TRUE, omega = TRUE)
# classical_table$omega$combined # omega values for combined scale
# classical_table$omega$`1`      # omega values for each scale, created when scalewise = TRUE
# classical_table$omega$`2`      # omega values for each scale, created when scalewise = TRUE

## ----eval = FALSE-------------------------------------------------------------
# classical_table <- runClassical(d, scalewise = TRUE, omega = TRUE, nfactors = 5) # not run

## ----cfa, results = "hide"----------------------------------------------------
out_cfa <- runCFA(d, scalewise = TRUE)

## ----eval = FALSE-------------------------------------------------------------
# out_cfa <- runCFA(d, scalewise = TRUE, std.lv = TRUE) # not run

## -----------------------------------------------------------------------------
out_cfa$combined
out_cfa$`1`
out_cfa$`2`

## -----------------------------------------------------------------------------
lavaan::summary(out_cfa$combined, fit.measures = TRUE, standardized = TRUE, estimates = FALSE)

## ----eval = FALSE-------------------------------------------------------------
# lavaan::summary(out_cfa$`1`, fit.measures = TRUE, standardized = TRUE, estimates = FALSE) # not run
# lavaan::summary(out_cfa$`2`, fit.measures = TRUE, standardized = TRUE, estimates = FALSE) # not run

## ----calib, results = "hide", error = TRUE, message = FALSE-------------------
try({
out_calib <- runCalibration(d, technical = list(NCYCLES = 1000))
})

## ----calib2, results = "hide", error = TRUE, message = FALSE------------------
try({
out_calib <- runCalibration(d, technical = list(NCYCLES = 10))
})

## ----mirt_plot, eval = FALSE, results = "hide"--------------------------------
# mirt::itemfit(out_calib, empirical.plot = 1)
# mirt::itemplot(out_calib, item = 1, type = "info")
# mirt::itemfit(out_calib, "S_X2", na.rm = TRUE)

## ----fig.align = "center", fig.width = 7, fig.height = 7----------------------
plotInfo(
  out_calib, d,
  scale_label = c("PROMIS Depression", "CES-D", "Combined"),
  color = c("blue", "red", "black"),
  lty = c(1, 2, 3)
)

## ----fixedpar, results = "hide", message = FALSE------------------------------
link_fixedpar <- runLinking(d, method = "FIXEDPAR")

## -----------------------------------------------------------------------------
head(link_fixedpar$ipar_linked)

## -----------------------------------------------------------------------------
link_fixedpar$mu_sigma

## ----sl, results = "hide", error = TRUE, message = FALSE----------------------
try({
link_sl <- runLinking(d, method = "SL", technical = list(NCYCLES = 1000))
link_sl
})

## -----------------------------------------------------------------------------
head(link_sl$ipar_linked)

## -----------------------------------------------------------------------------
link_sl$constants

## ----results = "hide", message = FALSE----------------------------------------
link_cp <- runLinking(d, method = "CP")

## -----------------------------------------------------------------------------
head(link_cp$ipar_linked)

## -----------------------------------------------------------------------------
link_cp$mu_sigma

## -----------------------------------------------------------------------------
rsss_fixedpar <- runRSSS(d, link_fixedpar)

## -----------------------------------------------------------------------------
head(round(rsss_fixedpar$`2`, 3))

## ----eqp_raw, results = "hide", message = FALSE-------------------------------
rsss_equate_raw <- runEquateObserved(
  d,
  scale_from = 2, # CES-D (scale to be linked)
  scale_to = 1,   # PROMIS Depression (anchor scale)
  eq_type = "equipercentile", smooth = "loglinear"
)

## -----------------------------------------------------------------------------
head(rsss_equate_raw$concordance)

## ----eqp_dir, results = "hide", message = FALSE-------------------------------
rsss_equate_tscore <- runEquateObserved(
  d,
  scale_from = 2, # CES-D (scale to be linked)
  scale_to = 1,   # PROMIS Depression (anchor scale)
  type_to = "tscore",
  rsss = rsss_fixedpar, # used to convert PROMIS Depression (anchor scale) raw to T
  eq_type = "equipercentile", smooth = "loglinear"
)

## -----------------------------------------------------------------------------
head(rsss_equate_tscore$concordance)

## ----fig.align = "center", fig.width = 7, fig.height = 7----------------------
plot(
  rsss_fixedpar$`2`$raw_2,
  rsss_fixedpar$`2`$tscore,
  xlab = "CES-D (Scale to be linked), raw score",
  ylab = "PROMIS Depression (Anchor scale), T-score",
  type = "l", col = "blue")
lines(
  rsss_equate_tscore$concordance$raw_2,
  rsss_equate_tscore$concordance$tscore_1,
  lty = 2, col = "red")
grid()
legend(
  "topleft",
  c("Fixed-Parameter Calibration", "Equipercentile Linking"),
  lty = 1:2, col = c("blue", "red"), bg = "white"
)

## -----------------------------------------------------------------------------
scores <- getScaleSum(d, 2)
head(scores)

## ----message = FALSE----------------------------------------------------------
theta_promis <- getTheta(data = d, ipar = link_fixedpar$ipar_linked, scale = 1)$theta
head(theta_promis)

## -----------------------------------------------------------------------------
t_promis_pattern <- data.frame(
  prosettaid = theta_promis$prosettaid,
  t_promis_pattern = round(theta_promis$theta_eap * 10 + 50, 1)
)
head(t_promis_pattern)

## -----------------------------------------------------------------------------
scores <- merge(scores, t_promis_pattern, by = "prosettaid")
head(scores)

## ----message = FALSE----------------------------------------------------------
rsss_fixedpar <- data.frame(
  raw_2 = rsss_fixedpar[["2"]]$raw_2,
  t_promis_rsss_fixedpar = round(rsss_fixedpar[["2"]]$tscore, 1)
)
scores <- merge(scores, rsss_fixedpar, by = "raw_2")
head(scores)

## ----message = FALSE----------------------------------------------------------
rsss_eqp <- data.frame(
  raw_2 = rsss_equate_tscore$concordance$raw_2,
  t_promis_rsss_eqp = round(rsss_equate_tscore$concordance$tscore_1, 1)
)
scores <- merge(scores, rsss_eqp, by = "raw_2")
head(scores)

## -----------------------------------------------------------------------------
c_fixedpar <- compareScores(
  scores$t_promis_pattern, scores$t_promis_rsss_fixedpar)
c_eqp <- compareScores(
  scores$t_promis_pattern, scores$t_promis_rsss_eqp)

stats           <- rbind(c_fixedpar, c_eqp)
rownames(stats) <- c("Fixed-parameter calibration", "Equipercentile")
stats