這份是入門篇（vignette("kits-flow", package = "kits")）的續篇。 入門篇是流程的第 0 到 6 步，這份是第 7 到 10 步：
| 步驟 | 函數 | 做什麼 |
|---|---|---|
| 7 題目分析與報告的圖 | kits_dif()、kits_dif_combine()、kits_item_table() 和它的 plot() |
性別 DIF 和 LR 檢定、通過月齡、題目的位置、成長曲線和階梯 |
| 8 新的波次 | kits_drift()、kits_item(fixed_params = )、plot(kits_vertical(), type = "median") |
偏移檢查、固定參數校準、各波次的中位數 |
| 9 演算法和程式 | kits_algorithm()、algo_ 函數、inst/scripts/ |
每一步的公式、函數和文獻；拆開的演算法；教學腳本 |
| 10 量尺的性質 | kits_convert() 的 se_link、kits_vertical()、kits_precision()、kits_expected_score()、kits_link_check()、summary(link, theta = )、coef(item, rp = )、kits_item_table(rp = ) |
連結誤差、垂直量尺、精準度、原始分對照、連結的檢查、反應機率、品質管制（Kolen & Brennan, 2014） |
看不懂的詞，查套件 README 的「名詞表」，或 ?"kits-package" 的 Glossary（英文）。
library(kits)
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(tidyr)
library(purrr)和入門篇相同的設計：語言 14 題、0–2 分；M03 答第 1–8 題（6、12、24 個月），M36 答第 5–12 題（36、48 個月）。 M03 另外有動作 6 題，用來一起接到 D。 這裡只重做一次，不印出來；每一步的說明在入門篇。 亂數的順序和入門篇不同，所以數字會差一點。
set.seed(11)
betas <- lapply(seq(-3, 5, length.out = 14), function(b) b + c(-0.6, 0.6))
names(betas) <- sprintf("lan%02d", 1:14)
persons <- tibble(cohort = rep(c("M03", "M36"), c(900, 600)),
wave_m = rep(c(6, 12, 24, 36, 48), each = 300)) |>
mutate(visit_id = sprintf("%s_%04d", cohort, row_number()), .by = cohort) |>
mutate(age_m = wave_m + runif(n(), -1, 1),
sex = sample(c("boy", "girl"), n(), replace = TRUE))
theta <- 2 * log(persons$age_m / 12) + rnorm(nrow(persons), 0, 0.9)
m03 <- persons$cohort == "M03"
X03 <- kits_sim(theta[m03], betas[1:8], prefix = "M03_")
X36 <- kits_sim(theta[!m03], betas[5:12], prefix = "M36_")
X03[sample(nrow(X03), 12), "lan08"] <- NA
motor_betas <- as.list(seq(-2.5, 2.5, length.out = 6))
names(motor_betas) <- sprintf("mot%02d", 1:6)
Xm03 <- kits_sim(theta[m03] + rnorm(sum(m03), 0, 0.5), motor_betas, prefix = "M03_")
# 1 校準
item03 <- kits_item(X03, scale = "lan M03")
item36 <- kits_item(X36, scale = "lan M36")
itemm03 <- kits_item(Xm03, scale = "motor M03")
# 2 世代連結
co <- kits_linking(item36, item03)
# 3 接到 D
pass <- function(item) filter(coef(item), item_score == 1)
milestone <- bind_rows(as.data.frame(pass(item03)), as.data.frame(pass(itemm03)))
anchors <- tibble(item_id = c("lan01", "lan02", "lan03", "lan06", "lan07", "mot01", "mot03", "mot05")) |>
left_join(select(milestone, item_id, thr), by = "item_id") |>
mutate(tau = 60 + 7 * thr + rnorm(n(), 0, 1.5)) |>
select(item_id, tau)
d_link <- kits_linking(list(lan = item03, motor = itemm03), anchors)
# 4 計分與換算
s03 <- kits_subj(X03, item03, persons = persons) |>
kits_convert(d_link)
s36 <- kits_subj(X36, item36, persons = persons) |>
kits_convert(co) |>
kits_convert(d_link)
lan <- rbind(s03, s36)
# 5 常模
use <- filter(lan, kits_norm_screen(lan)$normed)
norm <- kits_norms(use, degrees = c(3, 2, 1, 1))圖上的標題、軸、圖例都是英文。
每題 50% 通過的月齡，由第 7.3 節的 kits_item_table(data = ) 算。 兩個世代一起算，月齡範圍才夠寬。 先把兩個世代的作答疊成一張寬的表（沒答的題是 NA），再照 persons 的列接上月齡和波次。 kits_item_table() 收這種表（visit_id、age_m、wave_m，其他欄是題），也收 kits_resp()（入門篇第 0 節）。
answers <- bind_rows(as_tibble(X03, rownames = "visit_id"), as_tibble(X36, rownames = "visit_id"))
X_all <- persons |>
select(visit_id, age_m, wave_m) |>
left_join(answers, by = "visit_id")每一題做一個邏輯斯迴歸：分數至少 1（milestone step）算通過，通過對月齡。 age50 是 NA 的題：通過率在觀察到的月齡裡沒有跨過一半（太容易，或只有大孩子答）。 完全分開（小的孩子都沒通過、大的都通過）時，迴歸沒有有限的斜率；題目表的診斷列出這些題。
kits_dif()、kits_dif_combine()上面的資料沒有性別差異。 為了示範，另外模擬一份有性別差異的作答：同樣的題、同樣的能力，男孩一份、女孩一份。 女孩的步驟參數加上 shift：負的對女孩容易，正的對女孩難。 lan06 在 M03 對女孩容易，在 M36 對女孩難。
shift_betas <- function(b, shift) {
b[names(shift)] <- map2(b[names(shift)], shift, `+`)
b
}
girl03 <- shift_betas(betas[1:8], c(lan03 = -1, lan06 = -0.8, lan07 = 1))
girl36 <- shift_betas(betas[5:12], c(lan06 = 0.8, lan07 = 0.8, lan09 = -0.8, lan10 = 0.8, lan11 = -0.8))
Xs03 <- rbind(kits_sim(theta[m03], betas[1:8], prefix = "boy03_"),
kits_sim(theta[m03], girl03, prefix = "girl03_"))
Xs36 <- rbind(kits_sim(theta[!m03], betas[5:12], prefix = "boy36_"),
kits_sim(theta[!m03], girl36, prefix = "girl36_"))
sex03 <- rep(c("boy", "girl"), each = sum(m03))
sex36 <- rep(c("boy", "girl"), each = sum(!m03))kits_dif() 在一個世代裡比較男孩和女孩的題目參數（以這個領域的總分為條件）。 方法和 dexter 的 DIF() 相同（Bechger & Maris, 2015）：男女各自用 CML 校準，再比較相對的參數。 每題一個數字，單位是 logit；多點計分的題是各步驟的平均。 ref = "boy"（預設）時，正的對男孩難，也就是有利女孩。 第二個參數 group 可以是每一列一個值，或用 visit_id 當名字的向量。 性別不是寫成 boy、girl 時，用 ref 指定參照組，例如 kits_dif(X, sex, ref = "M")。
dif03 <- kits_dif(Xs03, sex03)
dif36 <- kits_dif(Xs36, sex36)
dif03
#> <kits_dif> groups by the vector `group`
#> Groups: boy (n = 900), girl (n = 900)
#> LR test (Andersen): LR = 294.97, df = 15, p = <2e-16
#> Largest |step difference|: 1.19 logits (lan07:1, girl minus boy, Wald z = 8.9)
#> Item DIF in logits: boy minus girl (positive = harder for boy, favors girl)
#> Overall DIF test: chi-square = 277.04, df = 15, p = <2e-16
#> item_id dif
#> 1 lan01 0.09
#> 2 lan02 0.05
#> 3 lan03 0.74
#> 4 lan04 -0.13
#> 5 lan05 0.01
#> 6 lan06 0.74
#> 7 lan07 -1.12
#> 8 lan08 -0.38
#> Checks: 2 notes (see summary() or kits_check())plot() 畫每題的 DIF 和 95% 信賴區間。 超過 cut（預設 0.3）的題上色：橘色有利女孩，藍色有利男孩。
plot(dif03)Overall test 檢定是不是有任何一題有 DIF。
同一次呼叫也做 Andersen（1973）的 LR 檢定：兩組各自校準，看步驟參數是不是一樣。 印出來的 LR test 一行就是它；數字在 attr(dif03, "lrtest")，每個步驟的差和 Wald z 在 attr(dif03, "diff")。 人數多時 p 值幾乎都很小，要看差的大小。 plot(type = "steps") 畫每個步驟在兩組的 beta，|z| 超過 1.96 的是橘色：
plot(dif03, type = "steps")group 也可以是 "raw"（每個作答組合裡，原始分的中位數分成高低兩組）或 "wave"（一個波次一組）。
kits_dif_combine() 把兩個世代合起來判斷（預設 cut = 0.3 logit）：
cut：none。cut，另一邊的方向相反：inconsistent。cohorts 是兩個世代的名字，用在欄名；沒給時用參數的名字（這裡會是 dif03、dif36）。
dif <- kits_dif_combine(dif03, dif36, cohorts = c("M03", "M36"))
dif
#> item_id dif_M03 dif_M36 dir_M03 dir_M36 DIF
#> 1 lan01 0.09 NA none <NA> none
#> 2 lan02 0.05 NA none <NA> none
#> 3 lan03 0.74 NA favors girls <NA> favors girls
#> 4 lan04 -0.13 NA none <NA> none
#> 5 lan05 0.01 0.34 none favors girls favors girls
#> 6 lan06 0.74 -0.70 favors girls favors boys inconsistent
#> 7 lan07 -1.12 -0.57 favors boys favors boys favors boys
#> 8 lan08 -0.38 0.02 favors boys none inconsistent
#> 9 lan09 NA 0.84 <NA> favors girls favors girls
#> 10 lan10 NA -0.78 <NA> favors boys favors boys
#> 11 lan11 NA 0.82 <NA> favors girls favors girls
#> 12 lan12 NA 0.03 <NA> none none模擬時 lan03、lan09、lan11 對女孩容易，lan07、lan10 對女孩難，lan06 在兩個世代方向相反。 這幾題都找到了。 lan05 和 lan08 沒有模擬 DIF，卻也超過了 cut。 lan05 是 M36 最容易的題，lan08 是 M03 最難的題：幾乎每個孩子都通過（或都沒通過），所以估得不準。 看 DIF 時，也要看這題在這個世代裡的位置和人數。 DIF 這一欄接到下面的題目表，圖用它上色。
kits_item_table()一題一列：milestone step 的 D、50% 通過的月齡（從 data）、常模中位數到達的月齡、DIF、題幹。 門檻表先用 kits_convert() 換到 D：
兩個世代都有的題，留第一個（這裡是 M03 的）。 表裡帶著常模，所以畫圖時不用再給。
thr03 <- coef(item03) |>
kits_convert(d_link)
thr36 <- coef(item36) |>
kits_convert(co) |>
kits_convert(d_link)
stems <- tibble(item_id = names(betas), stem = paste("Example item", names(betas)))
norm <- mutate(norm, domain = "lan")
item_tab <- kits_item_table(list(thr03, thr36), norm, data = X_all, dif = dif, stems = stems)
#> Item(s) given twice keep their first row: lan05, lan06, lan07, lan08.
item_tab |>
select(item_id, D, age50, age_reach, x, DIF) |>
head()
#> item_id D age50 age_reach x DIF
#> 1 lan01 39.84 NA 5.00 5.00 none
#> 2 lan02 43.81 NA 5.00 5.00 none
#> 3 lan03 47.13 NA 5.00 5.00 favors girls
#> 4 lan04 51.31 NA 5.00 5.00 none
#> 5 lan05 56.09 NA 6.00 6.00 favors girls
#> 6 lan06 60.42 8.32 8.50 8.32 inconsistentx 是題目在圖上的月齡。 有 50% 通過月齡就用它，沒有就用常模中位數到達這個 D 的月齡（age_reach）。 比最小月齡的 P50 還容易的題，below_range 是 TRUE；比最大月齡的 P50 還難的題，above_range 是 TRUE。 這種題沒有通過月齡時，沒有 x。
每題的通過不一樣時：step =。 例如認知在高年級有 4 點題（0–3，至少 1 分算通過）和反向的 5 點題（反向後 0–4，至少 3 分算通過）。 校準不需要「通過」的定義：PCM 讓每題有自己的步驟數，4 點題、5 點題放在同一個校準裡就好。 要「通過」的只有兩個地方：錨題接到 D（anchors 的 step 欄，一題一個），和題目表。 題目表用 step 給每題的通過步驟，沒列的題用 1：
item_tab <- kits_item_table(thr, norm, data = resp, step = c(cogD01 = 3, cogD02 = 3))這時 D 是那一步的門檻，通過月齡算「至少 step 分」，表裡多一欄 step。 不同 step 的題不是同一種里程碑，讀表和圖時要分開看。 反向計分在整理資料時做（例如 5 - x）；同一題在不同波次的點數不同時，要當成兩題（不同的 item_id）。
plot(item_tab)點的顏色是性別 DIF：灰色沒有、橘色有利女孩、藍色有利男孩、黑色兩個世代不一致。
plot(item_tab)plot(item_tab, type = "ladder")每題一條橫線，放在 milestone step 的 D；顏色和上圖相同。 右邊的刻度是常模中位數到達這個 D 的月齡，ages 給要標的月齡。
plot(item_tab, type = "ladder", ages = c(6, 12, 24, 36, 48))kits_norms(..., by = "sex")by 給組的欄，每一組配一個常模。 每一組的常模，和只用這一組的孩子跑 kits_norms() 完全相同。 結果是一張表，組的欄在最前面。
norm_sex <- kits_norms(use, degrees = c(3, 2, 1, 1), by = "sex")
norm_sex |>
select(sex, age_m, P10, P50, P90) |>
head()
#> <kits_norms> scale: D
#> norm fitted on 5 to 49 months (set by kits_norm_ages()); one norm per sex; showing 6 rows
#> sex age_m P10 P50 P90
#> 1 boy 5.00 45.92 54.82 63.16
#> 2 boy 5.50 46.92 55.90 64.33
#> 3 boy 6.00 47.86 56.90 65.43
#> 4 boy 6.50 48.74 57.84 66.46
#> 5 boy 7.00 49.56 58.73 67.43
#> 6 boy 7.50 50.34 59.57 68.35
#> Next: kits_norm_z(scores, x), plot(x), plot(x, scores = scores).plot() 一組畫一格。 好幾個領域疊在一起時（有 domain 欄），用 by = c("domain", "sex")：一列一個領域，一欄一個性別。
plot(norm_sex)同齡 z 用自己那一組的曲線：
use |>
kits_norm_z(norm_sex) |>
summarise(mean_z = mean(z, na.rm = TRUE), .by = sex)
#> sex mean_z
#> 1 boy 8.882258e-04
#> 2 girl -6.938507e-05plot(norm_sex, scores = use) 也把孩子畫在曲線後面，每一格是自己那一組的孩子。
其他報告的圖用法相同：
plot(kits_vertical(s), type = "median")（第 8.4 節）。kits_save_plot(p, "fig.html", "fig.png")：存成互動的 html 和靜態的 png。 html 需要 pandoc（RStudio 和 Quarto 都附了；在 RStudio 外面用時設 RSTUDIO_PANDOC）。 png 需要 webshot2 和一個 Chromium 系的瀏覽器。假設 M36 有兩個新的波次：60、72 個月，各 300 個孩子。 這兩波加了兩題新題 lan13、lan14。 另外 lan10 變容易了（偏移）。
new <- tibble(cohort = "M36", wave_m = rep(c(60, 72), each = 300)) |>
mutate(visit_id = sprintf("M36new_%04d", row_number()),
age_m = wave_m + runif(n(), -1, 1),
sex = sample(c("boy", "girl"), n(), replace = TRUE))
new_betas <- betas[5:14]
new_betas$lan10 <- new_betas$lan10 - 1.2
X_new <- kits_sim(2 * log(new$age_m / 12) + rnorm(600, 0, 0.9), new_betas, prefix = "M36new_")幾波一起進來時，疊成一個矩陣。 資料越多，偏移越容易看出來。
kits_drift()新資料自己校準一次，和 item36 的門檻比。
drift <- kits_drift(X_new, item36)
drift
#> <kits_drift> new data -> ref (scale lan M36): intercept 0.47, slope 1.01
#> items judged: 8 of 8; drifted: lan10
#> suggested for fixing: 7 item(s), in $fix (use kits_item(<new data>, fixed_params = x$fix))
#> item_id max_abs_resid drift suggest_fix reason
#> lan05 0.77 FALSE TRUE
#> lan06 0.93 FALSE TRUE
#> lan07 0.21 FALSE TRUE
#> lan08 0.17 FALSE TRUE
#> lan09 0.30 FALSE TRUE
#> lan10 1.33 TRUE FALSE
#> lan11 0.06 FALSE TRUE
#> lan12 0.11 FALSE TRUE
#> Checks: 1 note (see summary() or kits_check())plot() 畫新舊門檻和那條線；偏移的題是橘色，建議固定的題是藍色。 type = "difference" 畫每個步驟離線多遠。
plot(drift)怎麼讀：
drift = TRUE：這題的門檻離其他題的線太遠，在新資料裡變了。drift 是 NA 的題不判斷，reason 寫原因。suggest_fix：建議固定的題。item36 裡，不判斷。不判斷的題（drift 是 NA）：
min_n 個人。 例如幾乎每個孩子都通過的題，最低的步驟沒有資訊。item36 高（reason 是 “score range differs”）：多了一個步驟。 門檻由這題所有的步驟一起決定，所以舊步驟的門檻也會變，不能直接比。建議固定的題（suggest_fix = TRUE）：
reason 是 “no 0 score or one score only”）：只有 item36 的參數可用。item36 高的題：kits_item() 固定它舊的步驟，只估新的步驟。不建議固定：偏移的題（要重新估）。 最後固定哪些題，由你決定。
限制：
item36 不同。drift$fix 是建議固定的題在 item36 的參數，直接交給 kits_item()。 偏移的題和新題重新估。 新的校準自動在 lan M36 上（固定參數的量尺），不是資料的名字 X_new。
item_new <- kits_item(X_new, fixed_params = drift$fix)
item_new
#> <kits_item> PCM (kits, CML): 600 person-waves, 10 items, 20 steps, 1 item pattern(s)
#> Scale: lan M36
#> Fixed (fixed_params, SE 0): 14 step(s) of 7 item(s)
#> item_id item_score beta beta_se thr thr_se fixed
#> lan05 1 -2.85 0 -3.03 0 yes
#> lan05 2 -1.38 0 -1.20 0 yes
#> lan06 1 -1.86 0 -2.12 0 yes
#> lan06 2 -0.91 0 -0.65 0 yes
#> lan07 1 -1.77 0 -1.94 0 yes
#> lan07 2 -0.21 0 -0.05 0 yes
#> ... 14 more steps (more: print(x, n = 20); all steps: coef(x))
#> Item fit (conditional): infit and outfit within 0.6-1.4 for all 10 items.
#> Next: kits_subj(data, x), kits_linking(x, ...), plot(x, type = "fit"), plot(x), confint(x).
#> Checks: 1 note (see summary() or kits_check())固定的題 SE 是 0；印出來的表頭說固定了幾題，fixed 欄標出是哪幾題。 要看偏移的題和新題：
coef(item_new) |>
filter(item_id %in% c("lan10", "lan13", "lan14"))
#> item_id item_score beta beta_se thr thr_se
#> 1 lan10 1 -0.93 0.27 -1.14 0.22
#> 2 lan10 2 0.32 0.11 0.53 0.09
#> 3 lan13 1 2.02 0.10 1.83 0.08
#> 4 lan13 2 3.39 0.15 3.58 0.13
#> 5 lan14 1 2.81 0.11 2.57 0.09
#> 6 lan14 2 3.88 0.20 4.12 0.17想改名單時，直接改這個表，例如 filter(drift$fix, item_id != "lan06")。 固定參數要用某個世代自己的量尺（這裡是 M36），不能用換算過的參數。
新的校準在 lan M36 上，所以沿用原本的 co 和 d_link。
s_new <- kits_subj(X_new, item_new, persons = new) |>
kits_convert(co) |>
kits_convert(d_link)
s_new |>
summarise(median = median(theta), .by = wave_m)
#> wave_m median
#> 1 60 88.91341
#> 2 72 90.27956原本的常模只到 49 個月，新的波次超出這個範圍。 所以常模要用所有波次重配，步驟和入門篇第 5 節一樣。 兩個表都在 D 上，才疊得起來。
lan_all <- rbind(lan, s_new)
screen_all <- kits_norm_screen(lan_all)
use_all <- filter(lan_all, screen_all$normed)
norm_all <- kits_norms(use_all, degrees = c(3, 2, 1, 1))
range(norm_all$age_m)
#> [1] 5 73
s_new <- s_new |>
kits_norm_z(norm_all)
s_new |>
filter(!is.na(z)) |>
summarise(mean = mean(z), sd = sd(z), .by = wave_m)
#> wave_m mean sd
#> 1 60 0.004149484 0.9654797
#> 2 72 -0.011297320 1.0212105新波次也一起檢查：
kits_check(drift, item_new, s_new, norm_all)
#> object status check value message
#> drift note drifted items 1 not suggested for fixing: lan10; with no drift at all, about 0.5 of the 11 judged steps are flagged by chance (4.6%)
#> item_new note smallest category count 1 fewer than 5 persons in: lan05:0 (2), lan06:0 (1)
#> s_new note share of misfit 8 8.0% of person-waves (above 5%)plot(v, type = "median")kits_vertical() 給每個世代 × 波次一列。 裡面有中位數、P25、P75、人數，和實際月齡的中位數（age_med）。
v_all <- kits_vertical(lan_all)
v_all |>
select(cohort, wave_m, n, age_med, median, P25, P75)
#> <kits_vertical> 7 cohort x wave group(s); scale: D
#> cohort wave_m n age_med median P25 P75
#> M03 6 300 5.96 57.30 52.31 62.15
#> M03 12 300 12.00 64.61 59.72 69.90
#> M03 24 300 24.05 72.94 69.90 81.21
#> M36 36 300 35.96 82.20 77.02 87.79
#> M36 48 300 48.08 84.92 81.54 90.91
#> M36 60 300 59.94 88.91 84.89 93.10
#> M36 72 300 72.02 90.28 86.88 96.05plot(type = "median") 畫中位數，線段是 P25 到 P75，點在實際月齡的中位數。 月齡和波次差超過 1 個月時，標題下面會註明。
plot(v_all, type = "median")好幾個領域時，每個領域各做一次 kits_vertical()，再疊起來： list(lan = v_lan, motor = v_motor) |> list_rbind(names_to = "domain")。 一個領域畫一格。
資料一波一波進來時，用一個題庫（bank）。 目標是：舊的 beta 永遠不變，以前的分數也不用再估；每一波只估新加入的人。 每個步驟第一次估出來就放進題庫，之後不再動。 kits_drift() 的 ref 可以只有 item_id、item_score、beta，門檻由它自己算。
偏移的題在這一波當成新題：改一個題號（這裡加 _w2），自由估計。 其他共同題照 drift$fix 固定，當作錨。 題庫裡舊的 lan10 留給以前的波次，新的 lan10_w2 給這一波以後用。 偏移的題太多時（例如超過能判斷的共同題的 20%，或固定的題少於 5 題），就停下來，由人決定怎麼辦。
bank <- coef(item36)[c("item_id", "item_score", "beta", "beta_se")]
drift_1 <- kits_drift(X_new, bank)
moved <- drift_1$items$item_id[drift_1$items$drift %in% TRUE]
colnames(X_new)[colnames(X_new) %in% moved] <- paste0(moved, "_w2")
item_1 <- kits_item(X_new, fixed_params = drift_1$fix)
new_steps <- anti_join(coef(item_1), bank, by = c("item_id", "item_score"))
bank <- rbind(bank, new_steps[names(bank)])題庫只加進原本沒有的步驟：新題（lan13、lan14）、偏移題的新版本（lan10_w2），和舊題新出現的更高分。 下一波重複這幾行，X_new 換成新的資料，題號先換成每題最新的版本。 舊題多了一個步驟時，舊步驟的 beta 固定不變。 但舊步驟的門檻（RP50）會變，因為門檻由這題所有的步驟一起決定。 所以每一波的分數要在那一波算好、存下來；用後來的題庫重算，數字會稍微不同。 新步驟中間有缺口時（例如有人拿 3 分、沒人拿 2 分），kits_item() 會停下來。 這時先把這一波這題的分數往下併到已固定的最高步，等之後的波次再估新的步驟。 完整的例子在 ?kits_item。
kits_algorithm()每個結果都可以說出它背後的演算法。 kits_algorithm() 印大約 30 行：每一步一個公式、跑這些步驟的函數（|> 接起來）、主要的文獻、去哪裡讀更多。 公式裡的門檻寫成 thr_ik（欄名 thr）；tau 只用在錨題的 GSED 難度。 給結果時，This result 一行是這個結果自己的數字和量尺。
kits_algorithm(co)
#> <kits_algorithm> link: a link: mean/sigma with one round of screening, and the linking error
#> This result: lan M36 -> lan M03: 8 points, 0 dropped (|e| > 2 SD(e)); intercept 2.50, slope 1.05
#> Steps:
#> 1. Points: common steps of two calibrations, or pass step s of the anchors
#> (from_j, to_j) = (thr_ik, thr_ik) or (thr_is, τ_i), τ = the GSED difficulty
#> 2. mean/sigma line on all points (A = slope, B = intercept)
#> A = SD(to) / SD(from), B = mean(to) − A mean(from)
#> 3. One round of screening, then the line again on the kept points
#> e_j = to_j − (B0 + A0 from_j); drop j if |e_j| > cut SD(e)
#> 4. Linking error: R resamples of whole items, the screening redone in each
#> SE(B) = SD_r(B_r), SE(A) = SD_r(A_r), Cov(A, B) over the resamples
#> 5. Linking error at theta (a chain of links adds up the same way)
#> SE_link(θ) = √(Var(B) + θ² Var(A) + 2θ Cov(A, B))
#> 6. Domain check: one slope shared by the domains d, one intercept each
#> A_s = √(Σ (to − mean_d to)² / Σ (from − mean_d from)²), B_d = mean_d to − A_s mean_d from
#> Components:
#> pts <- kits_link_points(from, to) # from, to: two kits_item()
#> pts |> kits_link_screen(cut = 2) |> kits_link_line() |> kits_link_error(n_boot = 200)
#> screen <- algo_screen(pts$from, pts$to, pts$item, cut = 2)
#> screen |> algo_mean_sigma() # the line: B, A
#> screen |> algo_link_error(n_boot = 200) # SE(B), SE(A)
#> algo_line(link, theta) # new = B + A x old
#> References:
#> Marco, G. L. (1977). Journal of Educational Measurement, 14, 139-160.
#> Kolen, M. J., & Brennan, R. L. (2014). Test equating, scaling, and linking (3rd ed.). Springer.
#> Weber, A. M., Rubio-Codina, M., Walker, S. P., et al. (2019). BMJ Global Health, 4, e001724.
#> Read more:
#> tutorial_algorithms/03_linking.html
#> system.file("scripts", "link.R", package = "kits"), parts: points, line, error, domains, slope, rules給主題也可以："cml"、"wle"、"link"、"norms"、"norm_screen"、"dif"、"drift"、"information"。
algo_ 函數Components 裡的 algo_ 函數是演算法本身：給數字、回數字，不記量尺、不檢查輸入。 接起來就是 kits 的結果，例如 kits_stats(X) |> algo_cml() |> algo_thresholds() 是 kits_item(X) 的門檻（入門篇第 1 節）。 人的分數是 algo_jme()。 kits 叫它 WLE；對 PCM，WLE 和 Jeffreys 眾數估計相同，kits 算的是後者（Warm, 1989; Firth, 1993; Magis, 2015）。 它給每個原始分一個 θ 和 SE；答完 8 題的孩子，和 kits_subj() 相同：
cml03 <- kits_stats(X03) |>
algo_cml()
jme03 <- algo_jme(cml03$betas)
full03 <- filter(kits_subj(X03, item03), raw_max == 16)
all.equal(full03$theta, jme03$theta[full03$raw + 1])
#> [1] TRUE校準的結果也包著幾個物件： item03$stats（充分統計量）、 kits_information(item03)（測驗訊息）和 kits_expected_score(item03)（期望原始分），兩個都是 θ 的函數。
item03$stats
#> <kits_stats> sufficient statistics of the CML: 900 person-waves, 8 items, 16 steps, 2 item pattern(s)
#> n_ig, persons with a score of at least g on item i ($steps):
#> item_id item_score n
#> lan01 1 878
#> lan01 2 736
#> lan02 1 859
#> lan02 2 705
#> lan03 1 824
#> lan03 2 573
#> ... 10 more steps
#> Item patterns ($patterns); n_br, persons per raw score r of each item pattern, is in $raw:
#> pattern_id items n max_raw
#> bk1 lan01, lan02, lan03, lan04, lan05, lan06, lan07, lan08 888 16
#> bk2 lan01, lan02, lan03, lan04, lan05, lan06, lan07 12 14
#> Next: kits_cml(x). The algorithm: kits_algorithm(x).
information <- kits_information(item03)
information(c(-2, 0, 2))
#> [1] 2.002848 2.782118 1.933233這兩個函數一題一題拆開，就是 kits_pcm_prob(item03, theta)：每題的類別機率 P(X = x | θ)、期望分數和變異數（題目訊息）。 kits_esf() 正好反過來：給定原始分 r 的條件機率 P(X = x | r)，裡面沒有 θ。
pp <- kits_pcm_prob(item03, seq(-4, 4, by = 1))
pp
#> <kits_pcm_prob> PCM category probabilities: 8 item(s), 9 theta value(s); scale lan M03
#> Items: lan01, lan02, lan03, lan04, lan05, lan06, lan07, lan08
#> theta expected information slope
#> -4.00 0.577 0.534 0.457
#> -3.00 1.396 1.161 0.788
#> -2.00 2.977 2.003 0.801
#> -1.00 5.317 2.605 0.379
#> 0.00 8.043 2.782 -0.018
#> 1.00 10.746 2.552 -0.445
#> 2.00 13.015 1.933 -0.751
#> 3.00 14.560 1.159 -0.735
#> 4.00 15.393 0.554 -0.459
#> P(X = x | theta) of the PCM; information = Var(X | theta), slope = I'(theta).
#> expected and information are kits_expected_score(item)(theta) and
#> kits_information(item)(theta), from the same core.
#> kits_esf() is the other way round: P(X = x | r), conditional on the raw score r,
#> with theta gone; here the probabilities are conditional on theta, not on r.
#> Next: plot(x) (category curves), plot(x, type = "expected"), plot(x, type = "information"),
#> x$prob, x$moments, kits_esf(item).inst/scripts/kits_algorithm() 的最後一行指到一份教學腳本，和它分成的段落（parts）。 腳本是固定的 R 程式：用示範資料一步一步重算，模型用 RTMB 寫，最後和 kits 的結果對答案。 換成自己的資料，改開頭的「inputs」那一段。
list.files(system.file("scripts", package = "kits"))
#> [1] "cml.R" "dif.R" "drift.R" "information.R" "link.R" "norm_screen.R"
#> [7] "norms.R" "wle.R"例如校準的腳本是 system.file("scripts", "cml.R", package = "kits")，複製出來就能跑。
這一節的方法來自 Kolen 和 Brennan（2014）的 Test Equating, Scaling, and Linking（第 3 版）。 括號裡是書上印的頁碼。 連結誤差要用重抽，所以這裡的連結都用 n_boot = 200。
co_b <- kits_linking(item36, item03, n_boot = 200)
d_b <- kits_linking(list(lan = item03, motor = itemm03), anchors, n_boot = 200)n_boot > 0 的連結留下每次重抽的線（co_b$draws）。 在 θ 的連結誤差是 intercept + slope θ 在重抽之間的 SD（pp. 249-252）。 θ = 0 時就是截距的誤差 intercept_error；離 0 越遠，斜率的誤差影響越大。 一串連結（M36 → M03 → D）的兩條線用不同的題，所以重抽是獨立的（pp. 257-258、269-270）。 選定的 θ 也用 kits_convert() 換算：把它們放在一張有 theta 的表裡。 換算後的 theta 在 D 上，se_link 是這個位置的連結誤差。 這張表沒有記下量尺，所以用 from 說它在 lan M36 上。 例如 M36 各波次的中位數：
s36_raw <- kits_subj(X36, item36, persons = persons)
med36 <- s36_raw |>
summarise(theta = median(theta), .by = wave_m)
med36 |>
kits_convert(co_b, d_b, from = "lan M36")
#> wave_m theta se_link
#> 1 36 82.19987 2.360803
#> 2 48 84.91881 2.551101每條連結都有重抽時，kits_convert() 加上 se_link 和 se_total = √(se² + se_link²)。 有一條沒有重抽（例如上面的 co）時，不加欄。 se_link 只是連結題目抽樣的誤差。 同一個世代裡比孩子時，大家用同一條線，不需要它；跨世代比較或報告 D 的不確定性時才需要。
s36_b <- s36_raw |>
kits_convert(co_b) |>
kits_convert(d_b)
s03_b <- kits_subj(X03, item03, persons = persons) |>
kits_convert(d_b)
rbind(s03_b, s36_b) |>
summarise(across(c(se, se_link, se_total), median), .by = wave_m)
#> wave_m se se_link se_total
#> 1 6 4.220297 0.8372636 4.341000
#> 2 12 4.340466 1.2416013 4.449381
#> 3 24 4.866278 1.7143703 5.159431
#> 4 36 4.453088 2.3608035 5.040177
#> 5 48 4.569007 2.5511005 5.232966PCM 沒有鑑別度，所以 PCM 的連結只需要平移（Masters, 1984；p. 229）。 kits 保留自由的斜率，並在診斷裡比較斜率和 1：
summary(co_b)$diagnostics |>
filter(check == "slope vs 1")
#> status check value message
#> ok slope vs 1 0.971 (slope - 1) / slope_error: the slope is within 2 linking errors of 1, as the PCM expects(slope - 1) / slope_error 的絕對值超過 2 時是 note：量尺在兩次校準之間伸縮了，題目的鑑別力不同，值得在報告中說明。 接到錨題（D）的連結不比，因為兩邊的單位不同。
用錨題接到 D 是校準（calibration）或量尺對齊（scale aligning），不是等化（pp. 493-500）。 換算後的 D 能怎麼解釋，要看 KIT 和 GSED 的構念、母群和施測條件有多相似。
kits_vertical()分數都在 D 上之後，每個世代 × 波次算一列（pp. 461-469）。 sd_true 是扣掉測量誤差的 SD。 0 分或滿分的孩子很多時，sd_true 會低估，要當成下限來讀。 effect_size 是書上的 Eq. 9.47；p_sup 是這一波的孩子高過前一波孩子的機率，只看順序。 KIT 的 36 和 48 個月兩個世代都有；這種波次另外比較（attr(v, "cohorts")），分數有 se_link 時也算進連結誤差。 這份模擬資料的兩個世代沒有共同的波次，所以沒有這張表。
v <- kits_vertical(rbind(s03_b, s36_b))
v
#> <kits_vertical> 5 cohort x wave group(s); scale: D
#> cohort wave_m n age_med mean sd median P25 P75 sd_true growth effect_size p_sup se_link
#> M03 6 300 5.96 56.52 7.42 57.30 52.31 62.15 5.96 NA NA NA 0.80
#> M03 12 300 12.00 65.49 7.71 64.61 59.72 69.90 6.17 8.97 1.19 0.80 1.29
#> M03 24 300 24.05 74.93 7.67 72.94 69.90 81.21 4.68 9.44 1.23 0.80 1.83
#> M36 36 300 35.96 81.76 7.49 82.20 77.02 87.79 5.83 NA NA NA 2.33
#> M36 48 300 48.08 85.93 8.04 84.92 81.54 90.91 6.18 4.17 0.54 0.65 2.62
#> Checks: OK最大的波次的 sd_true 比前面小很多時，診斷會註記量尺縮小（scale shrinkage，p. 468）。
plot(v) 畫平均 ± SD。
plot(v)plot(v, type = "median") 畫中位數，線段是 P25 到 P75。
plot(v, type = "median")kits_precision()每個題本在量尺的每個位置量得多準（pp. 405、413、463-465）。 原始分的分佈用 Lord-Wingersky 遞迴（p. 232），原始分換成和 kits_subj_table() 相同的 WLE。 rmse 結合偏誤和 SD；可用範圍是 RMSE 小於 criterion 的區間。 這裡兩個世代的題本都換到 D，criterion 用 D 的單位（5 D）。
pr03 <- kits_precision(item03, list(M03 = colnames(X03)), link = d_b, criterion = 5)
pr36 <- kits_precision(item36, list(M36 = colnames(X36)), link = list(co_b, d_b), criterion = 5)
rbind(attr(pr03, "ranges"), attr(pr36, "ranges"))
#> form n_items max_raw from to min_rmse at
#> 1 M03 8 16 46.05 72.94 4.13 59.15
#> 2 M36 8 16 64.75 89.44 4.34 76.37給 kits_subj_table() 的結果時，每個波次的題本各一條曲線。
plot(pr03)kits_subj_table() 和 kits_expected_score()同一個校準裡的兩組題，一組的原始分在另一組大約是幾分？ 例如少答 lan08 的孩子（7 題），在完整的 8 題上的期望原始分。 先從 kits_subj_table() 拿這個作答組合每個原始分的 θ（WLE），再代進完整題本的期望分數：
tab7 <- kits_subj_table(item03) |>
filter(pattern_id == "bk2") # 少答 lan08 的作答組合（7 題）
mutate(tab7, full = kits_expected_score(item03)(theta)) # 完整 8 題的期望原始分
#> <kits_subj_table> every raw score of every item pattern, from the item parameters; scale: lan M03
#> pattern_id raw theta se full
#> 1 bk2 0 -4.27 1.55 0.45
#> 2 bk2 1 -3.04 0.94 1.35
#> 3 bk2 2 -2.38 0.78 2.27
#> 4 bk2 3 -1.89 0.70 3.21
#> 5 bk2 4 -1.46 0.65 4.17
#> 6 bk2 5 -1.07 0.63 5.14
#> 7 bk2 6 -0.69 0.62 6.14
#> 8 bk2 7 -0.32 0.62 7.15
#> 9 bk2 8 0.05 0.62 8.17
#> 10 bk2 9 0.42 0.63 9.22
#> 11 bk2 10 0.82 0.66 10.28
#> 12 bk2 11 1.26 0.70 11.38
#> 13 bk2 12 1.77 0.79 12.55
#> 14 bk2 13 2.45 0.96 13.81
#> 15 bk2 14 3.72 1.57 15.22PCM 的原始分是充分統計量，所以原始分 → θ 的表就是計分；兩組題不需要另外做原始分對原始分的等化（Kolen & Brennan, 2014, §6.4）。 固定的一組題 S 的分量表分數是 kits_expected_score(item, items = S)(theta)，SE 是 IS(θ) × SE(θ)，IS 來自 kits_information(item, items = S)。
kits_link_check()世代連結只靠共同題的門檻。 書上建議再用別的方法檢查同一條線（pp. 190、238、295-296）。 kits_link_check() 用兩個不靠篩過的共同題的方法：
這份模擬資料的兩個世代沒有共同的波次。 KIT 的 36 和 48 個月兩個世代都有。 所以這裡多模擬 M03 的一個 36 個月的波次。 這一波的題本是 lan03 到 lan10，所以 M03 的校準多了 lan09、lan10。
p36 <- tibble(cohort = "M03", wave_m = 36, visit_id = sprintf("M03x_%04d", 1:300)) |>
mutate(age_m = 36 + runif(n(), -1, 1))
X03_36 <- kits_sim(2 * log(p36$age_m / 12) + rnorm(300, 0, 0.9), betas[3:10], prefix = "M03x_")
items_x <- sprintf("lan%02d", 1:10)
X03x <- matrix(NA, nrow(X03) + 300, 10, dimnames = list(c(rownames(X03), p36$visit_id), items_x))
X03x[rownames(X03), colnames(X03)] <- X03
X03x[p36$visit_id, colnames(X03_36)] <- X03_36
item03x <- kits_item(X03x, scale = "lan M03x")
co_x <- kits_linking(item36, item03x, n_boot = 200)
persons_x <- bind_rows(select(persons, visit_id, wave_m), select(p36, visit_id, wave_m))
chk <- kits_link_check(co_x, item36, item03x, X36, X03x, persons = persons_x)
chk
#> <kits_link_check> lan M36 -> lan M03x: the link and independent lines, new value = intercept + slope x old value
#> method intercept slope
#> common items (the link) 1.81 1.02
#> same-age groups, wave 36 1.57 1.27
#> concurrent CML 1.81 1.01
#> common items (the link): 12 common steps of 6 items, 0 dropped (|e| > 2 SD(e))
#> same-age groups, wave 36: 300 person-waves of from, 300 of to (normal ability distributions matched)
#> concurrent CML: 12 items, 1800 person-waves calibrated together
#> Shared waves (each cohort on its own calibration):
#> wave_m n_from mean_from sd_from se_mean_from se_sd_from n_to mean_to sd_to se_mean_to se_sd_to converged
#> 36 300 0.64 0.82 0.06 0.06 300 2.38 1.04 0.07 0.07 TRUE
#> At theta (scale lan M36; value on lan M03x; diff = method minus the link):
#> theta method value diff se_method se_link z z_both
#> 0.70 common items (the link) 2.53 0.00 NA 0.07 0.00 0.00
#> 1.08 common items (the link) 2.91 0.00 NA 0.08 0.00 0.00
#> 0.70 same-age groups, wave 36 2.46 -0.07 0.11 0.07 -1.02 -0.56
#> 1.08 same-age groups, wave 36 2.93 0.02 0.12 0.08 0.27 0.15
#> 0.70 concurrent CML 2.52 -0.01 NA 0.07 -0.09 -0.09
#> 1.08 concurrent CML 2.90 -0.01 NA 0.08 -0.11 -0.11
#> Checks: OK怎麼讀：
Shared waves 是兩個世代在 36 個月的能力分佈，各自在自己的量尺上（M36 和 M03x）。At theta 是在 M36 各波次的中位數，三條線的換算值和差（diff）。z 是差除以連結的連結誤差；z_both 再加上同齡兩群的抽樣誤差（se_method）。同齡的兩群假設兩個世代同年齡的孩子來自同一個母群。 兩個世代出生的年份不同，測過的次數也不同（世代效果）。 所以同齡的線和連結不同時，不一定是連結錯了，也可能是世代真的不同。 合併校準用所有的共同題，包括篩點丟掉的題，所以它假設共同題沒有偏移。
spread、unit預設的規則：殘差超過 2 個 SD 的點丟掉，一個步驟一個點。 spread 換篩點的尺度（pp. 190、202-203、288-289）：
"sd"（預設）：|e| > cut × SD(e)。"mad"：|e| > cut × 1.4826 × MAD(e)，離群的點不會把邊界撐寬。"se"：殘差除以兩邊門檻的 SE；錨題（GSED 的 tau）沒有 SE，只用 KIT 這邊的。 這只適合小樣本：KIT 的 thr_se 很小，篩掉的點會太多。unit = "item" 時，一題有一個步驟被標記，整題都丟掉（書上篩的是整題）。
rules <- expand_grid(unit = c("step", "item"), spread = c("sd", "mad", "se")) |>
mutate(link = map2(spread, unit, \(sp, u) kits_linking(item36, item03, spread = sp, unit = u)))
rules |>
transmute(spread, unit,
intercept = map_dbl(link, "intercept"),
slope = map_dbl(link, "slope"),
dropped = map_chr(link, \(l) paste(l$dropped, collapse = " ")))
#> # A tibble: 6 × 5
#> spread unit intercept slope dropped
#> <chr> <chr> <dbl> <dbl> <chr>
#> 1 sd step 2.50 1.05 ""
#> 2 mad step 2.50 1.05 ""
#> 3 se step 2.50 1.05 ""
#> 4 sd item 2.50 1.05 ""
#> 5 mad item 2.50 1.05 ""
#> 6 se item 2.50 1.05 ""連結記下用的規則（settings$spread、settings$unit），kits_algorithm(link) 的 This result 一行寫出這個規則。
篩點把換算值移動多少？ summary(link, theta = ) 比較篩點前（所有點）和篩點後的線。 比較的表在 $screening，這裡用 M03 各波次的中位數：
s03_raw <- kits_subj(X03, item03, persons = persons)
med03 <- s03_raw |>
summarise(median = median(theta), .by = wave_m)
summary(d_b, theta = med03$median)$screening
#> theta before after shift se_link z
#> 1 -0.37 57.30 57.30 0 0.84 0
#> 2 0.69 64.61 64.61 0 1.24 0
#> 3 1.90 72.94 72.94 0 1.71 0shift 是篩點後減篩點前（D 的單位），z 是移動量除以這個位置的連結誤差。 這裡沒有篩掉錨題，所以 shift 是 0。
有領域的連結（D 連結），印出來多一行：留下的點在各領域的數目。
d_b
#> <kits_link> new value = intercept + slope x old value
#> scales: lan M03, motor M03 -> D
#> intercept = 59.84 (linking error 0.98)
#> slope = 6.90 (linking error 0.40)
#> points from: the pass step of each anchor item (score >= 1)
#> points: 8 (items: 8); dropped: 0
#> before screening: intercept 59.84, slope 6.90; screening: |e| > 2 SD(e)
#> kept points by domain: lan 5 (62%), motor 3 (38%)
#> domain check: one slope shared by the domains, fitted with domain intercepts (kept points);
#> it is not the slope of the line above:
#> domain n intercept slope
#> lan 5 60.66 7.02
#> motor 3 58.80 7.02
#> Next: kits_convert(scores, x), plot(x), summary(x, theta = ), confint(x).
#> Checks: OK一個領域超過 2/3 時有註記：共同題（錨題）應該是整份測驗的縮影（p. 239）。 KIT 的 D 連結留下的錨題大多是動作，會有這個註記。
兩個校準之間的連結，也看共同題佔各校準題數的比例（p. 288）：
summary(co_b)$diagnostics |>
filter(check == "common items / items")
#> status check value message
#> ok common items / items 50 4 kept common items: 50% of the 8 items of from, 50% of the 8 items of to低於 20% 時有註記。 書上這個經驗法則是給 40 題以上的測驗。
coef(item, rp = )、kits_item_table(rp = )門檻 thr 是 RP50：能力在這裡時，得到 k 分以上的機率是 .5。 題目地圖（item mapping）常用更高的 RP，例如 .8（pp. 414-419）。 coef(item, rp = 0.8) 給 RP80 的門檻：
coef(item03) |>
select(item_id, item_score, thr) |>
mutate(thr80 = coef(item03, rp = 0.8)$thr) |>
head()
#> item_id item_score thr thr80
#> 1 lan01 1 -2.90 -1.80
#> 2 lan01 2 -1.11 0.14
#> 3 lan02 1 -2.32 -1.28
#> 4 lan02 2 -0.83 0.38
#> 5 lan03 1 -1.84 -0.72
#> 6 lan03 2 0.06 1.33RP 越高，題目放在越高的位置。 接到 D 一直用 RP50，因為 GSED 的 tau 就是 RP50（一半的孩子通過的 D）。 RP80 的門檻換到 D，用同一條線：kits_convert(coef(item03, rp = 0.8), d_link)。
月齡也一樣。 kits_item_table(rp = 0.8) 用 RP80 的 D 放題目，月齡也換成 80% 的孩子通過的月齡（age80）：
thr80 <- coef(item03, rp = 0.8) |>
kits_convert(d_link)
item_tab80 <- kits_item_table(thr80, norm, data = X_all, rp = 0.8)
item_tab80 |>
select(item_id, D, age80, x) |>
head()
#> item_id D age80 x
#> 1 lan01 47.41 NA 5.00
#> 2 lan02 50.99 NA 5.00
#> 3 lan03 54.84 6.13 6.13
#> 4 lan04 58.73 8.61 8.61
#> 5 lan05 63.55 13.05 13.05
#> 6 lan06 67.77 17.92 17.92thr 的 RP 要和 rp 一樣，不一樣時函數會停下來。 每個波次都已經超過 80% 的孩子通過時，age80 是 NA。 兩種圖（plot(item_tab80)、type = "ladder"）會標出 RP。
每次新的波次進來，照這張表檢查（Kolen & Brennan, 2014，Table 8.8 p. 336、p. 333）：
kits_drift()（pp. 288-289）。drift 對第一次的校準做，不對上一次的更新做，誤差才不會累積（pp. 293-294）。kits_link_check()、kits_convert() 的 se_link（pp. 295-296）。kits_item() 的診斷 smallest main item pattern。common items / items、domain balance of kept points。kits_vertical() 的世代比較（z、z_link、sd_ratio）。kits_info()。?kits_check 也有這張表。 新的檢查都在 kits_check() 裡：
kits_check(item03x, co_x, chk, d_b)
#> object status check value message
#> item03x note smallest main item pattern 300 2 main item pattern(s) (5% or more of the person-waves, or the largest): 888, 300 person-waves; 1 other item pattern(s) hold 1%; 1 below 400, the rule of thumb for the Rasch model (Kolen & Brennan p. 304)citation("kits") 是 kits 本身。 用到的方法，原始文獻列在 ?"kits-package" 的 References，例如 PCM（Masters, 1982）、WLE（Warm, 1989）、 mean/sigma（Marco, 1977）、D score（Weber et al., 2019）、shash（Jones & Pewsey, 2009）、 cNORM（Lenhard & Lenhard, 2021）、DIF（Bechger & Maris, 2015）、 連結誤差、垂直量尺和原始分對照（Kolen & Brennan, 2014）。