R/optimal.DN.R

In datanugget: Create, Optimize, and Refine Data Nuggets

# This function finds the optimal number of datanuggets for the original data 
# using the propensity score index and creates them using create.DN() function.
# It returns a list of four objects:
# 1. optimal data nugget number.
# 2. optimal datanugget object of create.DN() function based on optimal number.
# 3. Elbow plot.
# 4. Relative second order differences plot.

# Function Inputs:
# x: original dataset.
# center.method: the method used for choosing data nugget centers. Must be 'mean' or 'random' or 'original'.
# dn.nos: vector of candidate datanugget numbers.
# R: the number of observations to sample from the data matrix when creating the initial data nugget centers.
# delete.percent: proportion of data points to be deleted at each iteration.
# DN.num1: the number of initial data nugget centers to create.
# eps: tolerance for stoppage when hit the elbow.
# dist.metric: pairwise distance measure (e.g., "euclidean" or "manhattan").
# seed: random seed for replication.
# no.cores: number of cores used for parallel processing.
# make.pbs: logical; whether to show a progress bar while the function runs.


optimal.DN <- function(x,
                       center.method = "mean", 
                       dn.nos, 
                       R = 5000,
                       delete.percent = .1,
                       DN.num1 = 10^4,
                       eps = 5e-3, 
                       dist.metric = "euclidean", 
                       seed = 291102,
                       no.cores = (parallel::detectCores() - 1),
                       make.pbs = FALSE){
  
  
  
  ## -------------------------------
  ## Argument checks
  ## -------------------------------
  
  # check x
  if (!any(class(x) %in% c("matrix", "data.frame", "data.table"))){
    stop('x must be of class "matrix", "data.frame", or "data.table"')
  }
  
  
  
  # check center.method
  if (!(center.method %in% c("mean", "random", "original"))){
    stop('center.method must be "mean" or "random" or "original"')
  }
  
  
  
  # check dn.nos
  if (!is.vector(dn.nos) || !(is.integer(dn.nos) || is.numeric(dn.nos))) {
    stop("dn.nos must be a vector of numeric or integer entries")
  }
  
  
  
  # make sure dn.nos must be >=3
  if (!(length(dn.nos) >= 3)) {
    stop("dn.nos must be a vector of length 3 or more")
  }
  
  
  
  # check R
  if (!(class(R) %in% c("numeric", "integer"))){
    stop('R must be of class "numeric" or "integer"')
  }
  
  
  
  # make sure R is between 100 and 10000
  if (R < 100 | R > 10000){
    stop("R must be within [100,10000]")
  }
  
  
  
  # check delete.percent 
  if (!is.numeric(delete.percent)){
    stop('delete.percent must be of class "numeric"')
  }
  
  
  
  # make sure delete.percent is between 0 and 1
  if (delete.percent <= 0 | delete.percent >= 1){
    stop("delete.percent must be within (0,1)")
  }
  
  
  
  # check DN.num1
  if (!(class(DN.num1) %in% c("numeric", "integer"))){
    stop('DN.num1 must be of class "numeric" or "integer"')
  }
  
  
  
  # make sure DN.num1 is greater than all candidate data nugget numbers
  if (!(DN.num1 >= max(dn.nos))){
    stop('DN.num1 must be larger than all dn.nos values')
  }
  
  
  
  # check eps
  if (!is.numeric(eps)){
    stop('eps must be of class "numeric"')
  }
  
  
  
  # check dist.metric 
  if (dist.metric != "euclidean" & dist.metric != "manhattan"){
    stop('dist.metric must be "euclidean" or "manhattan"')
  }
  
  
  
  # check seed 
  if (!is.numeric(seed)){
    stop('seed must be of class "numeric"')
  }
  
  
  
  # check no.cores 
  if (!(class(no.cores) %in% c("numeric", "integer"))){
    stop('no.cores must be of class "numeric" or "integer"')
  }
  
  
  
  # check make.pbs 
  if (!is.logical(make.pbs)){
    stop("make.pbs must be TRUE OR FALSE")
  }
  
  
  ## ----------------------------------
  ## Pre processing and Initialization
  ## ----------------------------------
  
  # convert the input to matrix
  x_mat <- as.matrix(x)
  
  
  # check x elements
  if (!is.numeric(x_mat)) stop("x must contain only numeric columns")
  
  
  # storage mode set to double
  storage.mode(x_mat) <- "double"
  
  
  # number of columns in the data
  n.col <- ncol(x_mat)
  
  
  # number of candidate datanugget numbers
  n.choices <- length(dn.nos)
  

  # Storing the datanuggets and their propensity score indices
  dn <- list()
  prop_score_idx <- numeric(n.choices)
  
  
  
  ## -----------------------
  ## Working function
  ## -----------------------
  
  for (i in 1:n.choices){
    
    number = dn.nos[i]
    
    cat("Processing data nugget number:", number, "\n")
    
    
    # Create data nugget
    dn[[i]] <- create.DN(x = x,
                         center.method = center.method, 
                         R = R,
                         delete.percent = delete.percent,
                         DN.num1 = DN.num1,
                         DN.num2 = number, 
                         dist.metric = dist.metric,
                         seed = seed,
                         no.cores = no.cores,
                         make.pbs = make.pbs)
    
    
    # candidate data nuggets information 
    candidate.DN.information <- dn[[i]][["Data Nuggets"]]
    candidate.DN.centers <- as.matrix(
      candidate.DN.information[, 1 + seq_len(n.col), drop = FALSE])
    candidate.DN.weight <- as.integer(candidate.DN.information[, "Weight"])
    
    
    # augmented dataset
    augment.df <- as.data.frame(rbind(candidate.DN.centers, x_mat))
    
    
    
    # weights of the augmented dataset
    augment.df$w <- 1
    augment.df$w[1:nrow(candidate.DN.centers)] <- candidate.DN.weight
    
    
    # Indicator variable
    augment.df$DN_flag <- 0
    augment.df$DN_flag[1:nrow(candidate.DN.centers)] <- 1
    
    
    
    # weighted propensity score model
    vars <- setdiff(names(augment.df), c("DN_flag", "w"))
    ps_formula <- as.formula(paste("DN_flag ~",
                                   paste0("s(", vars, ")", collapse = " + ")))
    suppressWarnings(
      ps_model <- mgcv::bam(ps_formula, family = binomial, 
                            data = augment.df, 
                            discrete = TRUE, weights = augment.df$w)
    )
    
    
    
    # Predicted propensity scores
    augment.df$prop_score <- predict(ps_model, type = "response")
    
    
    
    # weighted mean and variances of the estimated propensity scores
    tmp.w  <- augment.df$w / sum(augment.df$w)
    mean_dn <- sum(tmp.w * augment.df$prop_score)
    var_dn  <- sum(tmp.w * (augment.df$prop_score - mean_dn)^2)
    
    
    # propensity score indices
    prop_score_idx[i] = var_dn + (mean_dn - 0.5)^2
    
    cat("Propensity Score Index:", prop_score_idx[i], "\n")
    
    
    # stop early if run for at least 3 candidates and difference in propensity score indices < eps
    if ((i >= 3) && 
        (prop_score_idx[i] < eps * 10) && 
        (abs(prop_score_idx[i - 1] - prop_score_idx[i]) < eps)) {
      
      cat("No significant reduction in propensity score index. Stopping early.\n", sep = "")
      
      # shrink vectors to actual computed values
      dn.nos <- dn.nos[1:i]
      prop_score_idx <- prop_score_idx[1:i]
      dn <- dn[1:i]
      break
      
    }
  }

  ## ----------------
  ## Elbow Plot
  ## ----------------
  
  # data frame of the datanugget number candidates and the propensity score indices.
  ps_vs_dn <- data.frame(dn.nos, prop_score_idx)
  
  # Elbow Plot 
  elbow.plot <- ggplot2::ggplot(ps_vs_dn, aes(x = dn.nos, y = prop_score_idx)) +
    geom_point() +                 
    geom_line() +                   
    ggtitle("Elbow Plot") +                
    xlab("Data Nugget number") +                   
    ylab("Propensity Score Index") +                   
    theme_minimal()
  
  
  ## ----------------------------------------
  ## Relative Second Order differences plot
  ## ----------------------------------------
  
  
  # Optimal data nugget number based on the max relative second derivative
  dff <- diff(prop_score_idx)
  delta <- 1 - (dff[-1]/dff[-length(dff)])
  deltas <- c(NA, NA, delta)
  delta_vs_dn <- data.frame(dn.nos, deltas = deltas)
  
  
  # Difference plot
  diff.plot <- ggplot2::ggplot(delta_vs_dn, aes(x = dn.nos, y = deltas)) +
    geom_point() +                 
    geom_line() +                   
    ggtitle("Difference Plot") +                
    xlab("Data Nugget number") +                   
    ylab("Relative second order difference") +                   
    theme_minimal()   
  
  
  # Optimal datanugget number and create.DN output
  opt.idx <- which.max(deltas)
  opt.dn.no <- dn.nos[opt.idx]
  opt.dn <- dn[[opt.idx]]
  
  cat("Optimal Number of datanuggets:", opt.dn.no, "\n")
  
  return(
    list(
      opt.dn.no = opt.dn.no,
      opt.dn = opt.dn,
      elbow.plot = elbow.plot, 
      diff.plot = diff.plot)
  )
}