Penalized Gaussian copula with empirical marginals for data synthetization in R

glasso_emp_copula.R

# Graphical LASSO copula model
library(glasso)
library(MASS)
library(Matrix)
set.seed(45)

glasso_copula <- function(df, rho, penalize.diagonal = FALSE, ...) {
  N <- nrow(df)
  P <- ncol(df)
  if (missing(rho)) rho <- 2*sqrt(log(P)/N) # sls ch 9 p 252
  Z <- as.data.frame(lapply(df, function(x) {
    r <- rank(x, na.last = "keep", ties.method = "average")
    qnorm(r / (sum(!is.na(r)) + 1))
  }))
  Shat    <- glasso(cov(Z), rho = rho, penalize.diagonal = penalize.diagonal, ...)
  Shat$w  <- as(Shat$w, "dsCMatrix")
  structure(list(Shat = Shat, df = df, rho = rho), class = "glcopula")
}

simulate.glcopula <- function(object, nsim = 1) {
  N <- nrow(object$df)
  P <- ncol(object$df)
  Z <- mvrnorm(nsim, mu = rep(0, P), Sigma = object$Shat$w)
  X <- lapply(1:P, function(p) {
    if (!is.factor(object$df[[p]])) {
      out <- quantile(object$df[[p]], pnorm(Z[,p]), type = 1, names = FALSE)
    } else {
      qq <- quantile(as.integer(object$df[[p]]), pnorm(Z[,p]), type = 1, names = FALSE)
      ll <- levels(object$df[[p]])
      out <- factor(ll[qq], levels = ll)
    }
    out
  })
  as.data.frame(X, col.names = colnames(object$df))
}

print.glcopula <- function(x, ...) {
  cat(" Graphical LASSO copula with empirical marginals\n",
      "ϱ:", x$rho, "\n")
  if (ncol(x$Shat$w) > 8) {
    cat(" Σ[1:7, 1:7]:")
    m <- x$Shat$w[1:7, 1:7]
  } else {
    cat(" Σ:")
    m <- x$Shat$w
  }
  printSpMatrix(round(m, 3))
}

# example usage
fit <- glasso_copula(iris, rho = 0)
fit
new_iris <- simulate(fit, 150)


# penalization & compare to synthpop
fit_plo <- glasso_copula(iris, rho = .4)
fit_plo

fit_phi <- glasso_copula(iris, rho = .6)
fit_phi

plo_iris <- simulate(fit_plo, 150)
phi_iris <- simulate(fit_phi, 150)
syn_iris <- synthpop::syn(iris)$syn

# worse prediction as we move to higher penalization
# as it should be!
data.frame(
  true     = iris$Species,
  synthpop = predict(lda(Species~., syn_iris), iris)$class,
  copula   = predict(lda(Species~., new_iris), iris)$class,
  pen_lo   = predict(lda(Species~., plo_iris), iris)$class,
  pen_hi   = predict(lda(Species~., phi_iris), iris)$class
)