## ============================================================================
## colemanSDC_dpd2026.R
## Working with sdcMicro -- Data Protection Day 2026, University of Zurich.
## 30-minute LIVE DEMO companion (Workshop 1, KOL-G-222), repeated 4x.
##
## Author:    Matthias Templ (FHNW / SwissAnon), lead developer of sdcMicro.
## Based on:  earlier anonymization work by Roman Mueller; the Coleman demo
##            dataset re-analysed in Mueller, Templ & Strobl (2026).
## Copyright: (c) 2026 Matthias Templ. Reuse for teaching and research is
##            welcome with attribution (CC BY 4.0). Provided as-is, no warranty.
##
## Dataset: Coleman depression study (n = 107) -- 4 quasi-identifiers + 47
##          sensitive Likert items, already stripped of direct identifiers.
##          Colman et al. (2024), PsychArchives; re-analysed in Mueller,
##          Templ & Strobl (2026).
##
## STRUCTURE
##   PART A -- the 30-minute demo, five stages:
##             RISK -> RECODE -> SUPPRESS -> PERTURB -> DOCUMENT,
##             followed by a RESET block to re-run between the four sessions.
##   PART B -- EXTRAS: everything trimmed out of the live demo (SUDA,
##             sensitive-item rank swapping, the swap + PRAM alternative path,
##             fully synthetic data with synthpop, and the regression-based
##             utility comparison). Kept here for self-study and the
##             data-sharing Q&A. Run top-to-bottom -- later blocks depend on
##             earlier ones.
##
## HOW TO RUN
##   Open in RStudio and run line by line (Ctrl/Cmd+Enter). Every method below
##   can equally be done in the GUI: launch sdcApp() and use the matching
##   Risk / Recoding / Local suppression / PRAM tabs.
## ============================================================================


## 0. Make this script find its data no matter where R was started -----------
if (!file.exists("data/coleman_prepared.rda")) {
  if (requireNamespace("rstudioapi", quietly = TRUE) &&
      rstudioapi::isAvailable()) {
    setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
    message("Working directory set to: ", getwd())
  } else {
    stop("Cannot find 'data/coleman_prepared.rda'. Set the working directory ",
         "to this script's folder (Session -> Set Working Directory -> To ",
         "Source File Location) and re-run.")
  }
}

## 1. Packages -- the demo spine needs only sdcMicro -------------------------
library(sdcMicro)          # anonymization methods + the sdcApp() GUI

## 2. Load the prepared dataset ----------------------------------------------
## One load() gives: coleman (data frame), num_vars (47 Likert item names),
## scale_vars (4 scale scores), and the item-group vectors (depression_items,
## derailment_items, self_cri_items, self_re_items).
load("data/coleman_prepared.rda")
## CSV alternative:
## coleman <- read.csv("data/coleman_prepared.csv", stringsAsFactors = TRUE)

## create_scales(): recompute the four scale scores from the raw Likert items.
## Helper written by Roman Mueller (Anonymization-in-Psychology project) -- it is
## NOT part of sdcMicro, so we provide it here. A copy is also embedded in the
## .rda above; defining it explicitly keeps the script self-contained and the
## logic visible. Needs the item-group vectors loaded above. Used in Part B.
create_scales <- function(df) {
  df$Depression       <- rowSums (df[, depression_items], na.rm = TRUE)
  df$Derailment       <- rowMeans(df[, derailment_items], na.rm = TRUE)
  df$Self_Criticism   <- rowSums (df[, self_cri_items],   na.rm = TRUE)
  df$Self_Reassurance <- rowSums (df[, self_re_items],    na.rm = TRUE)
  return(df)
}


## ===========================================================================
## PART A -- THE 30-MINUTE LIVE DEMO
## ===========================================================================

## ---------------------------------------------------------------------------
## STAGE 1 -- RISK: scenario, build the object, measure exposure
## ---------------------------------------------------------------------------
## Disclosure scenario (write it down FIRST -- it tells you when to stop):
##   Release mode : Public Use File (open download, no data-use agreement).
##   Attacker     : knows the target's Age, Sex, Occupation, Ethnicity and
##                  that they took part in the study.
##   Harm         : identity disclosure => the attacker reads the depression score.
##   Utility goal : keep the original regression reproducible.

qi_vars <- c("Age", "Sex", "Occupation", "Ethnicity")   # the four quasi-identifiers

sdc <- createSdcObj(dat     = coleman,
                    keyVars = qi_vars,
                    numVars = num_vars)        # 47 Likert items (sensitive content)

print(sdc, type = "kAnon")
## Expected (n = 107): 60 records (56%) are SAMPLE-UNIQUE (violate 2-anonymity),
##                     88 (82%) violate 3-anonymity, 102 (95%) violate 5-anonymity.
## --> 56% can be re-identified by Age x Sex x Occupation x Ethnicity. THE PROBLEM.

riskyCells(sdc, maxDim = 4, threshold = 1)
## Age alone already produces many sample-uniques: a high-cardinality integer
## blows up the key space. Fix it first, in Stage 2.
## sdcApp: Risk tab -> k-Anonymity (red rows violate k).


## ---------------------------------------------------------------------------
## STAGE 2 -- RECODE: generalise before you suppress (cheaper in utility)
## ---------------------------------------------------------------------------
## Age -> 10-year bands (10-19, 20-29, ..., 70-79)
sdc <- globalRecode(sdc, column = "Age", breaks = seq(10, 80, 10))
print(sdc, type = "kAnon")          # ~60 -> ~16 sample-uniques

## Merge rare Ethnicity categories into "Other"
sdc <- groupAndRename(sdc, var = "Ethnicity",
                      before = c("Asian", "Black", "Mixed", "Undisclosed"),
                      after  = c("Other"))
print(sdc, type = "kAnon")          # small further gain

## undolast() rolls back the last step -- essential for iterative exploration.
## (Try the same merge on Occupation; if it buys nothing, undolast(sdc).)
## sdcApp: Anonymize -> Recoding -> Global recoding / Group categories.


## ---------------------------------------------------------------------------
## STAGE 3 -- SUPPRESS: close the remaining gaps to reach k-anonymity
## ---------------------------------------------------------------------------
## importance: one integer per key variable, SAME order as qi_vars.
## SMALLER number = higher priority = suppressed LAST. We protect Age the most.
sdc <- kAnon(sdc, k = 3,
             importance = c(1,   # Age        -> protect most
                            2,   # Sex
                            3,   # Occupation
                            4))  # Ethnicity  -> suppressed first if needed
print(sdc, type = "kAnon")          # all records now satisfy 3-anonymity

print(sdc, "ls")                    # suppressions per variable -- the cost
## Rule of thumb: > ~5% suppression on a variable distorts analysis. Then prefer
## a perturbative path (Stage 4) or a synthetic twin (Part B).
## sdcApp: Anonymize -> Local suppression (set k + importance, Apply).


## ---------------------------------------------------------------------------
## STAGE 4 -- PERTURB: PRAM, the no-NA alternative to suppression
## ---------------------------------------------------------------------------
## Recode + suppress protect IDENTITY but leave NAs and can still allow
## ATTRIBUTE disclosure (a 3-anonymous group sharing one depression level).
## PRAM (Post-RAndomization Method) swaps categories via an invariant
## transition matrix (marginal frequencies preserved), so even a re-identified
## record may carry a swapped value -- and it introduces no NAs.
##
## Shown on a small separate object so the categorical QIs can take the pram
## role. In sdcApp you simply open the PRAM tab on the same live session --
## no second object needed.
sdc_pram <- createSdcObj(dat      = coleman,
                         keyVars  = c("Sex", "Occupation", "Ethnicity"),
                         pramVars = c("Sex", "Occupation", "Ethnicity"))
sdc_pram <- pram(sdc_pram)
print(sdc_pram, "pram")             # original vs PRAMmed category counts
## sdcApp: Anonymize -> PRAM -> pick variables, accept/edit matrix, Apply.


## ---------------------------------------------------------------------------
## STAGE 5 -- DOCUMENT: utility note + the reports that ship WITH the data
## ---------------------------------------------------------------------------
## In this path we changed only the QIs; the 47 analysis items are untouched,
## so the depression regression is essentially unaffected -- the cost was local
## (a few suppressed QI cells). (Attribute-disclosure protection on the items
## itself = rank swap / synthetic data, see Part B.)
coleman_release <- extractManipData(sdc)
dim(coleman_release)
head(coleman_release[, qi_vars])

## Two reports -- the actual deliverable of an anonymization:
##   internal: full risk + methods + parameters + session info (custodian / DPO).
##   external: scenario + methods + how to treat NAs / recoded factors (data user).
report(sdc, filename = "ColemanSDC_internal", internal = TRUE)
report(sdc, filename = "ColemanSDC_external", internal = FALSE)


## ===========================================================================
## You can do the same in the App by point and click:
sdcApp() 
## ===========================================================================


## ===========================================================================
## PART B -- EXTRAS (NOT shown in the 30-min demo; self-study / Q&A)
## Run top-to-bottom: later blocks depend on earlier ones.
## ===========================================================================
library(synthpop)     # synthetic data generation
library(jtools)       # plot_coefs() for the utility comparison
library(ggplot2)
## source("helpers/association_heatmap.R")   # optional association-heatmap helper

## B1. SUDA -- rank records by exposure on SMALL subsets of QIs ---------------
## Run on a FRESH raw object (before any recoding). Records unique on just one
## or two QIs score highest; score > 0.1 => candidate for targeted protection.
sdc0 <- createSdcObj(coleman, keyVars = qi_vars, numVars = num_vars)
sdc0 <- suda2(sdc0)
print(sdc0@risk$suda2)

## B2. Rank-swap the 47 sensitive items on the demo object --------------------
## Guards against ATTRIBUTE disclosure within an equivalence class.
## R0 = 0.98 preserves ~98% of multivariate correlation (~8-9% of cells move).
sdc <- rankSwap(sdc, R0 = 0.98, seed = 34)
coleman_anon <- extractManipData(sdc)
coleman_anon[, num_vars] <- round(coleman_anon[, num_vars])   # restore integer Likert
coleman_anon  <- create_scales(coleman_anon)                  # recompute scale scores
data_original <- create_scales(coleman)                       # keep an original for B5

## B3. Alternative path -- swap + PRAM instead of recode + suppress -----------
## Age moves from keyVar (categorical) to numVar (numeric, rank-swapped); the
## categorical QIs are PRAMmed. Avoids NAs and factor conversion, but formal
## k-anonymity is harder to certify.
sdcSwap <- createSdcObj(dat      = coleman,
                        keyVars  = c("Sex", "Occupation", "Ethnicity"),
                        pramVars = c("Sex", "Occupation", "Ethnicity"),
                        numVars  = "Age")
sdcSwap <- rankSwap(sdcSwap, variables = "Age", seed = 32)
print(sdcSwap, "numrisk")
sdcSwap <- pram(sdcSwap)
print(sdcSwap, "pram")
data_Swap <- extractManipData(sdcSwap)
data_Swap <- rankSwap(data_Swap, variables = num_vars, R0 = 0.98, seed = 34)
data_Swap[, num_vars] <- round(data_Swap[, num_vars])
data_Swap <- create_scales(data_Swap)

## B4. Fully synthetic data with synthpop (CART) -----------------------------
## Synthesize only the 4 QIs + 4 scale scores (n = 107 is too small for the 47
## raw items). Synthesis order: QIs first, then the Colman et al. path-model order.
data_ori  <- data_original[, c(qi_vars, scale_vars)]
syn_order <- c("Age", "Sex", "Occupation", "Ethnicity",
               "Derailment", "Self_Criticism", "Self_Reassurance", "Depression")
synData   <- syn(data_ori, visit.sequence = syn_order, seed = 12345)
data_syn  <- synData$syn
head(data_syn)

## Synthetic != safe by construction -- evaluate the disclosure metrics:
risk_eval <- disclosure(synData, data_ori, keys = qi_vars, target = "Depression")
risk_eval$ident                       # UiO / UiS / UiOiS / repU
repu <- replicated.uniques(synData, data_ori, keys = qi_vars)
data_syn[repu$repU.rm, qi_vars]
multi.disclosure(synData, data_ori, keys = qi_vars, targets = scale_vars)  # DiSCO vs Dorig

## B5. Utility -- re-fit the original regressions on every candidate ----------
## A SDC method that breaks the science is worse than no release.
form1 <- Depression     ~ (Derailment + Self_Criticism) * Self_Reassurance
form2 <- Self_Criticism ~ Derailment
datasets <- list(original = data_original,
                 trad     = coleman_anon,     # recode + suppress + item swap
                 swap     = data_Swap,        # rank swap + PRAM
                 syn      = data_syn)         # fully synthetic
fit_reg  <- function(df) list(lm1 = lm(form1, data = df),
                              lm2 = lm(form2, data = df))
models   <- lapply(datasets, fit_reg)

plot_coefs(models$original$lm1, models$trad$lm1, models$swap$lm1, models$syn$lm1,
           model.names  = c("original", "trad", "swap", "synthetic"),
           legend.title = "Depression model")
plot_coefs(models$original$lm2, models$trad$lm2, models$swap$lm2, models$syn$lm2,
           model.names  = c("original", "trad", "swap", "synthetic"),
           legend.title = "Self-criticism model")
## Compare: coefficient signs, significance, point estimates (within ~10-15%),
## confidence-interval overlap, and adjusted R^2 across the four datasets.

## --- end of script ---------------------------------------------------------