R/Abso_diff_est.R

In quantoptr: Algorithms for Quantile- And Mean-Optimal Treatment Regimes

#' @title Estimate the Gini's mean difference/mean absolute difference(MAD) for a Given Treatment Regime
#' @description Estimate the MAD if the entire population follows a 
#' treatment regime indexed by the given parameters.
#' This function supports the \code{\link{IPWE_MADopt}} function.
#' 
#' @param x    a matrix of observed covariates from the sample. 
#' Notice that we assumed the class of treatment regimes is linear.
#' This is important that columns in \code{x} matches with \code{beta}.
#' @param beta a vector indexing the treatment regime.
#'             It indexes a linear treatment regime:  
#'             \deqn{ d(x)= I\{\beta_0 + \beta_1 x_1 + ... + \beta_k x_k > 0\}.
#'             }{ d(x)= I{\beta_0 + \beta_1*x_1 + ... + \beta_k*x_k > 0}.}
#' @param a    a vector of 0s and 1s, the observed treatments from a sample
#' @param y    a vector, the observed responses from a sample
#' @param prob a vector, the propensity scores of getting treatment 1 in the samples
#' @param Cnobs A matrix with two columns, enumerating all possible combinations of
#' pairs of indexes. This can be generated by \code{combn(1:n, 2)},
#' where \code{n} is the number of unique observations.
#'
#' @export
#' 
#' @seealso The function \code{\link{IPWE_MADopt}} is based on this function.
#' 
#' @references 
#' \insertRef{wang2017quantile}{quantoptr}
#' 
#' @examples
#' library(stats)
#' GenerateData.MAD <- function(n)
#' {
#'   x1 <- runif(n)
#'   x2 <- runif(n)
#'   tp <- exp(-1+1*(x1+x2))/(1+exp(-1+1*(x1+x2)))
#'   a<-rbinom(n = n, size = 1, prob=tp)
#'   error <- rnorm(length(x1))
#'   y <- (1 + a*0.3*(-1+x1+x2<0) +  a*-0.3*(-1+x1+x2>0)) * error
#'   return(data.frame(x1=x1,x2=x2,a=a,y=y))
#' }
#' \dontshow{
#' n <- 50
#' testData <- GenerateData.MAD(n)
#' logistic.model.tx <- stats::glm(formula = a~x1+x2, data = testData, family=binomial)
#' ph <- as.vector(logistic.model.tx$fit)
#' Cnobs <- combn(1:n, 2)
#' abso_diff_est(beta=c(1,2,-1), 
#'               x=model.matrix(a~x1+x2, testData),
#'               y=testData$y,
#'               a=testData$a,
#'               prob=ph,
#'               Cnobs = Cnobs)
#' 
#' }
#' 
#' \donttest{
#' n <- 500
#' testData <- GenerateData.MAD(n)
#' logistic.model.tx <- glm(formula = a~x1+x2, data = testData, family=binomial)
#' ph <- as.vector(logistic.model.tx$fit)
#' Cnobs <- combn(1:n, 2)
#' abso_diff_est(beta=c(1,2,-1), 
#'               x=model.matrix(a~x1+x2, testData),
#'               y=testData$y,
#'               a=testData$a,
#'               prob=ph,
#'               Cnobs = Cnobs)
#' }

abso_diff_est <- function(beta,x,y,a,prob, Cnobs ){
  g <- as.numeric(x%*%beta>0)
  c <- a*g+(1-a)*(1-g)
  wts <- g*1/prob+(1-g)*(1/(1-prob))

  wtsprod <- wts[Cnobs[1,]] * wts[Cnobs[2,]]
  cprod <- c[Cnobs[1,]] * c[Cnobs[2,]]
  absd <- abs(y[Cnobs[1,]] - y[Cnobs[2,]])
  S <- sum(cprod * wtsprod * absd)
  n <- length(y)
  val <- S/(n*(n-1)/2)
  return(val)
}