Statistical Methods
Univariate Analysis
The original RR design proposed by Warner in 1965 works as follows:
With randomization probability p, participants are supposed
to answer to the question ‘Have you ever used cocaine?’, and with
counter-probability 1-p, the question is reversed, i.e.,
‘Have you never used cocaine?’. As a randomization device, dices or coin
flips with known probabilities may be used. The model is depicted in the
following figure:
First, let’s load the package and generate a
single data frame using the function RRgen. We simulate
data for 1000 participants, a true prevalence of 30% cocaine users
(pi.true = .30; obviously, this rate is unrealistic, but it
should suffice for the sake of example), and a randomization probability
of p=1/6 (e.g., resulting from a dice roll):
library(RRreg)
set.seed(123)
data.W <- RRgen(
n = 1000,
pi.true = .30,
model = "Warner",
p = 1 / 6,
complyRates = c(1, 1),
sysBias = c(0, 0)
)
head(data.W)## true comply response
## 1 0 1 1
## 2 0 1 1
## 3 0 1 1
## 4 1 1 0
## 5 1 1 0
## 6 0 1 0In the function RRgen, the argument
complyRates gives the proportions of participants who do
adhere to the RR instructions for carriers and non-carriers of the
sensitive attribute, respectively. The argument sysBias
gives the probability of yes responses in case of
non-compliance for carriers and non-carriers, respectively. If
sysBias = c(0, 0), noncomplying individuals always respond
with self-protective no responses (also known as SP-‘no’
responses). In comparison, sysBias = c(0, 0.5) indicates
that noncomplying non-carriers answer yes or no by
chance. However, since we generated data assuming perfect compliance,
this argument can be neglected. These options can be used to test how
robust the results are against non-adherence.
Now that we have a data set, we can estimate the prevalence of
cocaine users by means of RRuni:
warner <- RRuni(
response = response,
data = data.W,
model = "Warner",
p = 1 / 6,
MLest = TRUE
)
summary(warner)## Warner Model with p = 0.166666666666667
## Sample size: 1000
##
## Estimate StdErr z Pr(>|z|)
## pi 0.300500 0.022874 13.137 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The response variable containing 1 for yes and
0 for no responses can either be given as a vector or as a
column name in the data set data. Moreover, the argument
model specifies the RR design with its randomization
probability given by p (see Sections 4 and 5 for details
about other RR models) . If MLest = FALSE is used, the
prevalence estimate is derived as a closed-form least-squares estimator,
which can assume values outside of the interval [0,1]. However, using
MLest = TRUE (default), estimates of pi are
restricted to the interval [0,1].
- Warner, S. L. (1965). Randomized response: A survey technique for eliminating evasive answer bias. Journal of the American Statistical Association, 60, 63-69.
Multivariate Analysis
First, we generate a continuous, non-RR variable. According to our simulation, respondents carrying the sensitive attribute (i.e., cocaine users) have higher scores on this covariate:
data.W$cov[data.W$true == 1] <- rnorm(sum(data.W$true == 1), 1)
data.W$cov[data.W$true == 0] <- rnorm(sum(data.W$true == 0))
RRcor: RR correlations
Now, we can estimate the bivariate (point-biserial) correlation between the dichotomous Warner RR variable and the continuous covariate:
RRcor(
x = data.W$response,
y = data.W$cov,
models = c("Warner", "direct"),
p.list = list(1 / 6),
bs.n = 0,
bs.type = c("se.n", "se.p", "pval"),
nCPU = 1
)## Randomized response variables:
## Variable RRmodel p
## 1 data.W$response Warner 0.1667
## 2 data.W$cov direct
##
## Sample size N = 1000
##
## Estimated correlation matrix:
## data.W$response data.W$cov
## data.W$response 1.000000 0.436066
## data.W$cov 0.436066 1.000000The argument models defines which variables given by
x and y are treated as RR variables.
x and y can be vectors and/or data
frames/matrices having one variable in each column (the order is then
given as following: First, columns of x from left to right,
and second, columns of y from left to right). The
randomization probabilities p.list have to be in the same
order as models. No values have to be specified for direct
(nonRR) variables. To obtain bootstrapped standard errors,
bs.n has to be set to a number of bootstrap replications
larger than 0. Additionally, bs.type defines which kind of
bootstrap should be performed ("se.n", "se.p":
nonparametric/parametric for standard errors; "pval"
parametric bootstrap assuming independence for p-values). For faster
processing, nCPUcan be increased to use parallel
computation. In our example, we also know the true states of the
participants in the data set and can check the results by:
cor(x = data.W$true, y = data.W$cov)## [1] 0.4491291- Fox, J. A., & Tracy, P. E. (1984). Measuring associations with randomized response. Social Science Research, 13, 188-197.
RRlog: Logistic RR Regression
Alternatively, we can run a logistic regression to predict the
probability of cocaine use. In the formula below, the
Warner RR variable is defined as criterion on the left side, whereas the
continuous covariate is used as predictor on the right. The interface is
similar to that of the standard linear regression in R (i.e.,
lm()).
log1 <- RRlog(
formula = response ~ cov,
data = data.W,
model = "Warner",
p = 1 / 6,
fit.n = 1
)
summary(log1)## Call:
## RRlog.formula(formula = response ~ cov, data = data.W, model = "Warner",
## p = 1/6, fit.n = 1)
##
## RR Model:
## Warner with p = 0.167
##
## Model fit:
## n logLik
## 1000 -617.2771
##
## Estimate StdErr Wald test Pr(>Chi2,df=1) deltaG2 Pr(>deltaG2)
## (Intercept) -1.35386 0.16999 63.43126 0.00000 116.013 < 2.2e-16 ***
## cov 1.01691 0.14569 48.72026 0.00000 80.125 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1To remove the intercept from the model, use
formula=response~cov - 1. The argument LR.test
determines whether the coefficients should be tested by means of a
likelihood ratio test, that is, by removing one predictor at a time.
Again, we can check the results by:
glm(
formula = true ~ cov,
data = data.W,
family = binomial(link = "logit")
)##
## Call: glm(formula = true ~ cov, family = binomial(link = "logit"),
## data = data.W)
##
## Coefficients:
## (Intercept) cov
## -1.432 1.113
##
## Degrees of Freedom: 999 Total (i.e. Null); 998 Residual
## Null Deviance: 1217
## Residual Deviance: 993.3 AIC: 997.3- van den Hout, A., van der Heijden, P. G., & Gilchrist, R. (2007). The logistic regression model with response variables subject to randomized response. Computational Statistics & Data Analysis, 51, 6060-6069.
RRmixed: Logistic RR Regression with Mixed Effects
If the data are nested (e.g., students within schools), the logistic
regression can be extended to include random effects. In R, the package
lme4 allows for flexible generalized linear mixed models
(GLMMs) which include RR logistic regression as a special case. The
function RRmixed uses the function glmerof
lme4 and an adjusted link function. However, at the moment
this can only be used with one-group RR models that assume the same link
function for all responses.
Some typical examples how to specify models fia the
formula are:
- random effects:
rrt ~ covariate + (1 | cluster) - random slopes:
rrt ~ covariate + (0+covariate | cluster) - both:
rrt ~ covariate +(covariate | cluster) - level-2 predictors:
rrt ~ covariate + lev2pred + (1 | cluster)(the predictor lev2pred must have constant values within clusters)
The function RRmixed is used as follows:
# make random cluster:
data.W$cluster <- c("a", "b")
mixmod <- RRmixed(
formula = response ~ cov + (1 | cluster),
data = data.W,
model = "Warner",
p = 1 / 6
)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00301051 (tol = 0.002, component 1)mixmod## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logRR(0.833333333333333,0.166666666666667) )
## Formula: response ~ cov + (1 | cluster)
## Data: data.W
## AIC BIC logLik deviance df.resid
## 1240.5541 1255.2774 -617.2771 1234.5541 997
## Random effects:
## Groups Name Std.Dev.
## cluster (Intercept) 0.0005283
## Number of obs: 1000, groups: cluster, 2
## Fixed Effects:
## (Intercept) cov
## -1.354 1.017
## optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings- van den Hout, A., van der Heijden, P. G., & Gilchrist, R. (2007). The Logistic Regression Model with Response Variables Subject to Randomized Response. Computational Statistics & Data Analysis, 51, 6060–6069.
RRlin: RR Predictors in a Linear Regression
Third, we can use the Warner RR variable as predictor in a linear regression:
lin1 <- RRlin(
formula = cov ~ response,
data = data.W,
models = "Warner",
p.list = 1 / 6,
bs.n = 0,
fit.n = 1
)
summary(lin1)## Call:
## RRlin(formula = cov ~ response, data = data.W, models = "Warner",
## p.list = 1/6, bs.n = 0, fit.n = 1)
##
## Randomized response variables:
## Variable Model p
## 1 response Warner 0.167
##
## Coefficients (beta):
## Estimate StdErr Wald test Pr(>Chi2,df=1)
## (Intercept) -0.01317 0.04512 0.0852 0.7704
## response 1.02931 0.10611 94.1002 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error (sigma): 0.995; N=1000As in RRcor, bootstrapped standard errors can be
obtained by setting bs.n to the number of bootstrap
replications. The order of the RR models is defined by
models and p.list. For instance, given three
RR variables rr1, rr2, and rr3
and two non-RR/direct variables x1 and x2, the
model could be specified by
formula = y ~ rr1 + x1 + rr2 + rr3 + x2 and
models=c('Warner', 'direct', 'Warner', 'Warner', 'direct').
Note that the predictors in formula can either be vectors
or column names of variables in the data frame data.
As a special case, if the number of predictors used in the formula
interface is larger than the number of RR variables specified in
models, the remaining predictors are treated as directly
measured, non-RR variables:
data.W$pred <- rnorm(1000)
lin2 <- RRlin(
formula = cov ~ response + pred,
data = data.W,
models = "Warner",
p.list = list(1 / 6),
bs.n = 0,
fit.n = 1
)
summary(lin2)## Call:
## RRlin(formula = cov ~ response + pred, data = data.W, models = "Warner",
## p.list = list(1/6), bs.n = 0, fit.n = 1)
##
## Randomized response variables:
## Variable Model p
## 1 response Warner 0.167
##
## Coefficients (beta):
## Estimate StdErr Wald test Pr(>Chi2,df=1)
## (Intercept) -0.01205 0.04507 0.0715 0.7892
## response 1.02376 0.10640 92.5788 <2e-16 ***
## pred -0.03897 0.03353 1.3509 0.2451
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error (sigma): 0.995; N=1000To check the results of the simulation, run:
lm(cov ~ true + pred, data = data.W)##
## Call:
## lm(formula = cov ~ true + pred, data = data.W)
##
## Coefficients:
## (Intercept) true pred
## -0.02671 1.07878 -0.03004It is also possible to include several RR variables as predictors if the randomization devices used are stochastically independent:
data.W2 <- RRgen(n = 1000, pi.true = .45, model = "Warner", p = .35)
data.W$cov <-
2 * as.numeric(data.W$true) -
3 * as.numeric(data.W2$true) +
rnorm(1000, mean = 1, sd = 5)
data.W$response2 <- data.W2$response
lin3 <- RRlin(
formula = cov ~ response + response2,
data = data.W,
models = c("Warner", "Warner"),
p.list = list(1 / 6, .35),
fit.n = 1
)
summary(lin3)## Call:
## RRlin(formula = cov ~ response + response2, data = data.W, models = c("Warner",
## "Warner"), p.list = list(1/6, 0.35), fit.n = 1)
##
## Randomized response variables:
## Variable Model p
## 1 response Warner 0.167
## 2 response2 Warner 0.35
##
## Coefficients (beta):
## Estimate StdErr Wald test Pr(>Chi2,df=1)
## (Intercept) 0.78721 0.48491 2.6354 0.10450
## response 2.07495 0.83232 6.2149 0.01267 *
## response2 -3.19676 1.22070 6.8581 0.00882 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error (sigma): 4.897; N=1000Besides the regression coefficients, RRlin() also
estimates the prevalences for all combinations of RR responses:
lin3$pi## 0:0 0:1 1:0 1:1
## 0.4317846 0.2688669 0.1557977 0.1435508In our example using two Warner variables, the combination 0:0 indicates the (latent) sub-group having neither of both sensitive attributes and shows the corresponding prevalence estimate. Similarly, the combination 0:1 indicates the sub-group having only the second sensitive attribute and so on. Note that the last combination 1:1 is not provided and can be calculated by summing up all other prevalence estimates and subtracting this sum from 1.
To include many RR variables as predictors, very large sample sizes
are required for reliable estimation. In such a case, it might be
necessary to adjust the optimization parameters maxit
(maximum number of iterations) and pibeta (relative ratio
of probability scale for pi to regression coefficients
beta; e.g., choose a smaller value for larger absolute beta
values).
- van den Hout, A., & Kooiman, P. (2006). Estimating the linear regression model with categorical covariates subject to randomized response. Computational Statistics & Data Analysis, 50, 3311-3323.
