Plausible values, survey weights and notation
JoonHo Lee
2026-09-11
Source:vignettes/m1-foundations-and-notation.Rmd
m1-foundations-and-notation.RmdAbstract
This article defines the quantities that pvstackr works with: plausible values, the final and replicate survey weights with Fay’s coefficient, the survey-weighted linear model whose fixed effects pvstackr reports, and the weighted likelihood of each plausible value. It lists the symbols and the main formulas of the Method articles and reads the design of the bundled synthetic data.
Plausible values and the replicate-weight design
Plausible values
Large-scale assessments such as PISA do not observe proficiency
directly. Each student receives \(M\)
plausible values: draws from the student’s posterior distribution of
proficiency given the item responses and background variables (Mislevy 1991; von Davier et
al. 2009). They are multiple imputations of a latent
variable, so an analysis fits its model once for each plausible value
and combines the \(M\) results with
Rubin’s rules (Rubin 1987). This
convention presumes that the analysis model is congenial to the model
that generated the plausible values, as the companion preprint notes
(Lee et al. 2026, sec. 2); pvstackr does
not check it. The variation of the \(M\) estimates, the between-imputation
variance, reflects the uncertainty that remains about proficiency; an
analysis that uses only PV1 omits it, so its standard
errors are too small. In pvstackr the plausible values are \(M\) columns of the data, such as
PV1READ and PV2READ, and the model formula
names the outcome OUTCOME, which pvstackr replaces by each
of these columns in turn.
Final and replicate weights
PISA selects schools and then students within schools, and the selection probabilities vary. The final weight \(W_i\) of student \(i\) combines the inverse selection probabilities of the school and of the student with adjustment factors for non-response and trimming (OECD 2024). The data also carry \(R\) replicate weights \(W_i^{(r)}\), \(r = 1, \dots, R\). An estimate is computed with the final weight and again with each replicate weight, and the spread of the \(R\) replicate estimates around the full-sample estimate estimates its sampling variance under the design.
Fay’s method
PISA builds its replicate weights by balanced repeated replication
(BRR) with Fay’s method (OECD 2024). In
BRR the sample of each variance stratum is split into two halves, and a
replicate multiplies the weights of one half by 2 and those of the other
by 0. Fay’s method multiplies them by \(2 -
k\) and \(k\) instead, where
\(0 \le k < 1\) is the Fay
coefficient, so each replicate changes the weights by \(1 - k\) times as much as in BRR. For linear
statistics the squared deviations of the replicate estimates from the
full-sample estimate are then smaller by the factor \((1 - k)^2\), and Fay’s variance estimate,
the mean of these squared deviations over the \(R\) replicates, is multiplied by \(1/(1 - k)^2\) (Judkins 1990). pvstackr applies this estimate
to the vector of regression coefficients, with the factor \(1/\{R(1 - k)^2\}\) in front of the sum of
the outer products of the replicate deviations. The target stores this
factor as fay_variance_multiplier, and The design-based
target: BRR–Fay replicate weights and Rubin’s rules gives the
formula.
In PISA 2022 the sampled schools are paired into variance strata (with a triple, which gets other factors, where a stratum has an odd number of schools), and each of the 80 replicates multiplies the weights of one school in each pair by 1.5 and those of the other by 0.5, so \(k = 0.5\). The non-response adjustments and the trimming are then repeated for each replicate, so the replicate weights differ slightly from these factors, and they are columns of the data file. PISA 2022 also provides ten plausible values per domain, and the factor \(1/\{R(1 - k)^2\} = 1/(80 \times 0.25) = 0.05\) is the one in its variance formula (OECD 2024). The bundled synthetic example data are much smaller:
| Data | \(M\) | \(R\) | \(k\) | \(1/\{R(1-k)^2\}\) |
|---|---|---|---|---|
| PISA 2022 | 10 | 80 | 0.5 | 0.05 |
Bundled example data, pisa_tiny
|
2 | 4 | 0.5 | 1 |
In the example the factor is exactly 1 because \(R(1 - k)^2 = 4 \times 0.25 = 1\), a property of these two numbers only.
pv_brr_target() takes the final weight and at least two
replicate weights as columns of the data (weight_col,
rep_weight_cols) and the Fay coefficient as
fay_k (default 0.5). It does not check fay_k
against the replicate weights, and it requires every weight to be
positive, so plain BRR weights, which are zero for half of the sample,
are not accepted.
The model pvstackr fits
pvstackr works with one model: a linear regression of the plausible value on the covariates, weighted by the survey weights. For plausible value \(m\) and student \(i = 1, \dots, N\),
\[ y_i^{(m)} = x_i^\top \beta + \varepsilon_i^{(m)}, \qquad \varepsilon_i^{(m)} \sim N(0, \sigma^2), \]
where \(x_i\) is the row of the
\(N \times p\) design matrix \(X\) for student \(i\) (the intercept and the covariates of
the formula), \(\beta\) holds the \(p\) fixed-effect coefficients, and \(\sigma\) is the residual standard
deviation. The covariates, and with them \(X\), are the same for every plausible
value; only the outcome changes. The formula of the examples,
OUTCOME ~ x + female, has \(p =
3\) coefficients, which pvstackr names b_Intercept,
b_x and b_female.
pvstackr uses this model in two places, both with the survey weights.
pv_brr_target() computes, for each plausible value, the
weighted least-squares estimate of \(\beta\) with the final weights
(lm.wfit() from the stats package) and its covariance from
the replicate weights; neither step uses the normal distribution of the
errors. The stacked fit of stack_direct fits the model,
normal errors included, once to all \(M\) plausible values together, with the
survey weights in its likelihood (below). It
uses a Bayesian engine, either the brms engine bundled with pvstackr or
fitting functions that you supply, and it accepts only the Gaussian
family with the identity link; any other family stops it with an error.
Only the fixed effects are reported: the stacked fit also estimates
\(\sigma\), which pvstackr neither
calibrates nor reports.
The model has no random effects. A formula with a random-effect term
such as (1 | CNTSCHID) stops pv_brr_target(),
stack_direct and stack_psis with an error, so
pvstackr 0.2.x fits no school-level variance. per_pv fits
no model itself: it pools, with Rubin’s rules, draws for each plausible
value that you supply or that your fitting function returns. It does not
check the formula for random-effect terms and does not add the survey
weights to it, so for per_pv the model is the one that your
function fits.
School means as covariates
pvstackr computes no school or cluster means. A school mean that you compute and store as a column, for example the weighted mean \(\bar x_j\) of a covariate \(x\) over the sampled students of school \(j\), is an ordinary covariate. When pvstackr stacks the data, it copies every covariate unchanged into each of the \(M\) copies, so such a column has the same value for every plausible value. The companion preprint (Lee et al. 2026) needs this invariance. Its stacked fixed-effect point identity (Theorem 4.1 of the preprint) requires, as condition (R2), a design matrix common to all plausible values, which holds because the covariates, the analyzed rows and the cluster membership do not change with the plausible value (Appendix A.2 of the preprint). One stacked fit and the Rubin point estimate lists the conditions.
With a school mean among the covariates, the regression can separate the association within schools from the association between schools. For student \(i\) in school \(j\), the within–between parameterisation enters the deviation \(x_{ij} - \bar x_j\) with slope \(\beta^W\) and the school mean \(\bar x_j\) with slope \(\beta^B\). It is a reparameterization of the model that enters \(x_{ij}\) and \(\bar x_j\) together, the form of Mundlak (1978), who added the group means of the covariates to an error-components model:
\[ \beta^W (x_{ij} - \bar x_j) + \beta^B \bar x_j = \beta^W x_{ij} + (\beta^B - \beta^W)\, \bar x_j . \]
In pvstackr both forms are covariates of the single-level model above; neither adds a school random effect.
The model of the companion preprint
The preprint develops the stacked fit and its calibration for a two-level model with a school random intercept (its Eq. (2)). In the notation of this article,
\[ y_{ij}^{(m)} = \beta_0 + \beta^W (x_{ij} - \bar x_j) + \beta^B \bar x_j + \zeta_j + \varepsilon_{ij}^{(m)}, \qquad \zeta_j \sim N(0, \tau^2), \quad \varepsilon_{ij}^{(m)} \sim N(0, \sigma^2), \]
where \(\bar x_j\) is the school mean of the covariate, which the preprint’s application estimates by the weighted mean of the sampled students; the preprint writes \(\theta_{ij}^{(m)}\) for the plausible value. Its weighted likelihood multiplies each student’s log-likelihood contribution by the student’s normalized weight inside the integral over the random intercept \(\zeta_j\) (its Eq. (3)), and the variance components \(\tau^2\) and \(\sigma^2\) are nuisance parameters. The preprint’s default target is computed from model-based fits of each plausible value, by restricted maximum likelihood, with a replicate-weight target as an alternative (its Appendices C.2 and C.3). pvstackr 0.2.x does not fit this model, and it computes only the replicate-weight target.
The estimand and the survey-weighted likelihood
The estimand is \(\beta\), the vector of the \(p\) fixed-effect coefficients (the intercept and the slopes) of the model above; \(\sigma\) is a nuisance parameter. For each coefficient \(\beta_\ell\), \(\ell = 1, \dots, p\), pvstackr reports an estimate, a standard error, degrees of freedom and an interval.
The survey weights enter through the likelihood. Let \(\overline W\) be the mean of the final weights and \(\tilde w_i = W_i / \overline W\) the normalized weight of student \(i\), which has mean 1. The survey-weighted likelihood for plausible value \(m\) is
\[ \tilde L_m(\beta, \sigma) = \prod_{i=1}^{N} \phi\big(y_i^{(m)};\, x_i^\top \beta,\, \sigma^2\big)^{\tilde w_i}, \]
where \(\phi(\cdot\,; \mu,
\sigma^2)\) is the normal density with mean \(\mu\) and variance \(\sigma^2\): each student’s log-likelihood
contribution is multiplied by \(\tilde
w_i\). This is the preprint’s weighted likelihood (its Eq. (3),
whose weights are also normalized to mean 1) for a model without the
random intercept. For any fixed \(\sigma\), \(\tilde L_m\) is largest, as a function of
\(\beta\), at the weighted
least-squares estimate \(\hat\beta_m\),
the estimate that pv_brr_target() computes for plausible
value \(m\) with the final weights; the
scale of the weights does not change it. pv_brr_target()
also estimates the covariance \(\hat
U_m\) of \(\hat\beta_m\) from
the replicate weights and combines the \(M\) results with Rubin’s rules into \(\bar\beta\), the average of the \(\hat\beta_m\), and its covariance \(T_{\text{MI}}\) (The design-based target: BRR–Fay
replicate weights and Rubin’s rules).
The stacked fit of stack_direct uses the stacked data:
one copy of the data for each plausible value, \(N \cdot M\) rows in all, where row \(i\) of copy \(m\) has the outcome \(y_i^{(m)}\) and the weight \(\tilde w_i / M\). The weights of each copy
sum to \(N/M\), and those of all rows
to \(N\). The bundled brms engine
multiplies the log-likelihood of each row by its weight, so the log
likelihood of the stacked fit is the average of the \(M\) survey-weighted log-likelihoods, the
stacked objective, and its posterior is the stacked fractional posterior
of the preprint (its Eq. (4)),
\[ q(\beta, \sigma) \propto p(\beta, \sigma) \prod_{m=1}^{M} \tilde L_m(\beta, \sigma)^{1/M}, \]
where \(p(\beta, \sigma)\) is the
prior. A fitting function of your own receives the weights in the
formula term weights(.pvstackr_weight) and the data column
.pvstackr_weight, and has to use them in the same way (The full analysis workflow).
The stacked fixed-effect point identity gives conditions under which the
maximizer of the stacked objective in \(\beta\) equals \(\bar\beta\); One stacked fit and the
Rubin point estimate states them and applies them to pvstackr’s
model. The identity concerns the point estimate only. The variance comes
from the target: the calibration gives the fixed-effect draws of the
stacked fit the target’s mean \(\bar\beta\) and covariance \(T_{\text{MI}}\) (Calibrating the fixed-effect draws
(CCC)).
Symbols and formulas
The Method articles use the symbols below, and each article defines them again where it first uses them. The last column gives the argument, field or column of pvstackr that holds the quantity.
Symbols
The indices:
| Symbol | Meaning | Range |
|---|---|---|
| \(i\) | student (a row of the data) | \(1, \dots, N\) |
| \(m\) | plausible value | \(1, \dots, M\) |
| \(r\) | replicate weight | \(1, \dots, R\) |
| \(\ell\) | fixed-effect coefficient | \(1, \dots, p\) |
| \(s\) | posterior draw | \(1, \dots, S\) |
| \(j\) | school, where a school mean is discussed | — |
The letter \(k\) is not an index: it denotes the Fay coefficient, and the Pareto shape estimate of importance sampling is always written Pareto \(\hat k\).
The data, the weights and the model:
| Symbol | Meaning | In pvstackr |
|---|---|---|
| \(y_i^{(m)}\) | outcome of student \(i\) under plausible value \(m\) | plausible-value columns such as PV1READ
(OUTCOME in the formula) |
| \(x_i\), \(X\) | covariate row of student \(i\); the \(N \times p\) design matrix, the same for every \(m\) | model matrix of the formula |
| \(\beta\), \(\beta_\ell\) | fixed-effect coefficients; coefficient \(\ell\) | rows of the estimate table (term:
b_Intercept, b_x, …) |
| \(\sigma\) | residual standard deviation, not reported | — |
| \(W_i\) | final survey weight of student \(i\) | weight_col |
| \(\overline W\) | mean of the final weights | — |
| \(\tilde w_i = W_i / \overline W\) | normalized final weight (mean 1) | — |
| \(\tilde w_i / M\) | weight of row \(i\) in each of the \(M\) stacked copies | column .pvstackr_weight of the stacked
data |
| \(W_i^{(r)}\) | replicate weight \(r\) of student \(i\) | rep_weight_cols |
| \(k\) | Fay coefficient |
fay_k (PISA: 0.5) |
| \(\tilde L_m(\beta, \sigma)\) | survey-weighted likelihood for plausible value \(m\) | — |
The target, the stacked fit and the calibration:
| Symbol | Meaning | In pvstackr |
|---|---|---|
| \(\hat\beta_m\) | weighted least-squares estimate for plausible value \(m\) with the final weights |
per_pv[[m]]$beta in the target |
| \(\hat\beta_m^{(r)}\) | the same with replicate weight \(r\) | columns of per_pv[[m]]$replicate_beta
|
| \(\hat U_m\) | BRR–Fay replicate covariance of \(\hat\beta_m\) | per_pv[[m]]$U |
| \(\bar\beta\) | combined estimate, the average of the \(\hat\beta_m\) |
beta, beta_bar
|
| \(\bar U\) | average covariance within plausible values | U_bar |
| \(B\) | covariance between plausible values | B |
| \(T_{\text{MI}}\) | total covariance of \(\bar\beta\) | T_MI |
| \(\lambda_\ell\) | fraction of missing information of coefficient \(\ell\) (Barnard and Rubin’s approximation) |
fmi, lambda (same
values) |
| \(\mathrm{riv}_\ell\) | relative increase in variance | riv |
| \(\nu_\ell\) | classic degrees of freedom |
df with
df_method = "classic"
|
| \(\nu_{\text{com}}\) | complete-data degrees of freedom, stated by the user | df_complete |
| \(\nu_{\text{obs},\ell}\), \(\nu_{\text{BR},\ell}\) | Barnard–Rubin observed-data and adjusted degrees of freedom |
df with
df_method = "barnard_rubin"
|
| \(q(\beta, \sigma)\) | stacked fractional posterior | — |
| \(\beta_s\) | fixed-effect draw \(s\) of the stacked fit, before calibration | — |
| \(\bar\beta^{\text{raw}}\), \(\Sigma_{\text{raw}}\) | mean and covariance of these draws | — |
| \(L_{\text{raw}}\), \(L_{\text{tgt}}\) | lower Cholesky factors of \(\Sigma_{\text{raw}}\) and \(T_{\text{MI}}\) | — |
| \(A\) | calibration matrix | — |
| \(\beta^{\text{cal}}_s\) | calibrated draw | rows of get_draws(fit)
|
| \(\kappa_A\) | condition number of \(A\), the ratio of its largest to its smallest singular value |
kappa_A in
get_diagnostics(fit)$ccc
|
| \(\Delta_c\) | center separation |
delta_c_max in
get_diagnostics(fit)$ccc
|
| \(\omega_s^{(m)}\) | importance weight of draw \(s\) for plausible value \(m\), supplied by the user; pvstackr divides each column by its sum |
psis_weights (one row per draw, one column
per plausible value) |
| \(\hat k_m\) | Pareto \(\hat k\) for plausible value \(m\), supplied by the user | pareto_k |
Choosing a
method and reading its intervals also uses \(\hat\beta_m\) and \(\hat U_m\) for the per-PV estimate and
covariance of per_pv and stack_psis.
Main formulas
The target of pv_brr_target() is derived in The design-based target: BRR–Fay replicate
weights and Rubin’s rules: the BRR–Fay replicate covariance of each
\(\hat\beta_m\),
\[ \hat U_m = \frac{1}{R(1-k)^2} \sum_{r=1}^{R} \big(\hat\beta_m^{(r)} - \hat\beta_m\big)\big(\hat\beta_m^{(r)} - \hat\beta_m\big)^\top , \]
Rubin’s rules,
\[ \bar\beta = \frac{1}{M} \sum_{m=1}^{M} \hat\beta_m, \qquad \bar U = \frac{1}{M} \sum_{m=1}^{M} \hat U_m, \]
\[ B = \frac{1}{M-1} \sum_{m=1}^{M} \big(\hat\beta_m - \bar\beta\big)\big(\hat\beta_m - \bar\beta\big)^\top, \qquad T_{\text{MI}} = \bar U + \Big(1 + \frac{1}{M}\Big) B , \]
and, for each coefficient \(\ell\), the fraction of missing information, the relative increase in variance, and the classic and Barnard–Rubin degrees of freedom:
\[ \lambda_\ell = \frac{(1 + 1/M)\, B_{\ell\ell}}{T_{\text{MI},\ell\ell}}, \qquad \mathrm{riv}_\ell = \frac{(1 + 1/M)\, B_{\ell\ell}}{\bar U_{\ell\ell}}, \qquad \nu_\ell = \frac{M - 1}{\lambda_\ell^2}, \]
\[ \nu_{\text{obs},\ell} = \frac{\nu_{\text{com}} + 1}{\nu_{\text{com}} + 3}\, \nu_{\text{com}} (1 - \lambda_\ell), \qquad \nu_{\text{BR},\ell} = \Big(\frac{1}{\nu_\ell} + \frac{1}{\nu_{\text{obs},\ell}}\Big)^{-1} . \]
The stacked fractional posterior \(q(\beta, \sigma)\) is given above; One stacked fit and the Rubin point estimate derives the stacked objective and the point identity.
The calibration map and the center separation are derived in Calibrating the fixed-effect draws
(CCC). With center = "target", which
stack_direct requires, each fixed-effect draw of the
stacked fit is mapped to
\[ \beta^{\text{cal}}_s = \bar\beta + A\big(\beta_s - \bar\beta^{\text{raw}}\big), \qquad A = L_{\text{tgt}} L_{\text{raw}}^{-1}, \]
where \(L_{\text{raw}} L_{\text{raw}}^\top = \Sigma_{\text{raw}}\) and \(L_{\text{tgt}} L_{\text{tgt}}^\top = T_{\text{MI}}\), so that the calibrated draws have mean \(\bar\beta\) and covariance \(T_{\text{MI}}\). Before the calibration, pvstackr measures how far the stacked fit is from the target by
\[ \Delta_c = \max_\ell \frac{\lvert \bar\beta_\ell - \bar\beta^{\text{raw}}_\ell \rvert} {\sqrt{T_{\text{MI},\ell\ell}}} , \]
stored as delta_c_max; delta_c_rel is the
root mean square of the same ratios.
The pooling of per_pv and stack_psis, with
the importance weights \(\omega_s^{(m)}\) and the Pareto \(\hat k_m\) of stack_psis, is
given in Choosing a
method and reading its intervals.
The symbols on the example data
pvstackr includes a small synthetic data set, pisa_tiny,
and a stack_direct fit made from it; Getting started with
pvstackr describes both. The chunks below only read the data and the
target of the fit; they estimate nothing.
pisa_tiny <- read.csv(
system.file("extdata", "pisa_tiny.csv", package = "pvstackr")
)
dim(pisa_tiny) # 12 students, 12 columns
#> [1] 12 12
names(pisa_tiny)
#> [1] "CNT" "CNTSCHID" "CNTSTUID" "x" "female"
#> [6] "PV1READ" "PV2READ" "W_FSTUWT" "W_FSTURWT1" "W_FSTURWT2"
#> [11] "W_FSTURWT3" "W_FSTURWT4"The rows are 12 students in three schools of one made-up country
(CNT is "SYN"), so \(N = 12\); CNTSCHID and
CNTSTUID identify the schools and the students. The
covariates are x and the 0/1 indicator female.
PV1READ and PV2READ are the \(M = 2\) plausible values \(y_i^{(1)}\) and \(y_i^{(2)}\), W_FSTUWT is the
final weight \(W_i\), and
W_FSTURWT1 to W_FSTURWT4 are the \(R = 4\) replicate weights \(W_i^{(r)}\). The replicate weights were
typed in for the example rather than formed from half-samples by Fay’s
method; fay_k = 0.5 is the value declared for them.
detect_pisa_pv_columns() and
detect_pisa_brr_replicate_weights() find the
plausible-value and replicate-weight columns by their names:
detect_pisa_pv_columns(pisa_tiny, suffix = "READ") # the M = 2 plausible values
#> [1] "PV1READ" "PV2READ"
detect_pisa_brr_replicate_weights(pisa_tiny) # the R = 4 replicate weights
#> [1] "W_FSTURWT1" "W_FSTURWT2" "W_FSTURWT3" "W_FSTURWT4"The plausible-value names carry the subject suffix READ,
as in PISA 2022 files, so the call gives suffix = "READ".
With the default suffix "",
detect_pisa_pv_columns() looks for bare names such as
PV1 and stops with an error on these data.
The target of the example fit records the same design and the names of the coefficients:
fit <- readRDS(
system.file("extdata", "examples", "pisa_tiny_stack_direct.rds",
package = "pvstackr")
)$fit
tg <- get_target(fit)
c(M = tg$M, R = tg$R, fay_k = tg$fay_k) # M, R and the Fay coefficient k
#> M R fay_k
#> 2.0 4.0 0.5
tg$fe_names # the p = 3 fixed-effect coefficients
#> [1] "b_Intercept" "b_x" "b_female"The three coefficients b_Intercept, b_x and
b_female of the formula OUTCOME ~ x + female
make up \(\beta\) in this example, and
the estimate table of the fit has one row for each.
The next Method article, The design-based target: BRR–Fay replicate weights and Rubin’s rules, derives the target from these quantities.