PSinference: Exact Inference for Multiply Released Plug-in Sampling Synthetic Data

Ricardo Moura, Mina Norouzirad, Vítor Augusto, Miguel Fonseca

2026-07-13

1 Background and Motivation

1.1 Statistical Disclosure Control and Synthetic Data

Statistical agencies collect sensitive individual-level data from surveys and administrative registers. Before releasing such data publicly they must apply statistical disclosure control (SDC) to prevent identification of individual respondents while preserving statistical utility.

Synthetic data replaces original confidential records with artificial values drawn from a model fitted to the original data. The released dataset contains no real observations, so individual privacy is protected by construction.

1.2 The Plug-in Sampling Framework

Plug-in sampling (PS) is a simple and tractable synthetic data mechanism. The unknown population parameters \(\boldsymbol{\mu}\) and \(\boldsymbol{\Sigma}\) are replaced by their sample counterparts \(\bar{\mathbf{x}}\) and \(\hat{\boldsymbol{\Sigma}}\), and synthetic observations are drawn from the resulting estimated model. Because the plug-in estimates are treated as fixed, the conditional distribution of the synthetic sufficient statistics given the original data is analytically tractable, enabling exact finite-sample inference.

1.3 Single versus Multiple Releases

When \(M \geq 1\) synthetic datasets are released, they may be stacked into: \[\mathbf{V}_{\mathrm{M}} = \begin{pmatrix} \mathbf{V}_1 \\ \vdots \\ \mathbf{V}_M \end{pmatrix} \in \mathbb{R}^{Mn \times p}\]

The stacked synthetic Wishart satisfies: \[(n-1)\mathbf{S}^\star_{\mathrm{M}} \mid \mathbf{S} \sim \mathcal{W}_p\!\left(\mathbf{S}, Mn-1\right)\]

This is the key result: the synthetic Wishart has degrees of freedom \(Mn-1\) while the original sample Wishart retains \(n-1\). Setting \(M=1\) recovers the single-release procedures of Klein et al. (2021).


2 The Brittany Soil Dataset

2.1 SDC Scenario

The package ships with brittany_soil_ps, a real dataset from a soil survey of cultivated fields on loess-silt parent material in Brittany, France (Morvan et al., 2023). It motivates a realistic SDC scenario:

This gives a natural \(p = 6\) dataset with a meaningful two-block structure (\(p_1 = p_2 = 3\)).

2.2 Loading and Inspecting the Dataset

data(brittany_soil_ps)
cat("Dimensions:", nrow(brittany_soil_ps), "x", ncol(brittany_soil_ps), "\n")
#> Dimensions: 37 x 6
cat("Variables :", colnames(brittany_soil_ps), "\n\n")
#> Variables : pH_water pH_KCl log_CEC_Metson log_Organic_C log_Total_N log_P_Olsen

round(head(brittany_soil_ps, 4), 3)
#>    pH_water pH_KCl log_CEC_Metson log_Organic_C log_Total_N log_P_Olsen
#> 2    -1.633 -1.829          0.762         0.872       1.207      -1.279
#> 3     0.756  0.655          0.762         0.036      -0.159       0.308
#> 4    -0.397 -0.442          0.335        -0.124      -0.504       0.071
#> 12   -1.130 -1.384         -0.062         0.047       0.092      -0.979

All variables are standardized. The correlation matrix reveals two highly correlated within-block structures and notable cross-block correlations through log_CEC_Metson:

round(cor(brittany_soil_ps), 2)
#>                pH_water pH_KCl log_CEC_Metson log_Organic_C log_Total_N
#> pH_water           1.00   0.98           0.34         -0.08       -0.02
#> pH_KCl             0.98   1.00           0.32         -0.04        0.00
#> log_CEC_Metson     0.34   0.32           1.00          0.62        0.73
#> log_Organic_C     -0.08  -0.04           0.62          1.00        0.95
#> log_Total_N       -0.02   0.00           0.73          0.95        1.00
#> log_P_Olsen        0.61   0.61           0.50          0.04        0.06
#>                log_P_Olsen
#> pH_water              0.61
#> pH_KCl                0.61
#> log_CEC_Metson        0.50
#> log_Organic_C         0.04
#> log_Total_N           0.06
#> log_P_Olsen           1.00

2.3 Checking Multivariate Normality

Before generating synthetic data it is essential to verify that the original data satisfy the multivariate normality assumption required by the PS procedures.

mvn_result <- mvn_test(brittany_soil_ps, hz_nsim = 1000)
Multivariate normality diagnostic panel for brittany_soil_ps. Histograms (blue = Shapiro-Wilk test not rejected) with fitted normal curves, and chi-square Q-Q plot.
Multivariate normality diagnostic panel for brittany_soil_ps. Histograms (blue = Shapiro-Wilk test not rejected) with fitted normal curves, and chi-square Q-Q plot.
#> 
#> Multivariate Normality Assessment
#> ------------------------------------------------------------ 
#>   n = 37 observations | p = 6 variables
#> ------------------------------------------------------------ 
#> 
#> 1. Univariate Shapiro-Wilk Tests
#>        variable statistic p.value          decision
#>        pH_water    0.9679  0.3550 Fail to Reject H0
#>          pH_KCl    0.9766  0.6142 Fail to Reject H0
#>  log_CEC_Metson    0.9549  0.1393 Fail to Reject H0
#>   log_Organic_C    0.9573  0.1653 Fail to Reject H0
#>     log_Total_N    0.9773  0.6378 Fail to Reject H0
#>     log_P_Olsen    0.9753  0.5683 Fail to Reject H0
#> 
#> 2. Mardia Skewness Test
#>    b1p = 7.6131 | chi-sq(df=56) = 46.9472 | p = 0.8002 | Fail to Reject H0
#> 
#> 3. Mardia Kurtosis Test
#>    b2p = 43.5375 | z = -1.3852 | p = 0.1660 | Fail to Reject H0
#> 
#> 4. Henze-Zirkler Omnibus Test
#>    HZ = 0.02569 | beta = 1.1415 | p = 0.3509 | Fail to Reject H0
#> 
#> 5. Royston H Test
#>    H = 1.0683 | eff. df = 2.57 | p = 0.7103 | Fail to Reject H0
#> 
#> ------------------------------------------------------------
#> Overall: Multivariate normality not rejected

All five tests fail to reject multivariate normality:

print(mvn_result)
#> MVN Assessment: Multivariate normality not rejected (n = 37, p = 6)

3 Generating Synthetic Releases

The function simSynthData(X, M) generates \(M\) PS synthetic releases stacked into an \(Mn \times p\) matrix. The agency decides to release \(M = 3\) independent synthetic copies:

set.seed(2026)
X <- brittany_soil_ps
n <- nrow(X) # 37
p <- ncol(X) # 6

# M = 1: single release
V1 <- simSynthData(X)
cat("Single release: ", nrow(V1), "x", ncol(V1), "\n")
#> Single release:  37 x 6

# M = 3: three stacked releases
V3 <- simSynthData(X, M = 3)
cat("Three releases: ", nrow(V3), "x", ncol(V3), "\n")
#> Three releases:  111 x 6

4 Exact Inferential Procedures

We apply all four exact PS procedures using the named block interface so that variable names appear directly in the output. For each test we show results for both \(M = 1\) (single release) and \(M = 3\) (three releases) side by side to illustrate how inference sharpens as more synthetic datasets are released.

4.1 Generalized Variance

Tests \(\mathcal{H}_0: |\Sigma| = |\Sigma_0|\) (two-sided), using \(\Sigma_0 = \hat{\Sigma}\) from the original data.

op <- par(mfrow = c(1, 2))
gv1 <- gv_test(V1, M = 1, Sigma = cov(X), iterations = 5000)
plot(gv1, main = "Gen. Variance  (M = 1)")
gv3 <- gv_test(V3, M = 3, Sigma = cov(X), iterations = 5000)
plot(gv3, main = "Gen. Variance  (M = 3)")
GV null distribution (log10 scale) for M = 1 (left) and M = 3 (right). The null distribution narrows with M; the observed statistic falls within the non-rejection region in both cases.
GV null distribution (log10 scale) for M = 1 (left) and M = 3 (right). The null distribution narrows with M; the observed statistic falls within the non-rejection region in both cases.
par(op)
cat("M = 1: p =", gv1$p.value, "|", gv1$decision, "\n")
#> M = 1: p = 0.4624 | Fail to reject H0
cat("M = 3: p =", gv3$p.value, "|", gv3$decision, "\n")
#> M = 3: p = 0.398 | Fail to reject H0

The test fails to reject for both \(M = 1\) and \(M = 3\): the PS mechanism preserves the overall covariance scale regardless of the number of releases.

4.2 Sphericity Test

Tests \(\mathcal{H}_0: \Sigma = \sigma^2 I_6\). Given the strong within-block correlations and different scales, strong rejection is expected.

op <- par(mfrow = c(1, 2))
sph1 <- sphericity_test(V1, M = 1, iterations = 5000)
plot(sph1, main = "Sphericity  (M = 1)")
sph3 <- sphericity_test(V3, M = 3, iterations = 5000)
plot(sph3, main = "Sphericity  (M = 3)")
Sphericity null distribution for M = 1 (left) and M = 3 (right). The observed statistic is far into the left rejection region in both cases; larger M concentrates the null distribution.
Sphericity null distribution for M = 1 (left) and M = 3 (right). The observed statistic is far into the left rejection region in both cases; larger M concentrates the null distribution.
par(op)
cat(
  "M = 1: stat =", round(sph1$statistic, 4),
  "| p =", sph1$p.value, "|", sph1$decision, "\n"
)
#> M = 1: stat = 0.2502 | p = 0 | Reject H0
cat(
  "M = 3: stat =", round(sph3$statistic, 4),
  "| p =", sph3$p.value, "|", sph3$decision, "\n"
)
#> M = 3: stat = 0.2417 | p = 0 | Reject H0

Sphericity is decisively rejected for both \(M = 1\) and \(M = 3\). The null distribution is visibly more concentrated for \(M = 3\), reflecting the increased effective sample size \(N = Mn = 111\).

4.3 Independence Test

Tests \(\mathcal{H}_0: \Sigma_{12} = \mathbf{0}\) — are the public soil indicators independent of the sensitive farm indicators?

op <- par(mfrow = c(1, 2))
ind1 <- independence_test(V1,
  M = 1,
  group_a = c("pH_water", "pH_KCl", "log_CEC_Metson"),
  group_b = c("log_Organic_C", "log_Total_N", "log_P_Olsen"),
  iterations = 5000
)
plot(ind1, main = "Independence  (M = 1)")
ind3 <- independence_test(V3,
  M = 3,
  group_a = c("pH_water", "pH_KCl", "log_CEC_Metson"),
  group_b = c("log_Organic_C", "log_Total_N", "log_P_Olsen"),
  iterations = 5000
)
plot(ind3, main = "Independence  (M = 3)")
Independence null distribution for M = 1 (left) and M = 3 (right). Rejection is clear in both cases; the critical value moves right as M increases.
Independence null distribution for M = 1 (left) and M = 3 (right). Rejection is clear in both cases; the critical value moves right as M increases.
par(op)
cat(
  "M = 1: stat =", round(ind1$statistic, 4),
  "| p =", ind1$p.value, "|", ind1$decision, "\n"
)
#> M = 1: stat = 0.1427 | p = 0 | Reject H0
cat(
  "M = 3: stat =", round(ind3$statistic, 4),
  "| p =", ind3$p.value, "|", ind3$decision, "\n"
)
#> M = 3: stat = 0.0741 | p = 0 | Reject H0

Independence is rejected in both cases. The stronger evidence at \(M = 3\) reflects the fact that more synthetic data narrows the null distribution and lowers the critical value, making it easier to detect the true cross-block correlations.

4.4 Regression Test

Tests \(\mathcal{H}_0: \Delta = \mathbf{0}\) — can the public indicators be predicted from the sensitive block?

op <- par(mfrow = c(1, 2))
reg1 <- regression_test(V1,
  M = 1,
  response = c("pH_water", "pH_KCl", "log_CEC_Metson"),
  predictors = c("log_Organic_C", "log_Total_N", "log_P_Olsen"),
  Delta0 = matrix(0, 3, 3),
  iterations = 5000
)
plot(reg1, main = "Regression  (M = 1)")
reg3 <- regression_test(V3,
  M = 3,
  response = c("pH_water", "pH_KCl", "log_CEC_Metson"),
  predictors = c("log_Organic_C", "log_Total_N", "log_P_Olsen"),
  Delta0 = matrix(0, 3, 3),
  iterations = 5000
)
plot(reg3, main = "Regression  (M = 3)")
Regression null distribution (log10 scale) for M = 1 (left) and M = 3 (right). Both reject the zero-regression null.
Regression null distribution (log10 scale) for M = 1 (left) and M = 3 (right). Both reject the zero-regression null.
par(op)
cat(
  "M = 1: stat =", round(reg1$statistic, 5),
  "| p =", reg1$p.value, "|", reg1$decision, "\n"
)
#> M = 1: stat = 0.02876 | p = 0.0124 | Reject H0
cat(
  "M = 3: stat =", round(reg3$statistic, 5),
  "| p =", reg3$p.value, "|", reg3$decision, "\n"
)
#> M = 3: stat = 0.73265 | p = 0 | Reject H0

The zero-regression null is rejected for both \(M = 1\) and \(M = 3\), confirming that partial reconstruction of the public indicators from the sensitive block is statistically detectable.


5 Effect of \(M\) on Inference

The side-by-side plots above already show the sharpening of inference as \(M\) increases from 1 to 3. The following panel extends this to \(M \in \{1, 2, 3, 5\}\) for the sphericity test:

op <- par(mfrow = c(2, 2))
for (m in c(1L, 2L, 3L, 5L)) {
  Vm <- simSynthData(X, M = m)
  res <- sphericity_test(Vm, M = m, iterations = 2000)
  plot(res, main = sprintf("Sphericity: M = %d  (N = %d)", m, m * n))
}
Sphericity null distribution for M = 1, 2, 3, 5. The distribution concentrates as M grows, lowering the critical value and strengthening evidence against the false null.
Sphericity null distribution for M = 1, 2, 3, 5. The distribution concentrates as M grows, lowering the critical value and strengthening evidence against the false null.
par(op)
# Numerical summary across M
results <- lapply(c(1L, 2L, 3L, 5L), function(m) {
  Vm <- simSynthData(X, M = m)
  s <- sphericity_test(Vm, M = m, iterations = 3000)
  i <- independence_test(Vm,
    M = m,
    group_a = c("pH_water", "pH_KCl", "log_CEC_Metson"),
    group_b = c("log_Organic_C", "log_Total_N", "log_P_Olsen"),
    iterations = 3000
  )
  data.frame(
    M = m,
    N = m * n,
    sph_stat = round(s$statistic, 5),
    sph_pval = round(s$p.value, 5),
    ind_stat = round(i$statistic, 5),
    ind_pval = round(i$p.value, 5)
  )
})
do.call(rbind, results)
#>   M   N sph_stat sph_pval ind_stat ind_pval
#> 1 1  37  0.24758        0  0.10415        0
#> 2 2  74  0.23357        0  0.05643        0
#> 3 3 111  0.24090        0  0.07909        0
#> 4 5 185  0.23820        0  0.09270        0

As \(M\) increases the evidence against both sphericity and independence strengthens monotonically, with \(p\)-values decreasing and test statistics moving further from the null.


6 Summary

What do the Results Tell the Agency?

Test \(M=1\) \(M=3\) Interpretation
Generalized variance Fail to reject Fail to reject PS preserves overall covariance scale
Sphericity Reject Reject Soil variables are highly correlated, unequal variance
Independence Reject Reject Public and sensitive blocks are correlated
Regression Reject Reject Public indicators are partially predictable from sensitive block

The consistency across \(M = 1\) and \(M = 3\) confirms that the conclusions are stable: they are a property of the original data structure, not of the specific synthetic realisation. The stronger \(p\)-values at \(M = 3\) reflect increased inferential precision from more synthetic information, not a change in the underlying finding.


7 References

Klein, M., Moura, R., and Sinha, B. (2021). Multivariate normal inference based on singly imputed synthetic data under plug-in sampling. Sankhya B, 83, 273–287.

Morvan, T., Lambert, Y., Germain, P., Lemercier, B., Moreira, M. and Beff, L. (2023). A dataset of physico-chemical properties, extractable organic N, N mineralisation and physical organic matter fractionation of soils. Data in Brief, 51, 109776.

Raghunathan, T. E., Reiter, J. P., and Rubin, D. B. (2003). Multiple imputation for statistical disclosure limitation. Journal of Official Statistics, 19, 1–16.

Rubin, D. B. (1993). Statistical disclosure limitation. Journal of Official Statistics, 9, 461–468.