Balancing Respondent Confidentiality with Variance Estimation Precision

# analyze survey data for free (http://asdfree.com) with the r language
# how to create de-identified replicate weights
# in less than ten steps

# anthony joseph damico
# ajdamico@gmail.com


# the institute for digital research and education at ucla
# hosts a very readable explanation of replicate weights,
# how they protect confidentiality, and why they're generally awesome
# http://www.ats.ucla.edu/stat/stata/library/replicate_weights.htm



# remove the # in order to run this install.packages line only once
# install.packages( c( "survey" , "sdcMicro" ) )


# load the r survey package
library(survey)

# load the sdcMicro package
library(sdcMicro)


# load some sample complex sample survey data
data(api)


# look at the first six records of the `apistrat` data.frame object
# so you have a sense of what you're working with in this example
# this particular example is student performance in california schools
# http://r-survey.r-forge.r-project.org/survey/html/api.html
head( apistrat )


###########################################
# # # # # # # # # # # # # # # # # # # # # #
# part one.  create the replicate weights #
# # # # # # # # # # # # # # # # # # # # # #
###########################################


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# step one: construct the taylor-series linearized design #
# that you use on your internal, confidential microdata   #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #


# there are many permutations of linearized designs.
# you can find well-documented examples of tsl design setups
# by searching for the text `svydesign` here:
# https://github.com/ajdamico/asdfree/search?q=svydesign&ref=cmdform


# does your tsl design have a strata argument? #

# this design does not have a `strata=` parameter
api.tsl.without.strata <-
    svydesign(
        id = ~dnum , 
        data = apistrat , 
        weights = ~pw
    )

# this design does have a `strata=` parameter
api.tsl.with.strata <- 
    svydesign(
        id = ~dnum , 
        strata = ~stype , 
        data = apistrat , 
        weights = ~pw ,
        nest = TRUE
    )


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# step two: convert this taylor-series linearized design  #
# to one of a few choices of replication-based designs    #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #


# for unstratified designs,
# use a "jackknife delete one" method
api.jk1 <- 
    as.svrepdesign( 
        api.tsl.without.strata , 
        type = "JK1" , 
        fay.rho = 0.5 ,
        mse = TRUE ,
        compress = FALSE # note: compressed replicate weights require less RAM but complicate the weight-extraction in step three
    )


# for stratified designs,
# you might first attempt to use a
# balanced repeated replication method
# with a fay's adjustment
# (this is the most common setup from the united states census bureau)
try( {

api.fay <- 
    as.svrepdesign( 
        api.tsl.with.strata , 
        type = "Fay" , 
        fay.rho = 0.5 ,
        mse = TRUE ,
        compress = FALSE # note: compressed replicate weights require less RAM but complicate the weight-extraction in step three
    )

} , silent = TRUE )
    
# however, if the sampling plan contains an
# odd number of clusters within any stratum
# you will hit this error
# Error in brrweights(design$strata[, 1], design$cluster[, 1], ..., fay.rho = fay.rho,  : 
  # Can't split with odd numbers of PSUs in a stratum


# for stratified designs with an
# odd number of clusters in any stratum,
# use a "jackknife delete n" method
api.jkn <- 
    as.svrepdesign( 
        api.tsl.with.strata , 
        type = "JKn" ,
        mse = TRUE ,
        compress = FALSE # note: compressed replicate weights require less RAM but complicate the weight-extraction in step three
    )


# now you have a matrix of replicate weights
# stored in your replication-based survey-design


# the purpose of this exercise is to produce comparable
# variance estimates without compromising confidentiality

# therefore, run a few standard error calculations
# using the original linearized designs,
# compared to the newly-created jackknife design

# calculate the mean 1999 and 2000 academic performance index scores
# using the original stratified taylor-series design..
svymean( ~ api99 + api00 , api.tsl.with.strata )
# ..and the newly-created jackknife replication-based design
svymean( ~ api99 + api00 , api.jkn )
# the standard errors for these estimates are nearly identical #

# run the same commands as above, broken down by award program eligibility
svyby( ~ api99 + api00 , ~ awards , api.tsl.with.strata , svymean )
svyby( ~ api99 + api00 , ~ awards , api.jkn , svymean )

# in each case, the replication-based design
# (that we created off of the linearized design)
# produced a comparable variance estimate



# # # # # # # # # # # # # # # # # # # # # # # # #
# step three: extract the replication weights   #
# from your newly created survey design object  #
# # # # # # # # # # # # # # # # # # # # # # # # #

# note that this example shows how to extract
# the weights from the `api.jkn` survey object
# however, this method would be identical
# for the `api.jk1` and `api.fay` objects
# displayed above as well
# despite `api.fay` throwing an error,
# because of an uneven number of clusters within strata


# look at your survey design
api.jkn

# look at the contents of your replication-based survey design
names( api.jkn )
    
# look at the first six replicate weight records within your survey design
head( api.jkn$repweights )
# note that these weights are not *combined* by default
# in other words, they still need to be multiplied by the original weight

# you can confirm this by looking at the flag that indicates:
# "have replicate weights been combined?"
api.jkn$combined.weights
# no.  they have not been combined


# therefore, these replication weights are `uncombined`
# and will need to be analyzed by the user as such
your.replicate.weights <- data.frame( unclass( api.jkn$repweights ) )



# # # # # # # # # # # # # # # # # # # # # # # #
# you've created the set of replicate-weights #

#############################################
# # # # # # # # # # # # # # # # # # # # # # #
# part two. mask the strata from evil users #
# # # # # # # # # # # # # # # # # # # # # # #
#############################################

# now prevent users from reidentifying strata #
# # # # # # # # # # # # # # # # # # # # # # # #


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# step four: understand how a malicious user might easily #
# identify clustering variables on un-obfuscated data     #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# start with the replicate weights object from the script
transposed.rw <- data.frame( t( your.replicate.weights ) )

# the first record of apistrat has dnum==401 and stype=='E'

# going back to the original microdata set,
# eleven records are in this cluster x stratum
which( apistrat$dnum == 401 & apistrat$stype == 'E' )

# without massaging your replicate weights at all,
# a malicious user could easily view the correlations
# between each record's replicate weights
# in order to determine what other records
# land in the same cluster x strata

# for example, here are all replicate-weight records that
# perfectly (after rounding) correlate with 
# the first record.
which( round( cor( transposed.rw )[ 1 , ] ) == 1 )

# same numbers!

# (wikipedia has a great definition: http://en.wikipedia.org/wiki/Obfuscation)
# if you do not obfuscate your data, a malicious user could
# identify unique clusters and strata even if only given replicate weights


# # # # # # # # # # # # # # # #
# step five: BRING THE NOISE  #
# # # # # # # # # # # # # # # #

# set a random seed.  this allows you to
# go back to this code and reproduce the
# exact same results in the future.

# the following steps include a random process
# setting a seed enforces the same random process every time
sum( utf8ToInt( "anthony is cool" ) )
# my current favorite number is 1482, so let's go with that.
set.seed( 1482 )


# figure out how much noise to add.
noisy.transposed.rw <- addNoise( transposed.rw , noise = 1 )$xm


# remember, on the original replicate weight objects,
# the correlations between the first record and
# other records within the same clusters and strata
# were perfect.
which( sapply( cor( transposed.rw )[ 1 , ] , function( z ) isTRUE( all.equal( z , 1 ) ) ) )

# run the same test on noisified weights
which( sapply( cor( noisy.transposed.rw )[ 1 , ] , function( z ) isTRUE( all.equal( z , 1 ) ) ) )
# and suddenly none of the other records are perfectly correlated


# but even moderately correlations might allow evildoers
# to identify clusters and strata
# these records have a 0.1 or higher correlation coefficient
which( cor( noisy.transposed.rw )[ 1 , ] > 0.1 )
# these records have a 0.2 or higher correlation coefficient
which( cor( noisy.transposed.rw )[ 1 , ] > 0.2 )
# whoops.  we did not add enough noise.
# records in the same cluster x strata still have
# too high of a correlation coefficient,
# relative to other records


# okay.  this will make a big difference
# on the size of the standard errors that
# your users will actually see.

# make a clutch decision:
# how much noise can you tolerate?
# hmncyt <- 1
hmncyt <- 3
# hmncyt <- 10

# you need to run this script a few times
# because the size of your standard errors
# are going to change, depending on your data set.


# essentially, choose the lowest noise value
# that does not lead you to uncomfortable levels
# of correlation within cluster/strata
# in your replicate weights columns


# crank up the random noise percentage to three
noisy.transposed.rw <- addNoise( transposed.rw , noise = hmncyt )$xm
# and suddenly..

# records 29, 73 and 158 have a correlation coefficient
# that's greater than zero point two.
which( cor( noisy.transposed.rw )[ 1 , ] > 0.2 )

# and of those, only record #29 is in the same cluster x strata
intersect(
    which( cor( noisy.transposed.rw )[ 1 , ] > 0.2 ) ,
    which( round( cor( transposed.rw )[ 1 , ] ) == 1 )
)
# that's awesome, because you have two false-positives here.
# the "73" and "158" will throw off a malicious user.


# these records have a 0.1 or higher correlation coefficient
which( cor( noisy.transposed.rw )[ 1 , ] > 0.1 )
# and that looks very good.
# because less than half of those records
# are actually in the same cluster x strata
intersect(
    which( cor( noisy.transposed.rw )[ 1 , ] > 0.1 ) ,
    which( round( cor( transposed.rw )[ 1 , ] ) == 1 )
)
# and there are lots of other records with high correlations (false positives)
# that in fact are not in the same strata.  great.  perf.  magnifique!

# this object `noisy.transposed.rw` is the set of replicate weights
# that you might now feel comfortable disclosing to your users.

# bee tee dubs

# you should un-transpose the weights rightaboutnow
noisy.rw <- t( noisy.transposed.rw )


# # # # # # # # # # # # # # # # # # # # #
# step six: check with your legal dept. #
# # # # # # # # # # # # # # # # # # # # #

# brilliant malicious users might still be able to identify certain records
# if they work really really really hard and have access to other information
# included in your survey microdata set.

# for example:
# if your technical documentation says that memphis was one of your sampled clusters, and
# if your microdata does have a state identifier, and
# if the content of your survey was about barbeque consumption

# it's possible that no amount of masking, obfuscating, massaging, noising, whathaveyouing
# will prevent malicious users from determining the geography of some of the records
# included in your public use file.  so don't mindlessly follow this example.

# consider exactly what you're disclosing, consider how someone might use it improperly.


# # # # # # # # # # # # # # # # # # # # # # #
# you've protected cluster confidentiality  #

#############################################
# # # # # # # # # # # # # # # # # # # # # # #
# part three. share these replicate weights #
# # # # # # # # # # # # # # # # # # # # # # #
#############################################

# now share these weights with your users.  #
# # # # # # # # # # # # # # # # # # # # # # #


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# step seven: remove the confidential fields from your data #
# and tack on the replicate weight columns created above.   #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #


# go back to our original data.frame object `apistrat`
# and store that entire table into `x`
x <- apistrat
# remove the columns that were deemed too confidential to release
x$dnum <- x$stype <- NULL
# in other words, delete the cluster and strata variables
# from this data.frame object so you have a "safe-for-the-public" data set


# look at the first six records
# to confirm they have been removed
head( x )

# merge on the replicate weight data.frame
# that we just created, which contains
# obfuscated-confidential information for users
y <- cbind( x , noisy.rw )

# look at the first six records
# to confirm the weights have been tacked on
head( y )
# the replicate weights are now stored
# as `X1` through `X162`

# this data.frame object `y`
# contains all of the information that
# a user needs to correctly calculate a
# variance, standard error, confidence interval

# uncomment this line to
# export `y` to a csv file
# write.csv( y , "C:/My Directory/your microdata.csv" )
# this csv file contains your full microdata set,
# except for the cluster and strata variables
# that you had determined to be confidential information


# in other words, now you've got a public-use file (puf)
# that no longer contains identifiable geographic information


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# step eight: determine the `svrepdesign` specification to  #
# match the survey object above, but without cluster info   #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #


# once users have a copy of this de-identified
# microdata file, they will need the
# replication-based complex sample design


# there are many permutations of replication-based designs.
# you can find well-documented examples of replicate weighted design setups
# by searching for the text `svrepdesign` here:
# https://github.com/ajdamico/asdfree/search?q=svrepdesign&type=Code


# for the object `y` built above,
# construct the replication-based
# "jackknife delete n" method
# complex sample survey design object
z <-
    svrepdesign(
        data = y ,
        type = "JKn" ,
        repweights = "X[1-9]+" ,
        weights = ~pw ,
        scale = 1 ,
        combined.weights = FALSE ,
        mse = TRUE
    )
# is it giving you a warning about calculating the rscales by itself?
# good.  then you're doing right by me.


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# step nine: confirm this new replication survey object   #
# matches the standard errors derived from as.svrepdesign #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# this `z` object can now be analyzed
# and the statistics but not the standard errors
# will match the `api.jkn` shown above.
svymean( ~ api99 + api00 , z )
svymean( ~ api99 + api00 , api.jkn )

svyby( ~ api99 + api00 , ~ awards , z , svymean )
svyby( ~ api99 + api00 , ~ awards , api.jkn , svymean )
# see how the standard errors have a-little-more-than-doubled?
# that's because we had to BRING THE NOISE to the replicate weights

# if you re-run this script but lower the value of `hmncyt`
# your standard errors will *decrease* and get closer to
# what they actually are when you have the confidential information.

# that's the trade-off.  re-run this entire script but set `hmncyt <- 0`
# and you'll see some almost-perfect standard errors..
# but doing that would allow malicious users to identify clusters easily

# then re-run it again and set `hmncyt <- 10` and suddenly
# no way in hell can a malicious user identify clustering information
# but your standard errors (and subsequent confidence intervals) are ginormous

# obfuscating your replicate weights is going to make it harder for users
# to detect statistically significant differences when analyzing your microdata.
# no way around that.
# just do your best to minimize the amount of obfuscation.

# dooooo it.  public use microdata are an indisputable good.