R/estimands.R

In bsaul/interferenceSim: interference Simulations

#-----------------------------------------------------------------------------#
#' Individual level causal estimand Y(a, alpha)
#' 
#' Takes n X n  x 2 array as input
#' @return n x 2 matrix with values of ybar_ij(a, alpha) 
#-----------------------------------------------------------------------------#


ybar_i_aalpha <- function(outcomes, alpha){
  
  ## Necessary bits ##
  nn <- dim(outcomes)[1]
  kk <- 0:(nn - 1)
  aa <- alpha
  
  ## Create array for output ##
  out <- array(dim = c(nn, 2), dimnames = list(id = dimnames(outcomes)$id,
                                              trt = dimnames(outcomes)$trt))
  
  ybar <- function(ind, alpha){
    sum(ind * choose(nn - 1, kk) * aa^(kk) * (1 - aa)^(nn - kk - 1))
  }
  out[ , 1] <- apply(outcomes[ , , 1], 1, ybar)
  out[ , 2] <- apply(outcomes[ , , 2], 1, ybar)
  
  return(out)
}

#-----------------------------------------------------------------------------#
#' Individual level marginal causal estimand Y(alpha)
#' Takes the n x 2 matrix returned by ybar_ij_aalpha
#' 
#' @return n x 1 vector of ybar_ij(alpha) 
#-----------------------------------------------------------------------------#
ybar_i_alpha <- function(input, alpha){
  alpha*input[, 2] + (1 - alpha) * input[, 1]
}

#-----------------------------------------------------------------------------#
#' Calculate estimands
#' 
#' Takes a set of potential outcomes generated by \code{\link{generate_po}} and 
#' returns a data frame containing the population level causal estimands
#' 
#' @param potential_outcomes a list of potential outcomes generated by 
#' \code{\link{generate_po}}
#' @param alphas a vector of 'allocation schemes' in (0, 1)
#' @return a data frame with 4 columns: 1) alpha (the allocation scheme) 2)
#' a0 - the outcome estimand for untreated, 3) a1 - the outcome estimand for
#' treated, 4) marg - the marginal outcome estimand
#' @examples 
#' calc_estimands(generate_po(sim_base_data(3, 40)), seq(.2, .8, by = .1))
#' @export
#-----------------------------------------------------------------------------#

calc_estimands <- function(potential_outcomes, alphas){
  
  ## Calculate population level outcomes based on 
  ## individual --> group --> population
  estimands <- function(potential_outcomes, alpha){
    
    # --- Individual average potential outcomes --- #  
    ind_avg_po <- lapply(potential_outcomes, function(x) ybar_i_aalpha(x, alpha))
    
    # --- Group average potential outcomes --- #                
    grp_avg_po <- lapply(ind_avg_po, function(x) apply(x, 2, mean))
    grp_avg_po <- matrix(unlist(grp_avg_po), ncol=2, byrow = T)
    # --- Population average potential outcomes --- # 
    pop_avg_po <- apply(grp_avg_po, 2, mean)
    
    # --- Marginal potential outcomes --- #
    ind_marg_po <- lapply(ind_avg_po, function(x) ybar_i_alpha(x, alpha))
    grp_marg_po <- unlist(lapply(ind_marg_po, mean))
    pop_marg_po <- mean(grp_marg_po)
    
    return(data.frame(alpha = alpha, 
                      a0    = pop_avg_po[1],
                      a1    = pop_avg_po[2],
                      marg  = pop_marg_po))
  }
  
  ## Apply the above function to each alpha level ##
  true_effects <- lapply(alphas, estimands, 
                         potential_outcomes = potential_outcomes)
  
  truth <- do.call(rbind, true_effects)
  
  return(truth)
}