#------------------------------------------------------------------------------#
# component of sandwich estimator for IPW estimator of value of a fixed regime
# with K > 1
#
# ASSUMPTIONS
# all K propensity score models are logistic models
# txs are denoted by Ak (k=1,...,K) and are binary coded as {0,1}
#
# INPUTS
# moPS : list of modeling objects for propensity score regressions
# must be logistic models; kth element corresponds to kth decision
# point
# data : data.frame containing all covariates and txs
# txs must be named 'Ak' (k=1,...,K) and coded 0/1
# y : outcome of interest
# regime : matrix of 0/1 containing recommended txs
# kth column corresponds to kth decision point
#
# RETURNS
# a list containing
# value : list returned by value_IPW()
# MM : M M^T matrix
# dMdG : matrix of derivatives of M wrt gamma
#------------------------------------------------------------------------------#
value_IPW_swv_md <- function(moPS, data, y, regime) {
# estimated value
value <- value_IPW_md(moPS = moPS, data = data, y = y, regime = regime)
# number of decision points
K <- length(x = moPS)
# number of patients in data set
n <- nrow(x = data)
# tx names in data set
nms <- paste0('A',1L:K)
piMM <- list()
pid <- rep(x = 1.0, times = n)
for (k in 1L:K) {
# pi(x_k; gamma_k)
pk <- modelObj::predict(object = value$fitPS[[ k ]], newdata = data)
if (is.matrix(x = pk)) pk <- pk[,ncol(x = pk), drop = TRUE]
# propensity for receiving recommended tx
pid <- pid * {pk*{regime[,k] == 1L} + {1.0 - pk}*{regime[,k] == 0L}}
# derivative
piMM[[ k ]] <- stats::model.matrix(object = modelObj::model(object = moPS[[ k ]]),
data = data)*{-regime[,k] + pk}
}
# indicator of tx received was that dictated by the regime
Cd <- as.integer(x = {rowSums(x = regime == data[,nms]) == K})
# value component of M MT
mmValue <- mean(x = {Cd * y / pid - value$valueHat}^2)
#### derivative w.r.t. gamma
dMdG <- NULL
for (k in 1L:K) {
dMdG <- c(dMdG, colMeans(x = Cd * y / pid * piMM[[ k ]]))
}
return( list("value" = value,
"MM" = mmValue,
"dMdG" = dMdG) )
}