Scale linking with PROsetta package

In PROsetta: Linking Patient-Reported Outcomes Measures

Introduction

This vignette explains how to perform scale linking with the PROsetta package. By way of illustration, we replicate the linking of the Center for Epidemiologic Studies Depression Scale (CES-D) to the PROMIS Depression metric as described in Choi, Schalet, Cook, and Cella (2014).

library(knitr)
library(kableExtra)
library(PROsetta)
library(dplyr)

Load datasets

First step is to load the input datasets comprised of three tables with loadData(). The PROMIS Depression – CES-D linking data are included in the PROsetta package directory under the folder labeled data-raw.

d <- data_dep

fp <- system.file("data-raw", package = "PROsetta")
d <- loadData(
  response  = "dat_DeCESD_v2.csv",
  itemmap   = "imap_DeCESD.csv",
  anchor    = "anchor_DeCESD.csv",
  input_dir = fp)

• response: Contains item response data from both instruments. You can supply a .csv filename or a data frame. In this example, we supply a .csv filename dat_DeCESD_v2.csv.
• itemmap: Specifies which items belong to which instruments. Can be a .csv filename or a data frame.
• anchor: Contains tem parameters for anchor items (e.g., PROMIS Depression). Can be a .csv filename or a data frame.
• input_dir: (Optional) The path of the directory to look for the input .csv files.

Response data

The response data contains individual item responses on both instruments (i.e., 28 PROMIS Depression items followed by 20 CES-D items). The data table should include the following columns.

• prosettaid: The person ID of the respondents (N = 747). This column does not have to be named prosettaid but should not conflict with other data tables (item map and anchor).
• Other columns should include the item response fields with their unique item IDs as column names. The item names should match the item_id column in both the item map and anchor files.

Run the following code, for example, to open the response data in edit mode.

file.edit(system.file("data-raw", "dat_DeCESD_v2.csv", package = "PROsetta"))

Item map data

The item map data requires the following columns.

• item_id: Contains the unique ID of the items. The name of this column does not have to be item_id but should be consistent with the item ID column in the anchor table. The IDs in this column should match the column names in the response data.
• instrument: Numerals (1 or 2) indicating to which of the two instruments the items belong (e.g., 1 = PROMIS Depression; 2 = CES-D)
• item_order: The sequential position of the items in the combined table (e.g., 1, 2, 3, ..., 28, ..., 48)
• item_name: Secondary labels for the items
• ncat: The number of response categories by item

Run the following code to open the item map data in edit mode.

file.edit(system.file("data-raw", "imap_DeCESD.csv", package = "PROsetta"))

Anchor data

The anchor data contains the item parameters for the anchor scale (e.g., PROMIS Depression) and requires the following columns.

• item_order: The sequential position of the items in the anchor scale (e.g., 1, 2, 3, ..., 28)
• item_id: The unique ID of the anchor items. The name of this column does not have to be item_id but should be consistent with the item ID column in the item map table The IDs in this column should refer to the specific column names in the response data.
• a: The slope parameter value for each anchor item
• cb1, cb2, ...: The category boundary parameter values for each anchor item

Run the following code to open the anchor data in edit mode.

file.edit(system.file("data-raw", "anchor_DeCESD.csv", package = "PROsetta"))

Descriptive analysis

Basic descriptive statistics

The frequency distribution of each item in the response data is obtained by runFrequency().

freq_table <- runFrequency(d)
head(freq_table)

The frequency distribution of the summed scores for the combined scale can be plotted as a histogram with plot(). The required argument is a PROsetta_data object created with loadData(). The optional scale argument specifies for which scale the summed score should be created. Setting scale = 'combined' plots the summed score distribution for the combined scale.

plot(d, scale = "combined", title = "Combined scale")

The user can also generate the summed score distribution for the first or second scale by specifying scale = 1 or scale = 2.

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

Basic descriptive statistics are obtained for each item by runDescriptive().

desc_table <- runDescriptive(d)
head(desc_table)

Classical reliability analysis

Classical reliability statistics can be obtained by runClassical(). By default, the analysis is performed for the combined scale.

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

The user can set scalewise = TRUE to request an analysis for each scale separately in addition to the combined scale.

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

Specifying omega = TRUE returns the McDonald's $\omega$ coefficients as well.

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

Additional arguments can be supplied to runClassical() to pass onto psych::omega().

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

Dimensionality analysis

Dimensionality analysis is performed with CFA by runCFA(). Setting scalewise = TRUE performs the dimensionality analysis for each scale separately in addition to the combined scale.

out_cfa <- runCFA(d, scalewise = TRUE)

runCFA() calls for lavaan::cfa() internally and can pass additional arguments onto it.

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

The CFA result for the combined scale is stored in the combined slot, and if scalewise = TRUE, the analysis for each scale is also stored in each numbered slot.

out_cfa$combined
out_cfa$`1`
out_cfa$`2`

CFA fit indices can be obtained by using summary() from the lavaan package. For the combined scale:

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

and also for each scale separately:

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

Item parameter calibration

runCalibration() performs IRT calibration without anchoring. runCalibration() calls mirt::mirt() internally. Additional arguments can be passed onto mirt, e.g., to increase the number of EM cycles to 1000, as follows:

out_calib <- runCalibration(d, technical = list(NCYCLES = 1000))

In case of nonconvergence, runCalibration() explicitly raises an error and does not return its results.

out_calib <- runCalibration(d, technical = list(NCYCLES = 10))

Also, specify fixedpar = TRUE to perform fixed parameter calibration using the anchor data.

out_calib <- runCalibration(d, fixedpar = TRUE)

The output object from runCalibration() can be used to generate additional output with functions from the mirt package.

Use coef() to extract item parameters:

mirt::coef(out_calib, IRTpars = TRUE, simplify = TRUE)

and also other commonly used functions:

mirt::itemfit(out_calib, empirical.plot = 1)
mirt::itemplot(out_calib, item = 1, type = "info")
mirt::itemfit(out_calib, "S_X2", na.rm = TRUE)

Scale information functions can be plotted with plotInfo. The two required arguments are an output object from runCalibration() and a PROsetta object from loadData(). The additional arguments specify the labels, colors, and line types for each scale and the combined scale. The last values in arguments scale_label, color, lty represent the values for the combined scale.

plotInfo(
  out_calib, d,
  scale_label = c("PROMIS Depression", "CES-D", "Combined"),
  color = c("blue", "red", "black"),
  lty = c(1, 2, 3))

Scale aligning

runLinking() performs item parameter linking based on the anchor item parameters supplied in the anchor table. Two linking, or more specifically scaling aligning, methods currently available are fixed-parameter calibration and linear transformation. Fixed-parameter calibration estimates the item parameters for the non-anchor items on the metric defined by the anchor items, while fixing the item parameters for the anchor items to their supplied anchor values. The linear transformation methods determine linear transformation constants, i.e., a slope and an intercept, to transform freely estimated item parameters to the metric defined by the anchor items.

Fixed parameter calibration method

Scale aligning through fixed parameter calibration is performed by setting method = "FIXEDPAR". The linked parameters are stored in the $ipar_linked slot.

out_link_fixedpar <- runLinking(d, method = "FIXEDPAR")

out_link_fixedpar$ipar_linked

Linear transformation methods

Scale aligning through linear transformation is performed by setting the method argument to one of the following options:

• MM (Mean-Mean)
• MS (Mean-Sigma)
• HB (Haebara)
• SL (Stocking-Lord)

Arguments supplied to runLinking are passed onto mirt::mirt() internally. In case of nonconvergence in the free calibration step, runLinking() explicitly raises an error and does not return its results.

out_link_sl <- runLinking(d, method = "SL", technical = list(NCYCLES = 1000))
out_link_sl

The item parameter estimates linked to the anchor metric are stored in the $ipar_linked slot.

out_link_sl$ipar_linked

Transformation constants (A = slope; B = intercept) for the specified linear transformation method are stored in the $constants slot.

out_link_sl$constants

Obtaining scaled scores

From the item parameter estimates transformed to the anchor metric, raw-score-to-scale-score (rsss) crosswalk tables can be generated by runRSSS().

The output from runRSSS() includes three crosswalk tables (labeled as 1, 2, and combined), one for each scale and the third one for the combined scale. Each table contains raw summed scores and corresponding scaled scores, including summed score EAP estimate, T-scores corresponding to the EAP estimates, as well as expected summed scores (i.e., true scores) for each scale from the EAP estimates.

rsss_fixedpar <- runRSSS(d, out_link_fixedpar)
rsss_sl       <- runRSSS(d, out_link_sl)
round(rsss_fixedpar$`2`, 3)

The columns in the crosswalk tables include:

• raw_1: raw summed score in Scale 1 (also raw_2 for Scale 2 and raw_3 for the combined)
• tscore: T-score corresponding to each summed score
• tscore_se: standard error associated with each T-score
• eap: summed score EAP equivalent for each raw summed score
• eap_se: standard error associated with each EAP estimate
• escore_1: expected summed score (true score) for Scale 1 given the EAP estimate
• escore_2: expected summed score (true score) for Scale 2 given the EAP estimate
• escore_combined: expected summed score (true score) for the combined scale given the EAP estimate

Equipercentile method: raw-raw

Equipercentile linking of observed summed scores is performed by runEquateObserved().

Cases with missing responses are removed to be able to generate correct summed scores in concordance tables.

This function requires four arguments:

• scale_from: numeric index of the scale (as specified in the item map) to be linked
• scale_to: numeric index of the scale (as specified in the item map) to serve as the anchor
• eq_type: the type of equating to be performed, equipercentile for this example. See ?equate::equate for details.
• smooth: the type of presmoothing to perform

By default, runEquateObserved() performs raw-raw equipercentile linking. In this example, each raw summed score in Scale 2 (CES-D, ranging from 20 to 80) is linked to a raw summed score equivalent in Scale 1 (PROMIS Depression, rangeing from 28 to 140) with loglinear presmoothing.

out_equate <- runEquateObserved(
  d, scale_from = 2, scale_to = 1,
  eq_type = "equipercentile", smooth = "loglinear")

The crosswalk table can be obtained from the concordance slot:

out_equate$concordance

Equipercentile method: raw-tscore

Raw summed scores can be linked to scaled scores (e.g., T-scores) directly by specifying type_to = 'tscore' in runEquateObserved(). In the following example, we map the raw summed scores from Scale 2 (CES-D, ranging from 20 to 80) onto the T-score equivalents in Scale 1 (PROMIS Depression, mean = 50 and SD = 10).

out_equate_tscore <- runEquateObserved(
  d, scale_from = 2, scale_to = 1,
  type_to = "tscore", rsss = rsss_fixedpar,
  eq_type = "equipercentile", smooth = "loglinear")

Again, the crosswalk table can be retrieved from the concordance slot:

out_equate_tscore$concordance

In what follows we display the linking relation obtained from the equipercentile method and compare it to that from the fixed-parameter calibration method.

plot(
  rsss_fixedpar$`2`$raw_2,
  rsss_fixedpar$`2`$tscore,
  xlab = "CES-D Summed Score",
  ylab = "PROMIS Depression T-score",
  type = "l", col = "blue")
lines(
  out_equate_tscore$concordance$raw_2,
  out_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"
)

Evaluation of linking results

The linking results produced so far are now evaluated. More specifically, we assess how closely the CES-D summed scores linked to the PROMIS Depression T-scores match the actual PROMIS Depression T-scores observed in the present linking sample. Should we have set aside a validation sample, we would have performed this evaluation on that sample.

Raw scores from Scale 2

To begin with, we create an object scores using getScaleSum() to contain raw summed scores on Scale 2 (i.e., CES-D). NA will result for any respondents with one or more missing responses on Scale 2. We could also create a summed score variable for Scale 1 using the same function, e.g., getScaleSum(d, 1).

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

EAP estimates based on item responses patterns on Scale 1

We obtain EAP estimates of theta on Scale 1 (i.e., PROMIS Depression) based on item response patterns using the getTheta() function. The first argument of the function is a data object of PROsetta class, which we created earlier with loadData(). The second argument specifies the item parameter estimates to be used for the EAP estimation. Here, we use the item parameter estimates previously obtained from the fixed-parameter calibration, out_link_fixedpar$ipar_linked. The third argument scale = 1 specifies the scale to be scored (i.e., PROMIS Depression). These EAP estimates are based on the item responses actually observed on PROMIS Depression and will serve as the reference when we assess the CES-D scores liked to PROMIS Depression derived from various methods.

eap_promis <- getTheta(d, out_link_fixedpar$ipar_linked, scale = 1)$theta
head(eap_promis)

The EAP estimates for PROMIS Depression will be converted to T-scores using a linear transformation.

t_promis <- data.frame(
  prosettaid = eap_promis$prosettaid,
  t_promis = round(eap_promis$theta_eap * 10 + 50, 1)
)
head(t_promis)

We then merge the PROMIS Depression T-scores with the raw summed scores for CES-D calculated in the previous step.

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

Now we are going to generate T-scores linked to PROMIS Depression using only item responses on Scale 2 (CES-D). These T-scores linked to PROMIS Depression can be generated in different ways as:

• EAP estimates based on item response patterns on Scale 2
• EAP estimates based on summed scores on Scale 2
• Equipercentile equivalents based on summed scores on Scale 2

The first two ways are based on the CES-D item parameters linked to the PROMIS Depression metric.

EAP estimates based on item responses patterns on Scale 2

First, we get EAP estimates based on item response patterns on Scale 2 using the CES-D item parameters linked to the PROMIS Depression metric (via fixed-parameter calibration). We then linearly transform the EAP estimates to T-scores and add the T-scores (t_cesd_pattern) to the data frame object scores.

eap_cesd <- getTheta(d, out_link_fixedpar$ipar_linked, scale = 2)$theta
t_cesd_pattern <- data.frame(
  prosettaid = eap_cesd$prosettaid,
  t_cesd_pattern = round(eap_cesd$theta_eap * 10 + 50, 1)
)
scores <- merge(scores, t_cesd_pattern, by = "prosettaid")
head(scores)

EAP estimates based on summed scores on Scale 2

Second, we use the raw-score-to-scale-score (RSSS) crosswalk table obtained above using summed score EAP estimation to map each raw summed score on Scale 2 onto a T-score on the PROMIS Depression metric, t_cesd_rsss.

rsss_eap <- data.frame(
  raw_2 = rsss_fixedpar$`2`$raw_2,
  t_cesd_rsss = round(rsss_fixedpar$`2`$tscore, 1)
)
scores <- merge(scores, rsss_eap, by = "raw_2")
head(scores)

Equipercentile linking of summed scores on Scale 2

Third, we use the concordance table from equipercentile linking to map each raw summed score on Scale 2 onto a T-score on the PROMIS Depression metric, t_cesd_eqp.

rsss_eqp <- data.frame(
  raw_2 = out_equate_tscore$concordance$raw_2,
  t_cesd_eqp = round(out_equate_tscore$concordance$tscore_1, 1)
)
scores <- merge(scores, rsss_eqp, by = "raw_2")
head(scores)

Comparison of equated and observed T-scores

Finally, use compareScores() to compare the obtained T-scores.

# Reference score: IRT pattern scoring of Scale 1
c_pattern <- compareScores(
  scores$t_promis, scores$t_cesd_pattern) ## IRT response pattern EAP to T-score
c_rsss <- compareScores(
  scores$t_promis, scores$t_cesd_rsss)    ## IRT summed score EAP to T-score
c_eqp <- compareScores(
  scores$t_promis, scores$t_cesd_eqp)     ## Equipercentile summed score to T-score

stats           <- rbind(c_pattern, c_rsss, c_eqp)
rownames(stats) <- c("IRT Pattern", "IRT RSSS", "Equipercentile")
stats

PROsetta documentation built on Feb. 16, 2023, 9:14 p.m.