Logistic Regression With Variables Subject To Post Randomization 

Method 

Yong Ming Jeffrey Woo ∗, Aleksandra B. Slavkovic ́

Department of Statistics, The Pennsylvania State University 

Abstract

An increase in quality and detail of publicly available databases increases the risk of disclosure of sensitive 

personal information contained in such databases. The goal of Statistical Disclosure Control (SDC) is to develop 

methodology that aims at minimizing disclosure risk while providing society with as much information as possible 

needed for valid statistical inference. The Post Randomization Method (PRAM) is a disclosure avoidance method, 

where values of categorical variables are perturbed via some known probability mechanism, and only the perturbed 

data are released thus raising issues regarding disclosure risk and data utility. In this paper, we propose a number of 

EM algorithms to obtain unbiased estimates of the logistic regression model with data subject to PRAM, and thus 

effectively account for the effects of PRAM and preserve data utility. The effect of the level of perturbation and 

sample size on the estimates are evaluated, and relevant standard error estimates are proposed. 

Keywords: Statistical disclosure control, logistic regression, generalized linear models, EM algorithm 

∗corresponding author email: yjw102@psu.edu 
1 Introduction

The goal of Statistical Disclosure Control (SDC), and related ®elds such as Privacy-Preserving Data-Mining, is to 

provide society with as much information as possible while individual information is suf®ciently protected against 

public disclosure. Development of SDC methodology and its practice are important to of®cial statistics for many 

reasons,a numberof which are discussed by Willenborg and de Waal (1996), Fienberg and Slavkovic ́ (2010), and 

Ramanayake and Zayatz (2010). For example, there are laws requiring the protectionof the con®dentialityof respondents’ information by statistical of®ces such as Title 13of the United States Code that outlines the roleof the U.S. 

Census. Disclosure control is also needed to maintain trust amongst respondents in order to obtain adequate response 

rates. Concerns about threats to data privacy have even led to the abolition of the census in the past (e.g., Netherlands 

in 1971; see van der Laan (2000)). 

Traditionally, publicly available data have been mostly inthe aggregate form, but nowadays there exists substantial demand for high-quality detailed data products. Microdata are sets of records containing detailed information on 

individual respondents and many SDC methods have been developed for microdata. When SDC methods are applied 

to a dataset they lead to a publication of altered datasets, which would be available for a wider use and have reduced 

disclosure risk of sensitive information but also reduced data utility. The Post Randomization Method (PRAM) is 

a disclosure avoidance method originally proposed by Gouweleeuw et al. (1998). The main idea behind PRAM is 

to publish redacted data after the values of categorical variables in the original data have been misclassi®ed by a 

known probability mechanism.This probabilitymechanismis described bya transitionmatrix,socalled aPRAMmatrix. While PRAM provides certain advantages compared to other SDC methodologies (e.g., unlike non-perturbative 

methods which lead to loss of detail, PRAM maintains detail in the variables; unlike most perturbative methods, 

the application of PRAM is probabilistic in nature and hence intruders cannot determine which records have been 

perturbed), it has seen a limited use in practice. One of the commonly-discussed issue in the literature is of model 

parameter estimation when data are subject to PRAM. Summary statistics from PRAMed data are typically biased 

and need to be adjusted to take the effects of PRAM into account. Current literature has focused on a speci®c subset 

of problems; for example, Gouweleeuw et al. (1998) proposed an unbiased moment estimator for frequency counts; 

van den Hout and van der Heijden (2002) proposed formulas to estimate odds ratios for data subject to PRAM; and 

vandenHoutand Kooiman(2006)proposedan EMalgorithmto estimatethelinearregressionmodelwithcovariates 

subject to PRAM. 

In this paper, we focus on measuring data utility from the standpoint of statistical inference (e.g., see Slavkovic ́ 

and Lee (2010)) and propose a way to obtain unbiased estimators of parameters in logistic regression models when 

data havebeen subjectedto PRAM. We developand implement EM-type algorithmsto obtain asymptotically unbiased 

estimators, that is the maximum likelihood estimators of parameters in logistic regression models, when variables are 

subject to PRAM. The basic ideas are based on the aEM by method of weights” developed by Ibrahim (1990) for 

generalized linear models(GLMs)withcovariates missingat random,andonthe approach proposed by vandenHout 

and Kooiman (2006) for linear regression with covariates subject to PRAM. There is an extensive literature for missing data with either missing covariate or missing response variable. We extend these ideas by developing an EM-type 

algorithm that obtains unbiased estimates of logistic regression when both covariate and response variables are subject 

to PRAM. This is a more dif®cult problem than either case when covariates or response variables are subject to PRAM, 

and has received little attention in both PRAM and missing data literatures. 

The outline of this paper is as follows. Section2 presents the EM-type methodology to obtain estimates of the 

logisticregressionmodelwhenvariablesaresubjecttoPRAMand reportsthe resultsof simulation studiestoevaluate 

the methodologyforthreedifferent cases:(1)covariatesubjecttoPRAM,(2) response variablesubjecttoPRAM,and 

(3)bothcovariateand response variables subjectto PRAM. Section 3 reportsthe resultsof simulation studies intended 

to evaluate the effects of varying the parameters of the logistic regression model on the proposed methodology. Section 

4 applies the proposed methodology to data from the 1993 Current Population Survey, and Section 5 contains some 

discussion. 

1 
 

 

 

 

 

 

 

 

 

 

 

 

  

 

 

 

 

 

 

 

 

2 Estimating Logistic Regression Model with Variables Subject to PRAM

We ®rst introduce the notation for the standard logistic regression model. Let x = (x0, x1, ..., xp)

t denote the vector 

with p covariates, and let Y be the binary response variable. The logistic regression model is written as 

exp(β

t

x)

E(Y |x) = (1)

1 + exp(β

t

x)

where β = (β0, β1, ..., βp) are the parameter coef®cients. We ®x x0 = 1 so β0 is the intercept. Let μ(x) = P(Y =

1|x, β), so 1 − μ(x) = P(Y = 0|x, β). 

The likelihood function is given by 

Y L(β) = [μ(xi)]yi

[1 − μ(xi)]1−yi, (2) 

i

the loglikelihood is 

Xn

`(β) = {yi

log[μ(xi)]] + (1 − yi)log[1 − μ(xi)]]}, (3) 

i

and the information matrix is given by 

∂

2

`(β) X I(β)jk = − = xijxikμ(xi)(1 − μ(xi)). (4) ∂βj∂βk

i

Similartothe ideasin Ibrahim (1990)and vanden Houtand Kooiman (2006)forcovariates missingat randomand 

for covariates subject to PRAM respectively, we develop and implement three EM-type algorithms to obtain unbiased 

estimators for parameters in the logistic regression model for three cases: (1) covariates subject to PRAM; (2) response 

variables subject to PRAM; and (3) both covariates and response variables subject to PRAM. The E-step consists of 

computing the weights with the derived formulas for the conditional distribution of the true data given the observed 

data and current parameter estimates. The M-step can be carried out using weighted regression with any standard 

statistical software, which makes implementation of the algorithms more convenient. Case 3 has not been studied as 

extensively as Cases (1) and (2) in either PRAM or missing data literature, and is more complex. 

2.1 Categorical Covariate Subject to PRAM

Let X denote the categorical covariate to which PRAM is applied, with X∗ denoting the observed, released versionof 

X. The levelsof X and X∗

are {x1, ..., xJ }. Let PX be the J x J PRAM transition matrix that contains the transition 

∗ probabilities, with pXjk = P(x = xk|X = xj ), π∗ = P(X∗ = xj ) and πj = P(X = xj ).

j

2.1.1 EM Algorithm I

We present an EM algorithm to obtain unbiased estimates of (1) for Case (1) when covariates are subject to PRAM. 

This method is similar to the aEM by method of weights” proposed by Ibrahim (1990), which is used to estimate 

parameters in GLMs with missing covariates. Consider the joint distribution of (xi

, yi), which can be speci®ed via 

the conditional distribution of yi given xi and the distribution of xi

. Then the complete data loglikelihood can be 

expressed as 

Xn

`(φ; X, y) = `(φ; xi, yi)

i=1

Xn

= {`yi|xi

(β) + `xi(π)} (5) 

i=1

where φ = (β, π), and the distribution of X is multinomial with parameter π. The E-step can be written as 

2 
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Xn

Q(φ|φ

(v)

) =

(v) E(`(φ; xi

, yi)|data, φ )

i

Xn X

J

= qj (i)`(φ; xi

, yi)

i=1 j=1

Xn X

J

= qj (i){`yi|xi

(β) + `xi(π)} (6) 

i=1 j=1

where )q

(v

j (i) is the conditional probability of X|observed data, φ for subject i.The ®rst part of (6) is the loglikelihood of the logistic regression model, and second part is the loglikelihood of a multinomial distribution. 

The M-step maximizes (6). This can be done via a weighted logistic regression, by creating a anew” dataset, 

with each subject i having (Xnew = x1),(Xnew = x2), ...,(Xnew = xJ ) with weights qj (i) = P(Xnew(i) =

∗ xj |y(i), x∗

(i), φ(v)), where x is the observed value of the PRAMed variable. Using Bayes’ rule, the conditional 

distribution of Xnew is 

(v)

pXjkP(Y |xj , β(v))πj

P(Xnew = xj |X∗ = xk, Y, φ

(v)

) = PJ(7) (v)

l=1 pXlkP(Y |xl

, β(v))πl

for k, j ∈ {1, ..., J}. An example on how the weighted logistic regressionis carried outis shown below in Table (1). 

Table 1: Methodology: Weighted logistic regression where Y is a response variable, X∗ denotes observed value of 

covariate and Xnew denotes a list of possible values of the covariate. 

Y X∗ Xnew Weight for Xnew

y(1) x

∗

(1) 0 

1 

q0(1)

q1(1)

y(2) x

∗

(2) 0 

1 

q0(2)

q1(2)

. 

. 

. 

y(n) x

∗

(n)

0 

1 

q0(n)

q1(n)

EM AlgorithmI runs as follows: 

EM Algorithm I: Initial values can be the estimates of β from the logistic regression on Y ∼ X∗

, where X∗is the 

PRAMed covariate. π

∗

can be used as the initial estimate of π. 

E-step: 

(v)

Compute q (i) for i = 1, ..., n and j = 1, ..., J.

j

M-step: 

(v)

Carry out weighted logistic regression with weights q (i), using standard software. j

Update φ

(v)

β

(v+1) = βˆ from weighted logistic regression 

(v) (v) (v) ( P v)

π

(v+1) = (q (+), ..., q (+))t

1 J

/n, where qj(+) = i

qi

With the updated φ

(v+1),a new dataset with new weights can be computed in the E-step, and the algorithm continues 

until convergence. 

3 
 

 

 

 

 

 

 

 

 

 

 

 

2.1.2 Comparison of Estimation With & Without Accounting for PRAM

To demonstrate the effect of PRAM on maximum likelihood estimates of parameter coef®cients in logistic regression model (1), a simulation study was carried out, with one binary covariate, β = (1, 2)t

, xi1 sampled from 

exp(β

t

x) Bernoulli(0.4), and yi sampled from Bernoulli(

exp t )1+ (β )

. The estimates for the logistic regression of Y on x

X are reportedin Table2, labeledasaβ

∗

real”. PRAM is then applied to x to obtain x , with the following transition 

matrix 

  p 1 − p

PX = ,

1 − p p

where we vary the value of p. 

Standard logistic regression is performed on y with x

∗

, while x is not released and is considered to be unobserved. We estimate and its approximate con®dence interval, β ˆ ± ˆ

1 95% β1 2SE(β1). We ran 500 simulations, with 

n = 100, 1000, 10000, p = 0.8, 0.9. The mean estimates of β1 and the coverage probabilities of the 95% CIs are 

computed and are reported in Table 2 ina row labeled aβnoadjust”. 

Next, using EM Algorithm I, we obtain the MLEs of model (1). The algorithm ran for 20 steps. The mean estimates 

of β1 and the coverage probabilities of the 95% CIs are computed and reported in Table 2 ina row labeled aβadjust”. 

The inverse of the observed information matrix can be used to estimate the covariance matrix. Following the method 

describedin Ibrahim etal. (2005)and Louis (1982),the estimated observed information matrix is given by 

X X n n

I(βˆ) = −Q ̈(βˆ|βˆv

) − {[ qˆj (i)Si(βˆ|xi, yi)Si(βˆ|xi, yi)

t

] − Q ̇

i(βˆ|βˆv

)Q ̇

i(βˆ|βˆv

)

t

} (8) 

i=1 i=1

Pn PJ ∂

2

 ̈ `(β|xi,yi) PJ ∂`(β|xi,yi) ∂`(β|xi,yi) where Q  ̇ (β|β

v

) = i=1 j=1 qj (i)

∂β∂βt

, Qi(β|β

v

) = j=1 qj (i)

∂β , and Si(β|xi

, yi) = ∂β . 

From Table 2, we can see that estimates of the logistic regression coef®cients without accounting for PRAM are 

biased (see βnoadjust). The bias appears to increase when p decreases (i.e. higher level of perturbation). Actual 

coverage probabilities also decreases as sample size increases and as p decreases. When accounting for PRAM, the 

estimates appear to be close to their true values (see βadjust). The estimates are less biased as sample size increases, 

for example, when p = 0.9, the bias is −0.0031, 0.0023 and 0.002 for n = 100, n = 1000 and n = 10000 respectively. Coverage probabilities decreases with higher level of perturbation, and larger sample sizes. 

Table 2: Simulated ML estimates of (1)with accounting forPRAM. Average ML estimates. Coverage probabilities in 

parentheses 

p = 0.9 p = 0.8

n = 100

βreal

βnoadjust

βadjust

0.6062 

0.4403 (0.986) 

0.6093 (1.000) 

0.6062 

0.2872 (0.872) 

0.6914 (0.998) 

n = 1000

βreal

βnoadjust

βadjust

0.4947 

0.3810 (0.960) 

0.4924 (1.000) 

0.4947 

0.2822 (0.720) 

0.5031 (0.926) 

n = 10000

βreal

βnoadjust

βadjust

0.4995 

0.3897 (0.350) 

0.5015 (0.890) 

0.4995 

0.2872 (0.02) 

0.5042 (0.790) 

Figures1 and2 display the plotsof the truevaluesof β1 (in black), along with the estimates of β1 (in red; left 

column without adjustment, right column using adjustment with EM algorithm), and the con®dence intervals for the 

estimates via the EM algorithm (in dotted green), for n = 10000, p = 0.90 and p = 0.80 respectively. We only show 

plots for the ®rst 50 simulations. When not adjusting for PRAM, the true values of β1 barely fall within the con®dence 

intervals, especially when p = 0.8. When the adjustment is made using EM algorithm I, the true values of β1 are more 

4 
 

 

 

00 

0 

<D 

6 

00 "' 0 

w 

a, .. 

10 

10 

Without Adjustment, p=0.9 

20 30 

Index 

Without Adjustment, p=0.8 

20 30 

Index 

40 

40 

Cl fCN"nolldjUst 

50 

real 

noad;Jst 

Cl l 0< no lldjUst 

50 

<D 

6 

00 

0 

<D 

6 

00 "' 0 

w 

a, .. 

10 

10 

With Adjustment, p=0.9 

20 30 

Index 

With Adjustment, p=0.8 

20 30 

Index 

40 

40 

EM esti'nal e 

Cl withEMestinate 

50 

EU ntinate 

Cl withEM estffllll e 

50 

likely to fall within the con®dence intervals. Plots for n = 100 and n = 1000 are shown in Appendix A, Figures 11 -

14. 

Figure 1: Plots of estimates with covariate subject to PRAM, with 95% CI, when n = 10000, p = 0.90

Figure 2: Plots of estimates with covariate subject to PRAM, with 95% CI, when n = 10000, p = 0.80

Additional simulations to study the effect of the proposed EM algorithm with varying π, β and multinomial 

covariates are discussed in Section 3. Next, we focus on presenting EM algorithms for cases (2) and (3). 

5 
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2.2 Response Variable Subject to PRAM

Let Y denote the binary response variable to which PRAM is applied, with Y

∗ denoting the observed, unreleased 

version of Y . Let PY be the PRAM transition matrix that contains the transition probabilities, with pY jk = P(Y

∗ =

k|Y = j), k, j ∈ {0, 1}. 

2.2.1 EM Algorithm II

Following the method proposed in Section 2.1.1, the parameter π in the complete data loglikelihood (5) can be estimated directly since X is not subject to PRAM. Thus, the E-step simpli®es to 

XX

n J

Q(φ|φ

(v)

) = rj (i){`yi|xi(β)} (9) 

i=1 j=1

where rj (i) is the conditional probability of Y |observed data, β(v)for subject i. 

The M-step maximizes (9). This can be done via a weighted logistic regression, by creating a anew” dataset, with 

each subject i having (Ynew = 0),(Ynew = 1) with weights rj (i) = P(Ynew(i) = j|y

∗

(i), x(i), β(v)). Using Bayes’ 

rule, the weights can be computed as 

pY jkP(Y = j|X, β

(v)

)

P(Ynew = j|Y

∗ = k, X, β(v)

) = P . (10) 

l

pY lkP(Y = l|X, β

(v))

EM Algorithm II runs as follows: 

EM Algorithm II: Initial values can be the estimates of β from the logistic regression on Y

∗ ∼ X, where Y∗

is the 

PRAMed response variable. 

E-step: 

(v)

Compute r (i) for i = 1, ..., n and j = 0, 1.

j

M-step: 

(v)

Carry out weighted logistic regression with weights r (i), using standard software. j

Update β

(v)

β

(v+1) = βˆ from weighted logistic regression 

With the updated β

(v+1), a newdataset with new weights can be computed in the E-step, and thealgorithm continues 

until convergence. 

2.2.2 Comparison of Estimation With & Without Accounting for PRAM

To demonstrate the effectof PRAM on maximum likelihood estimatesof parameter coef®cientsin logistic regression 

model (1), a simulation study was carried out, with one covariate, β = (1, 0.5)t, and xi1 sampled from N ∼ (1, 1), and 

exp(β

t

x)

yi sampled from Bernoulli( ). Like before, we report the estimated parameters of the logistic regression of 1+exp(βtx)

Y on X (see Table3, row labeled as aβreal”). PRAM is then applied to y to obtain y

∗

, with the following transition 

matrix 

 

PY =

p

1 − p

1 − p

p

,

where we vary the value of p. 

∗ Standard logistic regression is performed on y with x, while y is not released and is considered to be unobserved. We estimate β1 and its approximate 95% con®dence interval, βˆ

1 ± 2SE(βˆ1). We ran 500 simulations, with 

6 
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

n = 100, 1000, 10000, p = 0.8, 0.9. The mean estimates of β1 and the coverage probabilities of the 95% CIs are 

computed and are reported in Table 3 ina row labeled aβnoadjust”. 

Next, using EM Algorithm II, we obtain the MLEs of model (1). The algorithm ran for 20 steps. The mean 

estimates of β1 and the coverage probabilities of the 95% CIs are computedand reported in Table 3 in a row labeled 

aβadjust”. The inverse of the observed information matrix can be used to estimate the covariance matrix. Note that 

X I(β) = I(β)i

i

where I(β)iis the contribution of subject i tothe information matrix. Followingthe method describedin Louis (1982) 

IY ∗|X = IY |X − IY |Y ∗,X.

From (4), it follows that 

t

I(β)i = xixiμˆ(xi)(1 − μˆ(xi))

− xix

t

i

[P(Y = 1|Y

∗

, xi).(1 − P(Y = 1|Y

∗

, xi))]. (11) 

From Table 3, we can see that estimates of the logistic regression coef®cients without accounting for PRAM are 

biased (see βnoadjust). The bias appears to increase when p decreases. Actual coverage probabilities also decreases as 

sample size increases and as p decreases. When accounting for PRAM, the estimates appear to be close to their true 

values (see βadjust). The estimates are least biased for p = 0.9 and n = 10000. Coverage probabilities decreases with 

higher level of perturbation, and larger sample sizes. 

Table 3: Simulated ML estimates of (1)with accounting forPRAM. Average ML estimates. Coverage probabilities in 

parentheses 

p = 0.9 p = 0.8

n = 100

βreal

βnoadjust

βadjust

0.5361 

0.3550 (0.922) 

0.5777 (0.970) 

0.5361 

0.2259 (0.754) 

0.6077 (0.854) 

n = 1000

βreal

βnoadjust

βadjust

0.4962 

0.3216 (0.464) 

0.5015 (0.958) 

0.4962 

0.2117 (0.072) 

0.5101 (0.746) 

n = 10000

βreal

βnoadjust

βadjust

0.5016 

0.3229 (0.000) 

0.5007 (0.956) 

0.5016 

0.2118 (0.000) 

0.5007 (0.726) 

Figures 3 and 4 display the plotsof the true valuesof β1, along with the estimates of β1, and the con®dence intervals for the estimates via the EM algorithm, for n = 10000, p = 0.90 and p = 0.80 respectively. Similar to the results 

from Section 2.1.2, the true values of β1 barely fall within the con®dence intervals when not using theEM algorithm. 

When the adjustment is made using EM algorithm II, the true values of β1 are more likely to fall within the con®dence 

intervals. Plots for n = 100 and n = 1000 are shown in Appendix B, Figures 15 - 18. 

2.3 Covariate and Response Subject to PRAM

2.3.1 EM Algorithm III

The weighted logistic regressionis doneby creatinga anew” dataset with each subject i having (Ynew = 0),(Ynew =

1) and (Xnew = x1),(Xnew = x2), ...,(Xnew = xJ ) with weights wml(i) = P(Ynew(i) = m, Xnew = xl|Y

∗

(i) =

k, X∗(i) = xj , β(v)). The weights can be computed as 

7 
 

 

 

 

00 

w 

(Il 

1n 

w 

(Il 

"' 

"' 0 

"' 0 

"' 

"' 0 

N 

0 

10 

10 

Without Adjustment, p=0.9 

20 30 

Index 

Without Adjustment, p=0.8 

,, 

' ' 

40 

\ I \ , \,,,~ \,/ 

20 30 40 

Index 

real 

noadjust 

Cl for no adjust 

real 

noadjust 

50 

Cl for no adjust 

50 

w 

(Il 

1n 

w 

(Il 

N 

0 

"' 0 

"' 

"' 0 

N 

0 

10 

10 

With Adjustment, p=0.9 

20 30 

Index 

With Adjustment, p=0.8 

20 30 

Index 

40 

40 

real 

EMeslm'lle 

CJ with EM estinate 

real 

EMeslinllle 

50 

CJ with EM estinate 

50 

Figure 3: Plots of estimates with response subject to PRAM, with 95% CI, when n = 10000, p = 0.90

Figure 4: Plots of estimates with response subject to PRAM, with 95% CI, when n = 10000, p = 0.80

P pY mkP(Y = m|X = l, β

(v)

) pXljπ(l)

b

pY bkP(Y = b|X = l, β

(v)

)

wml = P P P . (12) 

pY akP(Y = a|X = l, β

(v)) pXcjπ(c) d

pY dkP(Y = d|X = c, β

(v)) a c

The weights (12) are more dif®cult to derive mathematically than the weights (7) and (10). See Appendix C for 

derivations of the weights. 

8 
 

 

 

 

 

 

 

 

 

 

 

 

 

EM Algorithm III runs as follows: 

EM Algorithm III: Initial values can be the estimates of β from the logistic regression on Y

∗ ∼ X∗

, where Y

∗

and 

X∗are the PRAMed response variable and covariate. π canbe estimatedby πˆ = (P

−1

)

tπ∗

X .

E-step: 

(v)

Compute w (i) for i = 1, ..., n and j = 0, 1.

j

M-step: 

(v)

Carry out weighted logistic regression with weights w (i), using standard software. j

Update β

(v)

β

(v+1) = βˆ from weighted logistic regression 

With the updated β

(v+1), a newdataset with new weights can be computed in the E-step, and thealgorithm continues 

until convergence. 

2.3.2 Comparison of Estimation With & Without Accounting for PRAM

To demonstrate the effect of PRAM on maximum likelihood estimates of parameter coef®cients in logistic regression model (1), a simulation study was carried out, with one binary covariate, β = (1, 2)t

, and xi1 sampled from 

exp(β

t

x) Bernoulli(0.4), and yi sampled from Bernoulli( ). The parameter estimates of the logistic regression of 1+exp(βtx)

∗ Y on X are reportedinTable4, labeled as aβreal”. PRAM is then applied to y and x to obtain y and x

∗

, with the 

following transition matrix 



P =

0.9

0.1

0.1

0.9



.

∗ Standard logistic regression is performed on y with x∗

, while y and x are not released and are considered to be 

unobserved. We estimate β1 and its approximate 95% con®dence interval, βˆ

1 ± 2SE(βˆ1). We ran 500 simulations, 

with n = 100, 1000, 10000. The mean estimates of β1 and the coverage probabilities of the 95% CIs are computed 

and are reported in Table 4 ina row labeled aβnoadjust”. 

Next, using EM Algorithm III, we obtain the MLEs of model (1). The algorithm ran for 20 steps. The mean 

estimates of β1 and the coverage probabilities of the 95% CIs are computed and reported in Table 4 in a row labeled aβadjust”. From Table4, we can see that estimatesof the logistic regression coef®cients without accounting 

for PRAM are biased (see βnoadjust). Similar to the results from the previous sections, the bias appears to increase 

when p decreases. Actual coverage probabilities also decreases as sample size increases and as p decreases. When 

accounting for PRAM, the estimates appear to be close to their true values (see βadjust). The estimates are less biased 

as sample size increases, for example, the bias is 1.0045, −0.0687 and 0.025 for n = 100, n = 1000 and n = 10000

respectively. Coverage probabilities decreases with higher level of perturbation, and larger sample sizes. 

Table 4: Simulated ML estimates of (1)with accounting forPRAM. Average ML estimates. Coverage probabilities in 

parentheses 

βreal 4.4725 

n = 100 βnoadjust 0.8872 (0.516) 

βadjust 3.4680 (1.000) 

βreal 2.0302 

n = 1000 βnoadjust 0.80912 (0.000) 

βadjust 2.0989 (0.922) 

βreal 1.9997 

n = 10000 βnoadjust 0.8003 (0.000) 

βadjust 1.9747 (0.822) 

9 
 

 

 

 

 

 

 

 

 

 

 

 

 

 

1n w 

l;J 

(D 

"' 

D 

D 

"' 

D 

10 

Without Adjustment, p=0.9 

20 30 

Index 

40 

real 

no adjust 

Cl for no adjust 

50 

1n 

w 

m 

., 

(D 

D 

q 

"' 

q 

10 

With Adjustment, p=0.9 

20 30 

Index 

40 

real 

EM estinate 

Cl withEMestinate 

50 

Figure 5 displays the plotsof the true valuesof β1, along with the estimates of β1, and the con®dence intervals 

for the estimates via the EM algorithm, for n = 10000, p = 0.90. Similar to the previous results, the true values of 

β1 fall outside the con®dence intervals when the EM algorithm is not used. When the adjustment is made using EM 

algorithm III, the true values of β1 are more likely to fall within the con®dence intervals. 

Figure 5: Plots of estimates with response and covariate subject to PRAM, with 95% CI, when n = 10000

3 Additional Simulations

In this section, we present additional results on performance of the proposed algorithms subject to varying probabilities 

of success for a binary response variable, varying distribution of a binary covariate, and perturbing multinomial rather 

than a binary covariate. 

3.1 Adjusting Proportions of Covariate and Regression Parameters

Another simulation study was carried out, with one binary covariate, β = (1, β)

t

, xi1 sampled from Bernoulli(π), 

exp(β

t

x)

and yi sampled from Bernoulli( ) with varying the valuesof β and π to assess their effect on EM Algorithm 1+exp(βtx)

I. The following values were used: β = (−2, −0.5, 0.5, 2), and π = (0.1, 0.2, 0.3, 0.4, 0.5). PRAM is then applied to 

x to obtain x

∗

, with the following transition matrix 

 

PX =

p

1 − p

1 − p

p

,

where we vary the value of p. 

Next, using EM Algorithm I, we obtain the MLEs of model (1). The algorithm ran for 20 steps. The mean estimates 

of β1 and the coverage probabilities of the 95% CIs are computedand reported in Table 5. We ran500 simulations. 

The results are displayed in Table 5 for p = 0.9 and p = 0.8. In terms of bias and coverage probabilities, the algorithm 

gives better estimates when the value of π goes closer to 0.5, when β = −0.5, 0.5, as well as when p = 0.9. For 

example, when p = 0.9, we get high probability coverages (greater than 0.95) when π = 0.4, 0.5 and β = −0.5, 0.5. 

This suggests that the EM algorithm works better when the distribution of the response variable and covariate is not 

skewed. 

10 
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Table 5: Simulated ML estimates of (1) with accounting for PRAM, coverage probabilities in parentheses. 

β

π p -2 -0.5 0.5 2 

0.1 0.9 

0.8 

-2.0766 (0.760) 

-1.819038 (0.622) 

-0.5000 (0.898) 

-0.4564 (0.690) 

0.5785 (0.828) 

0.5514 (0.642) 

2.4087 (0.630) 

1.6592 (0.594) 

0.2 0.9 

0.8 

-2.0344 (0.910) 

-1.9753 (0.700) 

-0.5070 (0.982) 

-0.5128 (0.760) 

0.5283 (0.954) 

0.5483 (0.692) 

2.2101 (0.724) 

1.9494 (0.640) 

0.3 0.9 

0.8 

-2.0201 (0.948) 

-1.9908 (0.742) 

-0.51071 (0.996) 

-0.5159 (0.856) 

0.5051 (0.970) 

0.5281 (0.780) 

2.1480 (0.788) 

2.1353 (0.618) 

0.4 0.9 

0.8 

-2.0174 (0.976) 

-2.0213 (0.784) 

-0.5037 (0.990) 

-0.5050 (0.866) 

0.5071 (0.986) 

0.5213 (0.806) 

2.0640 (0.882) 

2.1020 (0.642) 

0.5 0.9 

0.8 

-2.0214 (0.982) 

-2.0110 (0.822) 

-0.5165 (0.994) 

-0.5127 (0.844) 

0.5131 (0.992) 

0.5216 (0.852) 

2.0572(0.936) 

2.1020 (0.716) 

3.2 Covariate with More than Two Levels

Another simulation study was carried out, with one categorical covariate with 3 levels, β = (1, −1, 1)t, xi1 sampled 

exp(β

t

x)

from multinomial(0.4, 0.3, 0.3), and yi sampled from Bernoulli( ). The estimates of β1 and β2 from the 1+exp(βtx)

∗

logistic regression of Y on X are reportedin Table6, labeledasaβi,real”. PRAM is then applied to y to obtain y , 

with the following transition matrix 

 

PY =

p

1 − p

1 − p

p

,

where we vary the value of p. 

Standard logistic regression is performed on y

∗ with x, while y is not released and is considered to be unobserved. 

We estimate β1 and β2and their approximate 95% con®dence interval, βˆ

i ± 2SE(βˆi). We ran 500 simulations, with 

n = 100, 1000, 10000, p = 0.8, 0.9. The mean estimates of β1 and β2 and the coverage probabilities of the 95% CIs 

are computed and are reported in Table 6 ina row labeled aβi,noadjust”. 

Next, using EM Algorithm II, we obtain the MLEs of model (1). The algorithm ran for 20 steps. The mean 

estimates of β1 and β2 and the coverage probabilities of the 95% CIs are computedand reported in Table 6 in a row 

labeled aβi,adjust”. Similar to the previous results, the estimates from the algorithm are less biased as sample size 

increases, for both β1 and β2. For example, for p = 0.9, the bias for β2 is −0.1909, −0.0196, −0.0035 for n = 100, 

n = 1000 and n = 10000 respectively. 

Table 6: Simulated ML estimates of β1 & β2 for(1) with accounting for PRAM. Average ML estimates. Coverage 

probabilities in parentheses. 

n = 100

β1,real

β1,noadjust

β1,adjust

-1.0378 

-0.7762 (0.974) 

-1.0419 (0.972) 

β2,real

β2,noadjust

β2,adjust

1.5883 

0.7781 (0.894) 

1.7792 (0.850) 

n = 1000

β1,real

β1,noadjust

β1,adjust

-1.0043 

-0.7740 (0.800) 

-0.9988 (0.986) 

β2,real

β2,noadjust

β2,adjust

1.0137 

0.6513 (0.502) 

1.0333 (0.944) 

n = 10000

β1,real

β1,noadjust

β1,adjust

-0.9970 

-0.7739 (0.000) 

-0.9971 (0.984) 

β2,real

β2,noadjust

β2,adjust

1.0024 

0.6424 (0.000) 

1.0059 (0.952) 

Figure 6 displays the plotsof the true valuesof β1, along with the estimates of β1 and the con®dence intervals for 

the estimates via the EM algorithm for n = 10000, p = 0.90. Similar to the previous results, the algorithm works well 

since the true values of β1 are likely to fall inthe con®dence intervals forthe estimates from theEM algorithm. 

11 
 

 

 

 

 

 

 

 

 

 

 

Without Adjustment, p=0.9 

"' 9 

00 

'i' 

N 

'";-

10 20 30 

Index 

40 

real 

no adjust 

Cl for no adjust 

50 

With Adjustment, p=0.9 

"' 9 

real 

EM estinate 

Cl w~h EM estmate 

10 20 30 40 50 

Index 

Figure 6: Plots of estimates with covariate subject to PRAM, with 95% CI, when n = 10000, p = 0.90

In the next simulation study, we had one categorical covariate with 3 levels, β = (1, −1, 1)t

, xi1 sampled from 

exp(β

t

x) multinomial(0.4, 0.3, 0.3), and yi sampled from Bernoulli( ). The estimate of β1 and β2 from the logistic 1+exp(βtx)

regression of Y on X are reportedin Table7, labeledasaβi,real”. PRAM is then applied to x to obtain x

∗

, with the 

following transition matrix 

⎛

p

1−p

2

1−p

2

⎞

= ⎝

1−p PX 2

1−p

2

p

1−p

2

1−p

2

p

⎠ ,

where we vary the value of p. 

Standard logistic regression is performed on y with x

∗

, while x is not released and is considered to be unobserved. 

We estimate β1 and β2 and their approximate 95% con®dence interval, βˆ

i ± 2SE(βˆi).We ran 500 simulations, with 

n = 100, 1000, 10000, p = 0.8, 0.9. The mean estimates of β1 and β2 and the coverage probabilities of the 95% CIs 

are computed and are reported in Table 7 ina row labeled aβi,noadjust”. 

Next, using EM Algorithm I, we obtain the MLEs of model (1). The algorithm ran for 20 steps. The mean estimates of β1 and β2 and the coverage probabilities of the 95% CIs are computed and reportedin Table 7 ina row 

labeled aβi,adjust”. Similar to the previous results, the estimates from the algorithm are less biased as sample size 

increases, for both β1 and β2. 

Figure 7 displays the plotsof the true valuesof β1, along with the estimates of β1 and the con®dence intervals for 

the estimates via the EM algorithm for n = 10000, p = 0.90. Similar to the previous results, the algorithm appears 

to works well since the true values of β1 are likely to fall in the con®dence intervals for the estimates from the EM 

algorithm. 

4 Application to 1993 CPS Dataset

We implement the methodology describedin Section 2 on data from the 1993 Current Population Survey (CPS). The 

dataset contains 48842 records on 8 categorical variables. We perform logistic regression for salary (0 = <$50,000 

or 1 = >$50,000) on the covariates sex (0 = Female or 1 = Male), race (0 = Non White or 1 = White), and marital 

12 
 

 

 

 

 

 

 

1n w 

l,'l 

ID 

10 

Without Adjustment, p=0.9 

20 30 

Index 

40 

real 

no adjust 

Cl for no adjust 

50 

1n 

w 

l,'l 

ID 

10 

With Adjustment, p=0.9 

20 30 

Index 

real 

EM estinate 

Cl w~h EM estmllte 

40 50 

Table 7: Simulated ML estimates of β1 & β2 for(1) with accounting for PRAM. Average ML estimates. Coverage 

probabilities in parentheses. 

n = 100

β1,real

β1,noadjust

β1,adjust

-1.0378 

-0.7762 (0.974) 

-1.0419 (0.972) 

β2,real

β2,noadjust

β2,adjust

1.5883 

0.7781 (0.894) 

1.7792 (0.850) 

n = 1000

β1,real

β1,noadjust

β1,adjust

-1.0043 

-0.7740 (0.800) 

-0.9988 (0.986) 

β2,real

β2,noadjust

β2,adjust

1.0137 

0.6513 (0.502) 

1.0333 (0.944) 

n = 10000

β1,real

β1,noadjust

β1,adjust

-1.0034 

-0.8515 (0.038) 

-1.0040 (0.996) 

β2,real

β2,noadjust

β2,adjust

1.0021 

0.7706 (0.006) 

1.0038 (0.980) 

Figure 7: Plots of estimates with covariate subject to PRAM, with 95% CI, when n = 10000, p = 0.90

(0 = Married or 1 = Unmarried). The parameter estimates from ®tting the logistic regression with the original data 

are displayedin the ®rst line of Table8, labeled as aO.D.”. 

We consider the following three cases: 1) marital subject to PRAM; 2) salary subject to PRAM; and 3) both 

marital and salary subject to PRAM. In each case, the following PRAM matrix was applied to the variables that were 

subject to PRAM 

  0.9 0.1

P = .

0.1 0.9

Ineach case, standard logisticregressionis performedonthe PRAMed data. We estimate βi and their approximate 

95% con®dence intervals, βˆ

i ± 2SE(βˆi). We ran 500 simulations for each case. The mean estimates of βi and the 

coverage probabilities of the 95% CIs are computed and are reported in Table 8 in the rows labeled awithout”. 

Next, using the EM Algorithms I, II and III, we obtain the MLEs of model (1). The algorithm ran for 20 steps. 

The mean estimates of βi and the coverage probabilities of the 95% CIs are computed and reported in Table 8 in the 

rows labeled awith”. In terms ofbias, EM Algorithm II appears to workbest, followed by EM Algorithm I and EM 

Algorithm III. EM Algorithm III being least effective is no surprise, since PRAM is applied to more variables. For 

example, the bias of βˆ

3 is −0.0604, 0.0116, −0.1303 for algorithm I, II and III respectively. 

13 
 

 

' 

"' ' 

' 

0 

0 

10 

Without Adjustment 

Estrnate w/oadjust 

Cl of est 

real 

20 

Index 

30 

C) 

j 

a 

w 

40 50 

With Adjustment 

"" ' 

"' ' 

' 

0 

0 

N 

"I 

"I "" 

EM estimate 

"I Clol est 

real 

o 10 20 30 40 50 

Index 

Table 8: Parameter Estimates from Original Data (O.D.), PRAMed Data without EM algorithm, and PRAMed data 

withEM Algorithm. AverageML estimates.Coverage probabilitiesin parentheses Case1: marital subject to PRAM; 

Case 2: salary subject to PRAM; Case 3: both marital and salary subject to PRAM 

βˆ

0 βˆ1 βˆ2 βˆ3

O.D. -0.8585 0.2855 0.3925 -2.3166 

Case 1 without 

with 

-1.4458 (0) 

-1.047738 (0.290) 

0.7398 (0) 

0.3851 (0.124) 

0.4400 (0) 

0.4249 (0.734) 

-1.6009 (0) 

-2.2562 (0.510) 

Case 2 without 

with 

-0.5475 (0) 

-0.7785 (0.592) 

0.1550 (0) 

0.2138 (0.412) 

0.2323 (0) 

0.3745 (0.946) 

-1.4539(0) 

-2.3282 (0.842) 

Case 3 without 

with 

-1.2807 (0) 

-1.2469 (0.128) 

0.7283 (0) 

0.4372 (0.098) 

0.2721 (0) 

0.444 (0.468) 

-1.0581 (0) 

-2.1863 (0.262) 

Figures8,9 and10 display the plotsof the estimateof β3 on the original data, along with the estimates of β3, 

andthe con®dence intervalsforthe estimatesviatheEM algorithm,for case1, case 2 and case 3 respectively. When 

not adjusting for PRAM, the estimates of β3 from the original logistic regression fall outside the con®dence intervals. Indeed, for case 3, the estimates for β3 fall outside the range of the plot. When the adjustment is made using the 

EM algorithm,the valueof β3 from the original logistic regression is more likely to fall within the con®dence intervals. 

Figure 8: Plots of β3 with CPS data, Case 1: marital subject to PRAM 

4.1 Disclosure Risk Assessment

We need to evaluate both disclosure risk and data utility after the application of SDC methodology. This concept is 

illustrated via a risk-utility (R-U) map (for more details, see Duncan and Pearson (1991)). 

This paper focuses on preserving data utility of microdata subject to PRAM when ®tting a logistic regression on 

the PRAMed data. We now brie ̄y discuss disclosure risk assessment and how a PRAM matrix may be chosen. Ideally, the chosen PRAM matrix is one that maximizes data utility under some predetermined levels of disclosure risk 

14 
 

10 

" '7 

'7 

<D 

'7 

M 

j 

0 "I 

;rj 

N 

"I 

"I " 

"I 

0 10 

Without Adjustment 

Estinllte wlo adjusl 

Cl of est 

re.1! 

20 

Index 

30 

Without Adjustment 

Estimllte wloadJ)sl 

Cl ol est 

real 

20 30 

Index 

40 

40 

With Adjustment 

" '7 

"' '7 

<D 

'7 

M 

i 

0 "I 

1n 

w 

"I 

"I " 

EMestinllle 

"I Cl of est 

real 

50 0 10 20 30 40 50 

Index 

With Adjustment 

" '7 

'7 

<D 

'7 

M 

i 

0 "I 

;rj 

N 

"I 

"I " 

EMe stmllle 

"I Cl of est 

real 

50 0 10 20 30 40 50 

Index 

Figure 9: Plots of β3 with CPS data, Case 2: salary subject to PRAM 

Figure 10: Plots of β3 with CPS data, Case 3: marital and salary subject to PRAM 

set by the statistical agency. 

Since the estimates from the EM algorithm converge to the maximum likelihood estimates, data utility is preserved 

for all choices of PRAM matrices. Thus, the next step is to ®nd a PRAM matrix that satis®es some level of disclosure 

control or risk set by the statistical agency. A traditional measure of disclosure risk was proposed by de Wolf and 

van Gelder (2004), which involves calculating the conditional probability that given a score k in the perturbed ®le, 

the original score was k as well, P(X = k|X∗ = k). In the context of PRAM, this conditional probability can be 

estimated by 

15 
 

 

 

 

 

 

 

 

 

 

 

 

 

It------+-1------II t------+-1------

pkkTX(k)

Rˆ

P RAM(k) = P(X = k|X∗ = k) = P , (13) 

l

plkTX(l)

where TX(l) are the frequency counts in the sample for score l. The numerator of (13) estimates the number of 

scores k in the original ®le that will remain as k in the perturbed ®le, and the denominator estimates the number of 

scores k in the original ®le that remain as k in the perturbed ®le plus the number of scores that were not k in the 

original ®le that take on the score k in the perturbed ®le. According to a traditional threshold rule, a record is safe 

when its score occurs more than a certain threshold d, so a safe record can be linked with at least d records (d to be 

determinedby statistical of®ce orrules). de Wolf and vanGelder (2004)suggest that arecord is safewhenever 

Tξ(k)

Rˆ

P RAM(k) ≤ (14) d

We use case 1, marital subject to PRAM as an example. The three-way table of counts and Rˆ

P RAM(k) for marital, 

sex, and race is displayed in Table9. Rˆ

P RAM(k) for married females (both non-white and white) is much lower than 

unmarried females (non-white and white). This is expected since the number of unmarried females is greater than 

the number of married females in this example. Thus, even though both unmarried and married females had a 10% 

probability of being misclassi®ed as married and unmarried females, respectively, the actual number of unmarried 

females misclassi®ed as married females is much higher than married females being misclassi®ed. This leads to an 

observation that a higher proportion of married females in the perturbed ®le were originally unmarried females in the 

original ®le. Depending on the disclosure rules set by the agency, such values of Rˆ

P RAM(k) may provide suf®cient 

ˆ disclosure control. For example, if the agency decides that the threshold d = 800, RP RAM(k) is less than Tξ(k)

for d

all k, satisfying (14). 

Tξ(k)

Table9: Three Way Table for Marital, Sex and Race. Rˆ

P RAM(k), in parentheses d

Married 

Non-White White 

Female 521 (0.6394, 0.65) 2288 (0.6572, 2.86) 

Male 1990 (0.9029, 2.49) 18245 (0.9400, 22.81) 

Unmarried 

Non-White White 

Female 2644 (0.9786, 3.31) 10739 (0.9769, 13.42) 

Male 1925 (0.8970, 2.41) 10490 (0.8380, 13.11) 

5 Discussion

With increased computing power and built-in statistical packages, many researchers nowadays prefer to work with 

microdata instead of aggregated data. However, releasing microdata instead of aggregated data increases the risk of 

disclosure. The goal of SDC for microdata is that given an original microdata set V, a protected microdata set V’ is 

released in its place so that disclosure risk is low and data utility is high. Most measures of data utility measure the 

distortion in the distribution in V’, when compared to the distribution in V. For examples, see Domingo-Ferrer and 

Torra (2001). In our examples, we propose comparing the coef®cient estimates when ®tting the regression model on 

the original data, on the data that has been subject to PRAM, and using the EM-type methodology on the data that has 

been subject to PRAM. In general, the estimates are biased when ®tting the logistic regression model on the data that 

has been subject to PRAM; using the proposed EM-type methodology, we can obtain estimates that are much closer 

to their true values and asymptotically unbiased. The estimates are less biased with larger sample sizes, and when the 

distribution of the response variable and covariates are less skewed. In terms of coverage probabilities, the algorithm 

appears to work well since most of the coverage probabilities in our examples are high. The next step in this research 

is to extend these methodologies to handle a wider range of analyses including other models in the class of generalized 

linear models (GLMs). 

Shlomo and Skinner (2010) claim that combining sampling with a perturbation method like PRAM offers greater 

protection than using either method on its own. PRAM itself guarantees -differential privacy (see Dwork (2006)), as 

long as the PRAM matrix does not contain zero elements. We could next evaluate the effect of using both sampling 

and PRAM on the regression coef®cient estimates using our EM-algorithm. We can compare the estimates after both 

16 
sampling and PRAM has been applied to a dataset. 

There are other SDC methodologies that can be applied to categorical variables in microdata. Data swapping is 

used by agencies like the U.S. Census Bureau, and their approach guarantees that marginals involving the matching 

variables remain the same. However, the effect on regression analysis is ambiguous(e.g., see Fienberg and McIntyre 

(2004)). Reiter (2005) carried out an empirical study using fully synthetic data with the 2000 Current Population 

Study, and found that the coverage probabilities for the logistic regression are extremely low. An interesting next step 

would be to compare performance of synthetic data methodology to our proposed EM algorithms for PRAM. 

References

P. de Wolf and I. van Gelder.An Empirical Evaluation of PRAM. Technical Report Discussion Paper 04012, Statistics 

Netherlands, Voorburg/Heerlen, September 2004. 

J. Domingo-Ferrer and V. Torra. Disclosure Control Methods and Information Loss for Microdata. In P. Doyle, 

J. Lane, J. Theeuwes, and L. Zayatz, editors, Confidentiality, Disclosure and Data Access: Theory and Practical 

Applications for Statistical Agencies, chapter 5, pages 91±110. North-Holland: Elsevier, 2001. 

G.T. Duncan and R.W. Pearson. Enhancing Access to Microdata while Protecting Con®dentiality: Prospects for the 

Future. Statistical Science, 6:219±239, 1991. 

C. Dwork. Differential Privacy. In International Colloquium on Automata, Languages and Programming (ICALP), 

pages 1±12. Springer, 2006. 

S.E. Fienberg and J. McIntyre. Data Swapping: Variations on a Theme by Dalenius and Reiss. In Privacy in Statistical 

Databases, 3050, pages 14±29. Springer Berlin Heidelberg, 2004. 

S.E. Fienberg and A.B. Slavkovic.  ́ Data Privacy and Confidentiality. International Encyclopedia of Statistical Science. 

Springer-Verlag, 2010. 

J. Gouweleeuw, P. Kooiman, L. Willenborg, and P. de Wolf. Post Randomisation for Statistical Disclosure Control: 

Theory and Implementation. Journal of Official Statistics, 14(4):463±478, 1998. 

J. Ibrahim. Incomplete Data in Generalized Linear Models. Journal of the American Statistical Association, 85(411): 

765±769, 1990. 

J. Ibrahim, M.H. Chen, S.R Lipsitz, and A.H. Herring. Missing-Data Methods for Generalized Linear Models: A 

Comparitive Review. Journal of the American Statistical Association, 100(469):332±346, 2005. 

T.A. Louis. Finding the Observed Information Matrix when Using the EM Algorithm. Journal of the Royal Statistical 

Society, 44(2):226±233, 1982. 

A. Ramanayake and L. Zayatz. Balancing Disclosure Risk with Data Quality. Statistical Research Division Research 

Report Series 2010-04, U.S. Census Bureau, 2010. 

J.P. Reiter. Releasing Multiply-Imputed, Synthetic Public Use Microdata: An Illustration and Empirical Study. Journal 

of the Royal Statistical Society, 168:185±205, 2005. 

N. Shlomo and C.J. Skinner. Assessing the Protection Provided by Misclassi®cation-Based Disclosure Limitation 

Methods for Survey Microdata. Annals of Applied Statistics, 4(3):1291±1310, 2010. 

A.B.SlavkovicandJ.Lee. SyntheticTwo-Way Contingency Tablesthat Preserve Conditional Frequencies.  ́ Statistical 

Methodology, 7:225±239, 2010. 

A. vanden Houtand P. Kooiman. Estimatingthe Linear Regression Model with Categorical Covariates Subjectto 

Randomized Response. Computational Statistics & Data Analysis, 50:3311±3323, 2006. 

A. van den Hout and P. van der Heijden. Randomized Response, Statistical Disclosure Control and Misclassi®cation: 

A Review. International Statistical Review, 70(2):269±288, 2002. 

17 
P. van der Laan. The 2001 Census in the Netherlands: Integration of Registers and Surveys. In The Census of 

Population: 2000 and Beyond, Manchester, UK, 2000. Cathie Marsh Centre for Census and Survey Research, 

Faculty of Economics and Social Studies, University of Manchester. 

L. Willenborg and T. de Waal. Statistical Disclosure Control in Practice. Springer, New York, NY, 1996. 

18 
 

 

 

Without Adjustment, p=0.9 

10 20 30 40 

Incle)( 

Without Adjustment, p=0.8 

" 

10 20 30 40 

Index 

noadjJSI 

a tor no adjJSI 

50 

noadjJSI 

a tor noadjJSI 

50 

10 

10 

With Adjustment, p=0.9 

20 30 40 50 

Incle)( 

With Adjustment, p=0.8 

20 30 40 50 

Index 

Appendices

A More Figures from Section 2.1.2

Figure 11: Plots of estimates with covariate subject to PRAM, with 95% CI, when n = 1000, p = 0.90

Figure 12: Plots of estimates with covariate subject to PRAM, with 95% CI, when n = 1000, p = 0.80

19 
 

 

Without Adjustment, p=0.9 With Adjustment, p=0.9 

10 20 30 40 50 10 20 30 40 50 

Index Index 

Without Adjustment, p=0.8 With Adjustment, p=0.8 

10 20 30 40 50 10 20 30 40 50 

Index Index 

Figure 13: Plots of estimates with covariate subject to PRAM, with 95% CI, when n = 100, p = 0.90

Figure 14: Plots of estimates with covariate subject to PRAM, with 95% CI, when n = 100, p = 0.80

20 
 

 

 

0 

V 

0 

0 

0 

Without Adjustment, p=0.9 

llA .. A&,t~ 

l~JYY0V~ \~t·Y!~+f 

~ ~ \)ii \jV ' j J 

I 

10 20 30 

Index 

Without Adjustment, p=0.8 

I 

10 20 30 

Index 

40 

40 

50 

,., 

noltdjJSI 

Cl lorooad;Jsl 

50 

With Adjustment, p=0.9 

10 20 30 40 50 

Index 

With Adjustment, p=0.8 

10 20 30 40 50 

Index 

B More Figures from Section 2.2.2

Figure 15: Plots of estimates with response subject to PRAM, with 95% CI, when n = 1000, p = 0.90

Figure 16: Plots of estimates with response subject to PRAM, with 95% CI, when n = 1000, p = 0.80

21 
 

 

Without Adjustment, p=0.9 With Adjustment, p=0.9 

10 20 30 40 50 10 20 30 40 50 

Index Index 

Without Adjustment, p=0.8 With Adjustment, p=0.8 

10 20 30 40 50 10 20 30 40 50 

Index Index 

Figure 17: Plots of estimates with response subject to PRAM, with 95% CI, when n = 100, p = 0.90

Figure 18: Plots of estimates with response subject to PRAM, with CI, when n = 100, p = 0.80

22 
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

C Derivation of Weights for EM Algorithms I, II, III

C.1 Weights for EM Algorithm I

P(Y, Xnew = xj , X∗ = xk|φ

(v)

)

P(Xnew = xj |X∗ = xk, Y, φ

(v)

) =

P(Y, X∗ = xk|φ(v)))

P(Y |Xnew = xj , X∗ = xk, φ(v)

)P(Xnew = xj , X∗ = xk, φ(v))

= PJ

P(Y, X∗ = xk, X = xl|φ(v))) l

P(Y |Xnew = xj , φ(v))P(X∗ = xk|Xnew = xj )P(Xnew = xj )

= PJ

P(Y |Xnew = xl, X∗ = xk, φ(v))P(X∗ = xk|Xnew = xl)P(Xnew = xl)

l

(v)

pXjkP(Y |xj , β(v))πj = PJ (v)

l=1 pXlkP(Y |xl

, β(v))πl

C.2 Weights for EM Algorithm II

P(Ynew = j, Y

∗ = k|X, β(v)

)

P(Ynew = j|Y

∗ = k, X, β(v)

) =

P(Y ∗ = k|X, β

(v))

P(Ynew = j|X, β

(v)

)P(Y

∗ = k|Ynew = j)

= PJ

P(Y ∗ = k, Y = l|X, β

(v))

l

pY jkP(Y = j|X, β

(v)

)

= P P(Y ∗ = k|Y = l)P(Y = l|X, β

(v))

l

pY jkP(Y = j|X, β

(v)

)

= P

l

pY lkP(Y = l|X, β

(v))

C.3 Weights for EM Algorithm III

P(Y = m, X = l|Y

∗ = k, X∗ = j, β) = P(Y = m|X = l, Y∗ = k, X∗ = j, β)P(X = l|Y∗ = k, X∗ = j, β)

The ®rst part is 

P(Y = m, X = l, Y

∗ = k, X∗ = j|β)

P(Y = m|X = l, Y

∗ = k, X∗ = j, β) =

P(X = l, Y ∗ = k, X∗ = j|β)

P(Y

∗ = k|Y = m)P(Y = m|X = l, β)P(X∗ = j|X = l)P(X = l)

= P P(Y = a, X = l, Y ∗ = k, X∗ = j|β) a

PY mkP(Y = m|X = l, β)P(X∗ = j|X = l)P(X = l)

= P P(Y ∗ = k|Y = a)P(Y = a|X = l, β)P(X∗ = j|X = l)P(X = l) a

PY mkP(Y = m|X = l, β)

= P .

PY akP(Y = a|X = l, β) a

The second part is 

P(X = l, Y

∗ = k, X∗ = j|β)

P(X = l|Y

∗ = k, X∗ = j, β) =

P(Y ∗ = k, X∗ = j|β)

P P(Y = b, X = l, Y

∗ = k, X∗ = j|β)

b = P P P(Y = d, X = c, Y ∗ = k, X∗ = j|β)

c d

P P(Y

∗ = k|Y = b)P(Y = b|X = l, β)P(X∗ = j|X = l)P(X = l)

b = P P P(Y ∗ = k|Y = d)P(Y = d|X = c, β)P(X∗ = j|X = c)P(X = c)

c d

P pXljπ(l)

b

pY bkP(Y = b|X = l, β)

= P P .

c

pXcjπ(c) dpY dkP(Y = d|X = c, β)

23 
 

 

 

Hence 

P pY mkP(Y = m|X = l, β

(v)

) pXljπ(l)

b

pY bkP(Y = b|X = l, β

(v)

)

P(Y = m, X = l|Y

∗ = k, X∗ = j, β) = P P P .

pY akP(Y = a|X = l, β

(v)) pXcjπ(c) d

pY dkP(Y = d|X = c, β

(v)) a c

24