bspam: Passage-level Calibration & Scoring for ORF Data 1

Akihito Kamata

2026-10-27

Overview

This document is an in-depth tutorial for the basic use of the bspam R package for ORF data.

Installation

Please refer to the README document to install the bspam package and other related packages: https://github.com/kamataak/bspam

Load Packages

First, we load required packages, including the bspam and tidyverse packages. We also load the example data set passage2, which comes with the bspam package.

library(tidyverse)
library(bspam)
data("passage2")

Data Exploration

We use an example data included in the package for demonstration purposes.The data set passage2 is a passage-level student data on ORF assessment and consisted of reading accuracy and time data for 12 passages from 85 students.Although the 85 students were assigned to all 12 passages, the number of passages read by the 85 students varied from 2 to 12 passages. The number of students per passage were between 59 to 79. This is a small subset of the data collected by Nese and Kamata (2014-2018).

The seven variables in the data set includes: - id.student: individual ID - occasion: timing of the assessment - grade: grade level of the student - id.passage: passage ID - numwords.pass: the number of words in the passage - wrc: the number of correctly read words - sec: the time took to read the passage in seconds

Some descriptive statistics are as follows.

passage2 %>% 
  group_by(id.passage) %>% 
  summarise(n.students = n(),
            numwords.pass = unique(numwords.pass),
            mean.wrc = mean(wrc),
            mean.sec = mean(sec)) %>% 
  ungroup() %>% 
  mutate(rate.wrc = mean.wrc/numwords.pass,
         mean.sec.per10 = mean.sec/numwords.pass*10) %>% 
  knitr::kable(digits = 2)
id.passage n.students numwords.pass mean.wrc mean.sec rate.wrc mean.sec.per10
22078 74 50 45.76 30.35 0.92 6.07
23019 64 88 80.47 49.86 0.91 5.67
32004 71 52 47.28 32.54 0.91 6.26
32010 76 52 47.64 33.45 0.92 6.43
32012 76 49 44.24 29.25 0.90 5.97
32015 71 51 46.04 33.25 0.90 6.52
32016 71 50 45.17 40.48 0.90 8.10
33003 59 86 74.46 52.89 0.87 6.15
33037 72 73 66.03 42.39 0.90 5.81
33040 67 76 68.25 45.13 0.90 5.94
42055 79 53 46.47 34.77 0.88 6.56
43030 67 89 79.84 55.32 0.90 6.22 
passage2 %>% 
  group_by(id.student) %>% 
  summarise(n.passages = n(),
            tot.wrc = sum(wrc),
            tot.sec = sum(sec),
            obs.wcpm = tot.wrc/tot.sec*60) %>% 
  ungroup() %>% 
  knitr::kable(digits = 2)
id.student n.passages tot.wrc tot.sec obs.wcpm
1020 12 754 432.98 104.49
1037 12 653 413.74 94.70
1154 12 629 469.75 80.34
1164 9 526 303.00 104.16
1174 12 715 600.55 71.43
1184 12 703 541.44 77.90
1322 11 551 439.22 75.27
1332 10 596 385.35 92.80
1342 12 673 545.93 73.97
1492 11 677 312.24 130.09
1503 12 748 626.37 71.65
1513 12 753 353.77 127.71
1523 12 744 468.07 95.37
1565 12 657 342.02 115.26
157 12 757 445.69 101.91
1575 5 212 249.14 51.06
158 12 754 436.01 103.76
171 12 708 394.38 107.71
1737 10 632 435.88 87.00
1747 8 492 345.00 85.57
1757 2 67 180.24 22.30
1767 12 708 353.02 120.33
177 12 745 391.99 114.03
1777 11 565 375.02 90.40
181 10 626 493.86 76.05
1839 11 700 364.55 115.21
1849 7 324 287.02 67.73
1859 9 473 383.58 73.99
1890 2 81 75.69 64.21
1910 7 298 193.87 92.23
203 12 758 347.51 130.87
2033 11 680 358.08 113.94
2043 12 704 596.21 70.85
2079 12 744 386.08 115.62
2089 10 508 589.60 51.70
2125 8 399 338.33 70.76
2135 11 625 438.09 85.60
2145 12 676 553.79 73.24
2155 10 509 304.42 100.32
2185 9 432 561.72 46.14
2215 11 584 427.12 82.04
2235 7 310 202.04 92.06
2245 11 544 486.91 67.03
2255 10 510 414.87 73.76
2265 10 535 457.16 70.22
2275 12 716 396.16 108.44
2613 8 405 472.97 51.38
2648 12 714 444.02 96.48
2658 6 312 247.82 75.54
2699 11 627 497.16 75.67
2708 12 741 390.18 113.95
271 12 733 396.20 111.00
2718 12 715 578.50 74.16
281 12 750 318.16 141.44
289 12 749 391.22 114.87
291 11 578 278.10 124.70
305 12 754 357.13 126.68
309 11 665 441.40 90.39
335 7 320 452.56 42.43
369 6 251 238.58 63.12
389 12 741 310.34 143.26
396 12 703 353.77 119.23
409 11 516 343.47 90.14
412 7 368 397.09 55.60
58 5 243 311.21 46.85
592 3 129 226.55 34.16
60 3 161 179.21 53.90
629 5 265 214.50 74.13
633 2 117 127.09 55.24
645 12 746 495.97 90.25
66 9 423 504.04 50.35
698 12 661 542.18 73.15
700 11 675 415.52 97.47
705 12 669 659.79 60.84
709 12 726 249.25 174.76
715 12 707 258.01 164.41
830 11 662 320.51 123.93
838 11 574 318.96 107.98
848 11 596 313.08 114.22
909 11 680 362.76 112.47
919 7 404 400.43 60.53
954 12 745 397.61 112.42
964 8 315 667.88 28.30
968 12 747 469.16 95.53
992 9 450 557.02 48.47

Passage Calibration

In this section, ORF passage calibration, namely, model parameter estimation is demonstrated. We will use the fit.model() function with the default estimation method, which is the Monte Carlo EM algorithm described in Potgieter et al. (2025).

MCEM_run1 <- fit.model(data = passage2,
                       person.id = "id.student",
                       task.id = "id.passage",
                       max.counts = "numwords.pass",
                       obs.counts = "wrc",
                       time = "sec",
                       k.in = 5,
                       reps.in = 50,
                       est = "mcem")
MCEM_run1
#> $task.param
#> # A tibble: 12 × 6
#>        a     b alpha  beta task.id max.counts
#>    <dbl> <dbl> <dbl> <dbl>   <dbl>      <int>
#>  1 0.351 -1.42  6.64  1.74   22078         50
#>  2 0.302 -1.36  8.22  1.74   23019         88
#>  3 0.418 -1.37  6.98  1.80   32004         52
#>  4 0.490 -1.52  7.35  1.82   32010         52
#>  5 0.343 -1.35  5.45  1.74   32012         49
#>  6 0.605 -1.42  6.69  1.82   32015         51
#>  7 0.419 -1.42  2.55  1.98   32016         50
#>  8 0.496 -1.07  8.68  1.85   33003         86
#>  9 0.692 -1.42  7.87  1.74   33037         73
#> 10 0.634 -1.37  6.49  1.78   33040         76
#> 11 0.533 -1.28  5.59  1.79   42055         53
#> 12 0.513 -1.32  9.47  1.86   43030         89
#> 
#> $hyper.param
#> # A tibble: 1 × 2
#>   vartau   rho
#>    <dbl> <dbl>
#> 1 0.0979 0.461
#> 
#> attr(,"class")
#> [1] "fit.model"
MCEM_run1$task.param %>% knitr::kable(digits = 2)
a b alpha beta task.id max.counts
0.35 -1.42 6.64 1.74 22078 50
0.30 -1.36 8.22 1.74 23019 88
0.42 -1.37 6.98 1.80 32004 52
0.49 -1.52 7.35 1.82 32010 52
0.34 -1.35 5.45 1.74 32012 49
0.60 -1.42 6.69 1.82 32015 51
0.42 -1.42 2.55 1.98 32016 50
0.50 -1.07 8.68 1.85 33003 86
0.69 -1.42 7.87 1.74 33037 73
0.63 -1.37 6.49 1.78 33040 76
0.53 -1.28 5.59 1.79 42055 53
0.51 -1.32 9.47 1.86 43030 89 
MCEM_run1$hyper.param %>% knitr::kable(digits = 3)
vartau rho
0.098 0.461
Interpretations

[Add some interpretations of the estimated quantities.]

By default, the standard errors for the model parameters are not estimated. This will allow one to increase the number of Monte-Carlo iterations reps.in to improve the quality of the model parameter estimates, while minimizing the computation time. The number of reps.in should be 50 to 100 in realistic calibrations. Also note that k.in is the number of imputations for the MCEM algorithm, where the default is 5. Note that SE’s for model parameters are not required for running the scoring() function to estimate latent factor scores and WCPM scores in the next step. If standard errors for model parameters are desired, an additional argument se = "analytic" or se = "bootstrap" needs to be added to the fit.model() function as follows.

MCEM_run2 <- fit.model(data = passage2,
                       person.id = "id.student",
                       task.id = "id.passage",
                       max.counts = "numwords.pass",
                       obs.counts = "wrc",
                       time = "sec",
                       k.in = 5,
                       reps.in = 50,
                       est = "mcem",
                       se = "analytic")
MCEM_run2$task.param %>% knitr::kable(digits = 2)
MCEM_run2$hyper.param %>% knitr::kable(digits = 3)

Scoring: Estimating Latent Factor Scores and WCPM Scores 1

In order to estimate latent factor scores and/or WCPM scores, the scoring() function needs to be run.

Note that we use the output object MCEM_run from the previous passage calibration phase. By default (type = "general"), only factor scores (\(\theta\) and \(\tau\)) and their standard errors are reported. By providing an argument type = "orf", WCPM scores and their standard errors are also reported. Also, WCPM scores and/or factor scores, and their standard errors will be estimated for all cases in the data by default.

There are several estimator options and standard error estimation options: maximum likelihood (MLE), maximum a posteriori (MAP), expected a posteriori (EAP), and fully Bayesian approach. In this example, maximum a priori (MAP) estimators for scoring and analytic approach for estimating standard errors are used.

Scoring for All Cases

To estimate WCPM scores for all observations in the data set passage2, we use the scoring() function as follows. While there are four options for the estimator, this example uses the maximum a posteriori (MAP) estimator. This option is specified by the est = "map" argument. Also, the standard errors are estimated by the formula derived by the Delta method, which is specified by the se = "analytic" argument. Additionally, the type = "orf" argument will estimate the model-based WCPM scores and their standar errors, in addition to the speed and acuracy latent factor scores and their standard errors.

WCPM_all <- scoring(calib.data=MCEM_run, 
                    data = passage2,
                    person.id = "id.student",
                    occasion = "occasion",
                    group = "grade",
                    task.id = "id.passage",
                    max.counts = "numwords.pass",
                    obs.counts = "wrc",
                    time = "sec",
                    est = "map", 
                    se = "analytic",
                    type="orf")
summary(WCPM_all)
Interpretations

[Add some interpretations of the estimated quantities.]