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.
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.
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).
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\)).
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.979All 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.00Before generating synthetic data it is essential to verify that the original data satisfy the multivariate normality assumption required by the PS procedures.
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:
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 6We 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.
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)")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 H0The 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.
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)")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 H0Sphericity 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\).
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)")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 H0Independence 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.
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)")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 H0The 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.
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))
}# 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 0As \(M\) increases the evidence against both sphericity and independence strengthens monotonically, with \(p\)-values decreasing and test statistics moving further from the null.
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.
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.