Bayesian Plausible-Value Analysis with One Calibrated Stacked Fit • pvstackr

pvstackr fits survey-weighted linear regression models to large-scale assessment data such as PISA, with plausible values as the outcome, and reports estimates, standard errors and intervals for the fixed effects. The standard analysis of such data fits the model once for each plausible value, with a replicate-weight estimate of each fit’s sampling variance, and combines the results with Rubin’s rules for multiple imputation. An analysis of PV1 alone has standard errors that are too small, because it leaves out the between-imputation variance: the variation of the estimates from one plausible value to the next.

Installation

Install pvstackr from GitHub:

# install.packages("pak")
pak::pak("joonho112/pvstackr")

pvstackr installs and loads without brms or Stan, and its examples and default tests run without them. The engine that pvstackr bundles for stack_direct, selected with pv_control(backend = "brms"), needs the brms and posterior packages. It samples through cmdstanr if that package is installed, and CmdStan must then be configured or the fit stops with an error; without cmdstanr it samples through rstan. Getting started with pvstackr shows the call for a first fit and how to read its output.

What pvstackr computes

pv_brr_target() computes the design-based target with the standard analysis. For each plausible value it fits the regression by weighted least squares with the final survey weight, and estimates the sampling covariance of the coefficients from the replicate weights by balanced repeated replication with Fay’s method (BRR–Fay, as in PISA); Rubin’s rules then combine the plausible values. pvstackr calls the result, the combined estimates with their covariance matrix and degrees of freedom, the target.

The default method of pv_fit(), stack_direct, fits one Bayesian model, weighted by the survey weights, to the stacked data: one copy of the data for each plausible value, with that plausible value as the outcome. The Cholesky calibration correction (CCC) then transforms the fixed-effect draws of this fit to have the mean and covariance of the target. stack_direct requires center = "target" (the default in pv_control()), under which the reported estimates, standard errors and degrees of freedom are those of the target; the Bayesian fit does not change them. The stacked fit supplies the calibrated draws, whose mean and covariance are those of the target, and diagnostics that compare the fit with the target; quantiles of the draws do not carry the target’s degrees of freedom. The companion methods preprint (Lee, Williams and Savitsky 2026, https://doi.org/10.5281/zenodo.22407935) describes the stacked fit and its calibration.

The other two methods, per_pv and stack_psis, have no bundled engine and do not use the target:

Method	What it does	What you supply
`per_pv`	combines one Bayesian fit per plausible value	the posterior draws of each fit, or fitting functions that pvstackr calls once per plausible value (they do not receive the survey weights automatically)
`stack_psis`	reweights the draws of one stacked fit toward each plausible value	the stacked draws, or a fitting function that pvstackr calls once on the stacked data; importance weights and a Pareto k-hat value for each plausible value, computed outside pvstackr (for example by Pareto smoothed importance sampling with the loo package); the name and version of the program that computed them

Both combine the plausible values with Rubin’s rules, but their variance within each plausible value is model-based, taken from the posterior draws, so their intervals are always descriptive. By pvstackr’s reporting rule, only a stack_direct fit whose target uses Barnard–Rubin degrees of freedom (df_method = "barnard_rubin", not the default) has intervals that can be read as confidence intervals with nominal coverage; Reading and reporting the results explains the interval columns and Comparing the three fitting methods compares the methods.

Scope

Only the fixed effects (the intercept and the slopes) are reported; other parameters, such as the residual standard deviation sigma, are neither calibrated nor reported. pv_brr_target(), stack_direct and stack_psis stop with an error on a random-effect term such as (1 | school).

An example on synthetic data

pvstackr includes pisa_tiny.csv, a small synthetic data set with PISA column names (12 students, two plausible values in reading, a final weight and four replicate weights), and a stack_direct fit of these data that was made without a sampler (Getting started with pvstackr explains how).

The data were made up for the examples and the package tests. They contain no real PISA records and are not suitable for real inference.

pv_design() finds and checks the plausible-value and weight columns; the two functions it uses, detect_pisa_pv_columns() and detect_pisa_brr_replicate_weights(), can also be called on their own.

library(pvstackr)

pisa_tiny <- read.csv(
  system.file("extdata", "pisa_tiny.csv", package = "pvstackr")
)

design <- pv_design(
  data = pisa_tiny,
  formula = OUTCOME ~ x + female,  # OUTCOME stands for each plausible value
  pv_suffix = "READ",              # finds PV1READ and PV2READ
  expected_M = 2L,                 # stop unless 2 plausible values are found
  expected_R = 4L,                 # stop unless 4 replicate weights are found
  id_cols = "CNTSTUID"
)

design
#> pvstackr design
#>   rows: 12
#>   formula: OUTCOME ~ x + female
#>   plausible values: 2
#>   final weight: W_FSTUWT
#>   replicate weights: 4
#>   fay_k: 0.5
#>   design hash: fa78c04b

pv_brr_target() computes the target from the same data and formula:

target <- pv_brr_target(
  data = pisa_tiny,
  formula = OUTCOME ~ x + female,
  pv_cols = design$pv_cols,
  weight_col = design$weight_col,
  rep_weight_cols = design$rep_weight_cols,
  fay_k = design$fay_k,
  id_cols = design$id_cols
)

target
#> pvstackr BRR-Fay target
#>   fixed effects: 3
#>   plausible values: 2
#>   replicate weights: 4
#>   fay_k: 0.5
#>   df method: classic
#>   interval role: descriptive_classic_rubin
#>   source: external_brr_fay_rubin

Fitting the model needs a Bayesian sampler, so the example reads the fit that comes with pvstackr:

fit <- readRDS(
  system.file("extdata", "examples", "pisa_tiny_stack_direct.rds",
              package = "pvstackr")
)$fit

fit
#> pvstackr fit
#>   method: stack_direct
#>   status: ok
#>   fixed effects: 3
#>   target: external_brr_fay_rubin
#>   draws: not retained
#>   diagnostics: preflight, sampler, sampler_gate, stack_fit, stack_fit_warnings, ccc
#>   interval note: intervals are descriptive rather than coverage-claimable.

get_estimates(fit)[, c("term", "estimate", "se", "df", "conf_low", "conf_high")]
#>          term   estimate        se       df  conf_low conf_high
#> 1 b_Intercept 457.894088 1.2873118 1.021194 442.31804 473.47013
#> 2         b_x  46.883361 0.3717929 1.402308  44.41457  49.35215
#> 3    b_female   2.143702 3.5550309 1.013730 -41.60687  45.89428

The intervals are wide because the degrees of freedom are between 1.0 and 1.4 in this example, and descriptive because the target uses the default, classic degrees of freedom; Getting started with pvstackr explains both. get_target(), get_draws() and get_diagnostics() return the target, the calibrated draws and the diagnostics of a fit (Reading and reporting the results describes them); this fit was saved without its draws, so get_draws(fit) returns NULL.

Fitting real data

On real data the steps are the same, except that the model is fitted with pv_fit() and a Bayesian engine instead of read from a file. The code below is not run here: pvstackr includes no real data, and the fit needs a sampler. pisa_country stands for the student data of one country in PISA 2022, with ten plausible values in mathematics and 80 replicate weights, and the fit uses the engine bundled with pvstackr:

design <- pv_design(
  data = pisa_country,
  formula = OUTCOME ~ ESCS + ST004D01T,  # ST004D01T is the student's gender
  pv_suffix = "MATH",                    # finds PV1MATH to PV10MATH
  expected_M = 10L,
  expected_R = 80L,
  id_cols = "CNTSTUID"
)

target <- pv_brr_target(
  data = design$data,
  formula = design$formula,
  pv_cols = design$pv_cols,
  weight_col = design$weight_col,
  rep_weight_cols = design$rep_weight_cols,
  fay_k = design$fay_k,
  id_cols = design$id_cols
)

fit <- pv_fit(
  data = design$data,
  formula = design$formula,
  target = target,
  method = "stack_direct",
  control = pv_control(
    method = "stack_direct",
    backend = "brms",  # the engine bundled with pvstackr
    seed = 20260607
  )
)

summary(fit)

The bundled engine saves the brms fit in a folder cache in the working directory and reuses it for later calls with the same model and data; read cache_dir in ?pv_fit_direct before you change the sampler settings.

Instead of the bundled engine you can pass to pv_fit() three functions of your own, fit_function, draws_function and diagnose_function, which The full analysis workflow describes.

PISA 2022 names the plausible values by subject, so for columns such as PV1MATH set pv_suffix = "MATH". The default, pv_suffix = "", matches only bare plausible-value names such as PV1.

Two counts give the size of an analysis. The stacked data have N * M rows for N students and M plausible values, and pv_brr_target() makes M * (R + 1) weighted least-squares fits for R replicate weights, 810 for PISA 2022, before any Bayesian sampling. Real PISA data: access, memory and reproducibility, also available as vignette("a5-real-pisa-guidance", package = "pvstackr"), covers obtaining the files under their terms, where to keep the data and the engine’s cache files, memory and running time, and what to record so that the analysis can be repeated.

Articles

The package website has ten articles. The Applied articles show how to run an analysis and read its result; the Method articles give the formulas, their sources and the assumptions under which they hold.

Applied article	What it covers
Getting started with pvstackr	Installing pvstackr, the `pv_fit()` call for `stack_direct`, and reading the example fit and its intervals
The full analysis workflow	The steps of an analysis in order, from checking the columns with `pv_design()` to reading the result; why `stack_direct` requires `center = "target"`; fitting functions of your own
Reading and reporting the results	The functions that read a fit, the estimate table and its interval columns, the fraction of missing information, and what to report
Comparing the three fitting methods	How `stack_direct`, `per_pv` and `stack_psis` differ, a comparison of three example fits with `pv_compare_methods()`, two cautions in reading it, and pvstackr’s recommendation for choosing a method
Real PISA data: access, memory and reproducibility	Obtaining PISA files under their terms, declaring the PISA 2022 design, memory and running time, a reproducibility checklist, and what the companion preprint shows about PISA

Method article	What it covers
Plausible values, survey weights and notation	Plausible values, final and replicate weights with Fay’s coefficient, the survey-weighted model that pvstackr fits, and the symbols and main formulas of the Method articles
The design-based target: BRR–Fay replicate weights and Rubin’s rules	The BRR–Fay replicate covariance, Rubin’s rules, the classic and Barnard–Rubin degrees of freedom and the fraction of missing information, with their sources and a check on the example target
One stacked fit and the Rubin point estimate	The weighted stacked data, the stacked fractional posterior, and the stacked fixed-effect point identity of the companion preprint: when the stacked estimate equals the Rubin point estimate, and what the identity does not cover
Calibrating the fixed-effect draws (CCC)	The Cholesky calibration correction, its exact moment properties, the algorithm that pvstackr uses, and the diagnostics `delta_c_max`, `delta_c_rel` and `kappa_A` with pvstackr’s thresholds
Choosing a method and reading its intervals	The formulas of the three methods, the importance weights and Pareto k-hat values of `stack_psis`, and why, by pvstackr’s reporting rule, only `stack_direct` fits whose target uses Barnard–Rubin degrees of freedom have intervals that can be read as confidence intervals

Citation

citation("pvstackr") gives the reference for the package:

citation("pvstackr")

Author and maintainer: JoonHo Lee (jlee296@ua.edu, ORCID 0009-0006-4019-8703).

For the method, cite the companion methods preprint, which describes the stacked fit and its calibration:

Lee, J., Williams, M. R., and Savitsky, T. D. (2026). One Markov Chain Monte Carlo Fit for Many Plausible Values: A Calibrated Stacked Posterior Workflow for Bayesian Multilevel Models of Large-Scale Assessment Data. Zenodo preprint, version 1. https://doi.org/10.5281/zenodo.22407935

The preprint’s default target comes from model-based fits of each plausible value, and its PISA application has a random intercept for schools; pvstackr 0.2.x computes only a target from replicate weights, for models without random effects.

pvstackr computes the standard analysis itself, without the survey package: pv_brr_target() fits the weighted least-squares regressions with lm.wfit() from the stats package. The optional development tests compare the results with those of svyglm() from the survey package on a Fay replicate-weight design. pvstackr combines the plausible values with Rubin’s rules, as the multiple-imputation packages mitools and mice combine the results of imputed data sets. For other analyses of large-scale assessment data with plausible values and replicate weights, see intsvy and BIFIEsurvey.

Getting help

Report bugs and ask questions at https://github.com/joonho112/pvstackr/issues.

Status

The lifecycle stage of pvstackr is experimental: its functions and arguments may change in later versions. In version 0.2.1 the documentation was rewritten; the computations are those of 0.2.0. The package includes no real PISA data. pvstackr is an independent research package and is not affiliated with or endorsed by the OECD or the PISA programme.

pvstackr