Regression Based Estimator





The regression-based approach to estimate the value of \(\Psi\)-specific fixed regime \(d = (d_{1}(h_{1}), d_{2}(h_{2}), \dots, d_{K}(h_{K}))\), \(\widehat{\mathcal{V}}_{Q}(d)\), is as follows.

Start at Decision \(K\). For \((h_{K}, a_{K}) = (\overline{x}_{k}, \overline{a}_{K})\) define

\[ Q_{K}^{d}(h_{K},a_{K}) = E(Y|H_{K} = h_{K}, A_{K} = a_{K}), \] and let

\[ V_{K}^{d}(h_{K}) = Q_{K}^{d}\left\{h_{K}, d_{K}(h_{K})\right\}, \] which is \(Q_{K}^{d}(h_{K},a_{K})\) evaluated at \(a_{K} = d_{K}(h_{K})\).

We posit a model for \(Q_{K}^{d}(h_{K},a_{K})\), \(Q_{K}^{d}(h_{K},a_{K};\beta_{K})\) and solve a suitable M-estimating equation to estimate \(\beta_{K}\), \(\widehat{\beta}_{K}\). Once we have a fitted model, \(Q_{K}^{d}(h_{K},a_{K}; \widehat{\beta}_{K})\), form pseudo outcomes

\[ \widetilde{V}_{Ki}^{d} = Q_{K}^{d}\left\{H_{Ki}, d_{K}(H_{Ki}); \widehat{\beta}_{K}\right\}, \quad i = 1, \dots, n, \] which are the predicted outcomes based on the fitted model if each individual were to have received treatment according to rule \(d_{K}\) at Decision \(K\).

Now consider Decision \(K-1\). Posit a parametric model for \(Q_{K-1}^{d}(h_{K-1},a_{K-1})\), \(Q_{K-1}^{d}(h_{K-1},a_{K-1};\beta_{K-1})\) for

\[ E\left\{V_{K}^{d}(H_{K})|H_{K-1} = h_{K-1}, A_{K-1} = a_{K-1}\right\}, \]

and solve a suitable M-estimating equation to obtain the estimator \(\widehat{\beta}_{K-1}\) for \(\beta_{K-1}\). Once we have a fitted model, \(Q_{K-1}^{d}(h_{K-1},a_{K-1}; \widehat{\beta}_{K-1})\), form pseudo outcomes

\[ \widetilde{V}_{K-1,i}^{d} = Q_{K-1}^{d}\left\{H_{K-1,i}, d_{K-1}(H_{K-1,i}); \widehat{\beta}_{K-1}\right\}, \quad i = 1, \dots, n. \]

Continue in this fashion for each \(k=K-2,\dots,1\). Finally, the estimator for the value \(\mathcal{V}(d)\) is

\[ \widehat{\mathcal{V}}(d) = n^{-1} \sum_{i=1}^{n} \widetilde{V}_{1i}^{d}. \]

This approach is the “backward induction” approach discussed in the book. An alternative approach is that of “monotone coarsening,” which we do not discuss here.



We define function value_Q_md() to calculate the value of a \(\Psi\)-specific fixed regime using the regression-based approach.

#------------------------------------------------------------------------------#
# regression based estimator for value of a fixed tx regime with feasible tx 
#   sets at each of the K>1 decision points
#
# ASSUMPTIONS
#   each tx is coded as {0,1}
#
# INPUTS
#  moQ    : list of modeling objects for Q-function regressions
#           kth element corresponds to kth decision point
#  data   : data.frame containing all covariates and txs
#  y      : outcome of interest
#  txName : vector of tx names (each tx must be coded as 0/1)
#  regime : matrix of 0/1 containing recommended txs
#           kth column corresponds to kth decision point
#
# RETURNS
#  a list containing
#      fitQ : list of value objects returned by Q-function regressions
#  valueHat : estimated value
#------------------------------------------------------------------------------#
value_Q_md <- function(moQ, data, y, txName, regime) {
  
  # number of decision points
  K <- length(x = moQ)

  #### Q-function Regressions

  qStep <- function(moQ, data, y, txName, regime) {
    # fit Q-function model
    fitObj <- modelObj::fit(object = moQ, data = data, response = y)

    # set tx received to recommended tx
    data[,txName] <- regime

    # predict Q-function for recommended tx
    qk <- modelObj::predict(object = fitObj, newdata = data)

    return (list("fit" = fitObj, "qk" = qk))
  }

  fitQ <- list()
  response <- y

  for (k in K:2L) {
    
    ### consider only individual with more than 1 tx option
    one <- data[,txName[k-1L]] == 1L

    ## kth Q-function regression

    sk <- qStep(moQ = moQ[[ k ]],
                data = data[!one,,drop = FALSE],
                y = response[!one],
                txName = txName[k],
                regime = regime[!one,k])

    fitQ[[ k ]] <- sk$fit
    response[!one] <- sk$qk

  }

  ## k=1 Q-function regression

  s1 <- qStep(moQ = moQ[[ 1L ]],
              data = data,
              y = response,
              txName = txName[1L],
              regime = regime[,1L])

  fitQ[[ 1L ]] <- s1$fit
  response <- s1$qk

  value <- mean(x = response)

  return( list("fitQ" = fitQ, "valueHat" = value) )
}


With value_Q_md() defined (see Implementation tab), it is straightforward to estimate the value. See \(Q_{k}(h_{k},a_{k};\beta_{k})\) in the sidebar for a review of the models and their basic diagnostics.


As was done in Chapter 5, we compare the estimated value under a variety of fixed regimes. Each tab below illustrates how to define the fixed regime under consideration and complete an analysis using the posited Q-function models.





The first fixed regime we consider is the static regime \(d = \{d_{1}(x_{1}) \equiv 1, d_{2}(\overline{x}_{2}) \equiv 1, d_{3}(\overline{x}_{3}) \equiv 1\}\). That is, all individuals are recommended to receive treatment \(A_{k} = 1\), \(k = 1, 2, 3\).

The estimated value of this fixed regime is obtained as follows:

q1 <- modelObj::buildModelObj(model = ~ CD4_0*A1,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')
q2 <- modelObj::buildModelObj(model = ~ CD4_0 + CD4_6*A2,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')
q3 <- modelObj::buildModelObj(model = ~ CD4_0 + CD4_6 + CD4_12*A3,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')
regime <- matrix(data = 1L, nrow = nrow(x = dataMDPF), ncol = 3L)
value_Q_md(moQ = list(q1, q2, q3), 
           data = dataMDPF,  
           y = dataMDPF$Y,  
           txName = c("A1","A2","A3"),
           regime = regime)
$fitQ
$fitQ[[1]]

Call:
lm(formula = YinternalY ~ CD4_0 + A1 + CD4_0:A1, data = data)

Coefficients:
(Intercept)        CD4_0           A1     CD4_0:A1  
 846.727405    -0.305528    10.459837     0.006195  


$fitQ[[2]]

Call:
lm(formula = YinternalY ~ CD4_0 + CD4_6 + A2 + CD4_6:A2, data = data)

Coefficients:
(Intercept)        CD4_0        CD4_6           A2     CD4_6:A2  
   923.4224       2.0372      -1.8156     -78.5422      -0.0542  


$fitQ[[3]]

Call:
lm(formula = YinternalY ~ CD4_0 + CD4_6 + CD4_12 + A3 + CD4_12:A3, 
    data = data)

Coefficients:
(Intercept)        CD4_0        CD4_6       CD4_12           A3    CD4_12:A3  
   317.3743       2.0326       0.1371      -0.4478     603.5614      -1.9814  



$valueHat
[1] 723.2833

We see that the regression based estimator of the value under \(d = \{d_{1}(x_{1}) \equiv 1, d_{2}(\overline{x}_{2}) \equiv 1, d_{3}(\overline{x}_{3}) \equiv 1\}\) is 723.28 cells/mm\(^3\).



The second fixed regime we consider is the static regime \(d = \{d_{1}(x_{1}) \equiv 0, d_{2}(\overline{x}_{2}) \equiv 0, d_{3}(\overline{x}_{3}) \equiv 0\}\). That is, all individuals are recommended to receive treatment \(A_{k} = 0\), \(k = 1, 2, 3\).

The estimated value of this fixed regime and its standard error are obtained as follows:

q1 <- modelObj::buildModelObj(model = ~ CD4_0*A1,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')
q2 <- modelObj::buildModelObj(model = ~ CD4_0 + CD4_6*A2,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')
q3 <- modelObj::buildModelObj(model = ~ CD4_0 + CD4_6 + CD4_12*A3,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')
regime <- matrix(data = 0L, nrow = nrow(x = dataMDPF), ncol = 3L)
value_Q_md(moQ = list(q1, q2, q3), 
           data = dataMDPF,  
           y = dataMDPF$Y,  
           txName = c("A1","A2","A3"),
           regime = regime)
$fitQ
$fitQ[[1]]

Call:
lm(formula = YinternalY ~ CD4_0 + A1 + CD4_0:A1, data = data)

Coefficients:
(Intercept)        CD4_0           A1     CD4_0:A1  
    318.206        1.754      538.981       -2.053  


$fitQ[[2]]

Call:
lm(formula = YinternalY ~ CD4_0 + CD4_6 + A2 + CD4_6:A2, data = data)

Coefficients:
(Intercept)        CD4_0        CD4_6           A2     CD4_6:A2  
   318.0672       1.9303      -0.1408     527.2873      -1.6443  


$fitQ[[3]]

Call:
lm(formula = YinternalY ~ CD4_0 + CD4_6 + CD4_12 + A3 + CD4_12:A3, 
    data = data)

Coefficients:
(Intercept)        CD4_0        CD4_6       CD4_12           A3    CD4_12:A3  
   317.3743       2.0326       0.1371      -0.4478     603.5614      -1.9814  



$valueHat
[1] 1102.807

We see that the regression based estimator of the value under \(d = \{d_{1}(x_{1}) \equiv 0, d_{2}(\overline{x}_{2}) \equiv 0, d_{3}(\overline{x}_{3}) \equiv 0\}\) is 1102.81 cells/mm\(^3\).



The final fixed regime we consider is a covariate dependent regime

\[ d = \left\{ \begin{array}{l} d_{1}(x_{1}) = \text{I}(\text{CD4_0} < 250 ~\text{cells}/\text{mm}^{3}) \\ d_{2}(\overline{x}_{2}) = d_{1}(x_{1}) + \{1-d_{1}(x_{1})\} ~ \text{I}(\text{CD4_6} < 360 ~\text{cells}/\text{mm}^{3}) \\ d_{3}(\overline{x}_{2}) = d_{2}(\overline{x}_{2}) + \{1-d_{2}(\overline{x}_{2})\} ~ \text{I}(\text{CD4_12} < 300 ~\text{cells}/\text{mm}^{3}) \end{array} ~~~. \right. \]

The estimated value of this fixed regime and its standard error are obtained as follows:

q1 <- modelObj::buildModelObj(model = ~ CD4_0*A1,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')
q2 <- modelObj::buildModelObj(model = ~ CD4_0 + CD4_6*A2,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')
q3 <- modelObj::buildModelObj(model = ~ CD4_0 + CD4_6 + CD4_12*A3,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')
d1 <- as.integer(x = {dataMDPF$CD4_0 < 250.0})
d2 <- d1 + {1L-d1}*as.integer(x = {dataMDPF$CD4_6 < 360.0})
d3 <- d2 + {1L-d2}*as.integer(x = {dataMDPF$CD4_12 < 300.0})
regime <- cbind(d1, d2, d3)
value_Q_md(moQ = list(q1, q2, q3), 
           data = dataMDPF,  
           y = dataMDPF$Y,  
           txName = c("A1","A2","A3"),
           regime = regime)
$fitQ
$fitQ[[1]]

Call:
lm(formula = YinternalY ~ CD4_0 + A1 + CD4_0:A1, data = data)

Coefficients:
(Intercept)        CD4_0           A1     CD4_0:A1  
    372.657        1.646      484.531       -1.945  


$fitQ[[2]]

Call:
lm(formula = YinternalY ~ CD4_0 + CD4_6 + A2 + CD4_6:A2, data = data)

Coefficients:
(Intercept)        CD4_0        CD4_6           A2     CD4_6:A2  
   344.8571       1.8542      -0.1218     500.8351      -1.6030  


$fitQ[[3]]

Call:
lm(formula = YinternalY ~ CD4_0 + CD4_6 + CD4_12 + A3 + CD4_12:A3, 
    data = data)

Coefficients:
(Intercept)        CD4_0        CD4_6       CD4_12           A3    CD4_12:A3  
   317.3743       2.0326       0.1371      -0.4478     603.5614      -1.9814  



$valueHat
[1] 1111.321

We see that the regression based estimator of the value under this regime is 1111.32 cells/mm\(^3\).



Inverse Propensity Weighted Estimator





The inverse probability weighted (IPW) estimator of the value of a fixed regime, \(\widehat{\mathcal{V}}_{IPW}(d)\), is defined as

\[ \widehat{\mathcal{V}}_{IPW}(d) = n^{-1} \sum_{i=1}^{n} \frac{\mathcal{C}_{d,i}Y_{i}}{\left\{\prod_{k=2}^{K}\pi_{d,k}(\overline{X}_{ki});\widehat{\gamma}_{k})\right\}\pi_{d,1}(X_{1i};\widehat{\gamma}_{1})}. \] where

\[ \mathcal{C}_{d,i} = \text{I}\{\overline{A}_{i} = \overline{d}(\overline{X}_{i})\} \]

is the indicator of whether or not all treatments actually received coicides with the options dictated by \(d\). For binary treatments \(\mathcal{A}_{k} = \{0,1\}\), \(k=1, \dots, K\)

\[ \pi_{d,k}(h_{k}; \gamma_{k}) = \sum_{l=1}^{\ell_{k}}\text{I}\{s_{k}(h_{k}) = l\} \prod_{a_{k} = 1}^{m_{kl}} \omega_{k,l}(h_{k},a_{k};\gamma_{k})\text{I}\{d_{k}(h_{k}) = a_{k}\} \]

Recall that \(\omega_{k}(h_{k},a_{k}; \gamma_{k})\) is a model for

\[ \omega_{k}(h_{k},a_{k}) = P(A_{k} = a_{k} | H_{k} = h_{k}), \quad k = 1, \ldots, K. \]



We define function value_IPW_md() to calculate the inverse probability weighted estimator of the value of a fixed regime.



#------------------------------------------------------------------------------#
# IPW estimator for value of a fixed tx regime with binary tx at each of the
#   K>1 decision points
#
# ASSUMPTIONS
#   each tx is binary coded as {0,1}
#
# INPUTS
#  moPS   : list of modeling objects for propensity score regressions
#           kth element corresponds to kth decision point
#  data   : data.frame containing all covariates and txs
#           *** each tx must be coded 0/1 ***
#  y      : outcome of interest
#  txName : vector of tx names (each tx must be coded as 0/1)
#  regime : matrix of 0/1 containing recommended txs
#           kth column corresponds to kth decision point
#
# RETURNS
#  a list containing
#     fitPS : list of value objects returned by propensity score regression
#  valueHat : estimated value
#------------------------------------------------------------------------------#
value_IPW_md <- function(moPS, data, y, txName, regime) {

  # number of decision points
  K <- length(x = moPS)

  # number of patients in data set
  n <- nrow(x = data)

  #### Propensity Score Regressions

  piStep <- function(moPS, data, y, regime) {
    # fit propensity score model
    fitObj <- modelObj::fit(object = moPS, data = data, response = y)

    # predict propensity score
    pk <- modelObj::predict(object = fitObj, newdata = data)
    if (is.matrix(x = pk)) pk <- pk[,ncol(x = pk),drop = TRUE]

    # pull propensity score corresponding to recommended treatment
    pidk <- pk*{regime == 1L} + {1.0 - pk}*{regime == 0L}

    return (list("fit" = fitObj, "pidk" = pidk))
  }

  fitPS <- list()
  pid <- matrix(data = 1.0, nrow = n, ncol = K)

  for (k in K:2L) {

    # use only individuals with more than 1 tx option
    one <- data[,txName[k-1L]] == 1L

    ## kth propensity score regression

    sk <- piStep(moPS = moPS[[ k ]], 
                 data = data[!one,,drop = FALSE],  
                 y = data[!one,txName[k]],
                 regime = regime[!one,k])

    fitPS[[ k ]] <- sk$fit
    pid[!one,k] <- sk$pidk

  }

  ## k=1 propensity score regression

  s1 <- piStep(moPS = moPS[[ 1L ]], 
               data = data,  
               y = data[,txName[1L]], 
               regime = regime[,1L])

  fitPS[[ 1L ]] <- s1$fit
  pid[,1L] <- s1$pidk

  # indicator of all txs received being those recommended
  cee <- apply(X = {data[,txName] == regime}, MARGIN = 1L, FUN = prod)

  # P(C == 1|H)
  pid <- apply(X = pid, MARGIN = 1L, FUN = prod)

  #### Value Estimate

  value <- mean(x = cee * y / pid)

  return( list("fitPS" = fitPS, "valueHat" = value) )
}


With value_IPW_md() defined (see Implementation tab), it is straightforward to estimate the value. See \(\pi_{k}(h_{k};\gamma_{k})\) in the sidebar for a review of the models and their basic diagnostics.


We compare the estimated value under the fixed regimes considered for the regression based estimator. Each tab below illustrates how to define the fixed regime under consideration and complete an analysis using the posited propensity score models.





The first fixed regime we consider is the static regime \(d = \{d_{1}(x_{1}) \equiv 1, d_{2}(\overline{x}_{2}) \equiv 1, d_{3}(\overline{x}_{3}) \equiv 1\}\). That is, all individuals are recommended to receive treatment \(A_{k} = 1\), \(k = 1, 2, 3\).

The estimated value of this fixed regime and its standard error are obtained as follows:

p1 <- modelObj::buildModelObj(model = ~ CD4_0,
                              solver.method = 'glm',
                              solver.args = list(family='binomial'),
                              predict.method = 'predict.glm',
                              predict.args = list(type='response'))
p2 <- modelObj::buildModelObj(model = ~ CD4_6,
                              solver.method = 'glm',
                              solver.args = list(family='binomial'),
                              predict.method = 'predict.glm',
                              predict.args = list(type='response'))
p3 <- modelObj::buildModelObj(model = ~ CD4_12,
                              solver.method = 'glm',
                              solver.args = list(family='binomial'),
                              predict.method = 'predict.glm',
                              predict.args = list(type='response'))
regime <- matrix(data = 1L, nrow = nrow(x = dataMDPF), ncol = 3L)
value_IPW_md(moPS = list(p1, p2, p3), 
             data = dataMDPF,  
             y = dataMDPF$Y,  
             txName = c("A1","A2","A3"),
             regime = regime)
$fitPS
$fitPS[[1]]

Call:  glm(formula = YinternalY ~ CD4_0, family = "binomial", data = data)

Coefficients:
(Intercept)        CD4_0  
   2.385956    -0.006661  

Degrees of Freedom: 999 Total (i.e. Null);  998 Residual
Null Deviance:      1316 
Residual Deviance: 1227     AIC: 1231

$fitPS[[2]]

Call:  glm(formula = YinternalY ~ CD4_6, family = "binomial", data = data)

Coefficients:
(Intercept)        CD4_6  
   1.234491    -0.004781  

Degrees of Freedom: 631 Total (i.e. Null);  630 Residual
Null Deviance:      608.5 
Residual Deviance: 577.8    AIC: 581.8

$fitPS[[3]]

Call:  glm(formula = YinternalY ~ CD4_12, family = "binomial", data = data)

Coefficients:
(Intercept)       CD4_12  
   0.929409    -0.003816  

Degrees of Freedom: 513 Total (i.e. Null);  512 Residual
Null Deviance:      622.5 
Residual Deviance: 609.1    AIC: 613.1


$valueHat
[1] 716.5866

We see that the inverse probability weighted estimator of the value under this regime is 716.59 cells/mm\(^3\). Note that 368 individuals in our data set followed this regime and thus contribute to the estimator.



The second fixed regime we consider is the static regime \(d = \{d_{1}(x_{1}) \equiv 0, d_{2}(\overline{x}_{2}) \equiv 0, d_{3}(\overline{x}_{3}) \equiv 0\}\). That is, all individuals are recommended to receive treatment \(A_{k} = 0\), \(k = 1, 2, 3\).

The estimated value of this fixed regime and its standard error are obtained as follows:

p1 <- modelObj::buildModelObj(model = ~ CD4_0,
                              solver.method = 'glm',
                              solver.args = list(family='binomial'),
                              predict.method = 'predict.glm',
                              predict.args = list(type='response'))
p2 <- modelObj::buildModelObj(model = ~ CD4_6,
                              solver.method = 'glm',
                              solver.args = list(family='binomial'),
                              predict.method = 'predict.glm',
                              predict.args = list(type='response'))
p3 <- modelObj::buildModelObj(model = ~ CD4_12,
                              solver.method = 'glm',
                              solver.args = list(family='binomial'),
                              predict.method = 'predict.glm',
                              predict.args = list(type='response'))
regime <- matrix(data = 0L, nrow = nrow(x = dataMDPF), ncol = 3L)
value_IPW_md(moPS = list(p1, p2, p3), 
             data = dataMDPF,  
             y = dataMDPF$Y,  
             txName = c("A1","A2","A3"),
             regime = regime)
$fitPS
$fitPS[[1]]

Call:  glm(formula = YinternalY ~ CD4_0, family = "binomial", data = data)

Coefficients:
(Intercept)        CD4_0  
   2.385956    -0.006661  

Degrees of Freedom: 999 Total (i.e. Null);  998 Residual
Null Deviance:      1316 
Residual Deviance: 1227     AIC: 1231

$fitPS[[2]]

Call:  glm(formula = YinternalY ~ CD4_6, family = "binomial", data = data)

Coefficients:
(Intercept)        CD4_6  
   1.234491    -0.004781  

Degrees of Freedom: 631 Total (i.e. Null);  630 Residual
Null Deviance:      608.5 
Residual Deviance: 577.8    AIC: 581.8

$fitPS[[3]]

Call:  glm(formula = YinternalY ~ CD4_12, family = "binomial", data = data)

Coefficients:
(Intercept)       CD4_12  
   0.929409    -0.003816  

Degrees of Freedom: 513 Total (i.e. Null);  512 Residual
Null Deviance:      622.5 
Residual Deviance: 609.1    AIC: 613.1


$valueHat
[1] 1094.317

We see that the inverse probability weighted estimator of the value under this regime is 1094.32 cells/mm\(^3\). Note that 118 individuals in our data set followed this regime and thus contribute to the estimator.



The final fixed regime we consider is a covariate dependent regime

\[ d = \left\{ \begin{array}{l} d_{1}(x_{1}) = \text{I}(\text{CD4_0} < 250 ~\text{cells}/\text{mm}^{3}) \\ d_{2}(\overline{x}_{2}) = d_{1}(x_{1}) + \{1-d_{1}(x_{1})\} ~ \text{I}(\text{CD4_6} < 360 ~\text{cells}/\text{mm}^{3}) \\ d_{3}(\overline{x}_{2}) = d_{2}(\overline{x}_{2}) + \{1-d_{2}(\overline{x}_{2})\} ~ \text{I}(\text{CD4_12} < 300 ~\text{cells}/\text{mm}^{3}) \end{array} ~~~. \right. \]

The estimated value of this fixed regime and its standard error are obtained as follows:

p1 <- modelObj::buildModelObj(model = ~ CD4_0,
                              solver.method = 'glm',
                              solver.args = list(family='binomial'),
                              predict.method = 'predict.glm',
                              predict.args = list(type='response'))
p2 <- modelObj::buildModelObj(model = ~ CD4_6,
                              solver.method = 'glm',
                              solver.args = list(family='binomial'),
                              predict.method = 'predict.glm',
                              predict.args = list(type='response'))
p3 <- modelObj::buildModelObj(model = ~ CD4_12,
                              solver.method = 'glm',
                              solver.args = list(family='binomial'),
                              predict.method = 'predict.glm',
                              predict.args = list(type='response'))
d1 <- as.integer(x = {dataMDPF$CD4_0 < 250.0})
d2 <- d1 + {1L-d1}*as.integer(x = {dataMDPF$CD4_6 < 360.0})
d3 <- d2 + {1L-d2}*as.integer(x = {dataMDPF$CD4_12 < 300.0})
regime <- cbind(d1, d2, d3)
value_IPW_md(moPS = list(p1, p2, p3), 
             data = dataMDPF,  
             y = dataMDPF$Y,  
             txName = c("A1","A2","A3"),
             regime = regime)
$fitPS
$fitPS[[1]]

Call:  glm(formula = YinternalY ~ CD4_0, family = "binomial", data = data)

Coefficients:
(Intercept)        CD4_0  
   2.385956    -0.006661  

Degrees of Freedom: 999 Total (i.e. Null);  998 Residual
Null Deviance:      1316 
Residual Deviance: 1227     AIC: 1231

$fitPS[[2]]

Call:  glm(formula = YinternalY ~ CD4_6, family = "binomial", data = data)

Coefficients:
(Intercept)        CD4_6  
   1.234491    -0.004781  

Degrees of Freedom: 631 Total (i.e. Null);  630 Residual
Null Deviance:      608.5 
Residual Deviance: 577.8    AIC: 581.8

$fitPS[[3]]

Call:  glm(formula = YinternalY ~ CD4_12, family = "binomial", data = data)

Coefficients:
(Intercept)       CD4_12  
   0.929409    -0.003816  

Degrees of Freedom: 513 Total (i.e. Null);  512 Residual
Null Deviance:      622.5 
Residual Deviance: 609.1    AIC: 613.1


$valueHat
[1] 1102.451

We see that the inverse probability weighted estimator of the value under this regime is 1102.45 cells/mm\(^3\). Note that 375 individuals in our data set followed this regime and thus contribute to the estimator.



Augmented Inverse Propensity Weighted Estimator





The augmented inverse probability weighted (AIPW) estimator of the value of a fixed regime, \(\widehat{\mathcal{V}}_{AIPW}(d)\), is defined as

\[ \begin{align} \widehat{\mathcal{V}}_{AIPW}(d) = n^{-1} \sum_{i=1}^{n} \Bigg[ &\left. \frac{\mathcal{C}_{d,i}Y_{i}}{\left\{\prod_{k=2}^{K}\pi_{d,k}(\overline{X}_{ki});\widehat{\gamma}_{k})\right\}\pi_{d,1}(X_{1i};\widehat{\gamma}_{1})} \right. \\ &\left. + \sum_{k=1}^{K} \left\{ \frac{\mathcal{C}_{\bar{d}_{k-1},i}}{\overline{\pi}_{d,k-1}(\overline{X}_{k-1,i};\widehat{\overline{\gamma}}_{k-1})} - \frac{\mathcal{C}_{\bar{d}_{k},i}}{\overline{\pi}_{d,k}(\overline{X}_{ki};\widehat{\overline{\gamma}}_{k})} \right\} Q_{k}^{d}\left\{\overline{X}_{ki},\overline{d}_{k-1}(\overline{X}_{k-1,i});\widehat{\beta}_{k}\right\}\right]. \end{align} \] where

\[ \mathcal{C}_{d} = \text{I}\{\overline{A} = \overline{d}(\overline{X})\}, \] is the indicator of whether or not all treatments actually received coicides with the options dictated by \(d\). And

\[ \mathcal{C}_{\bar{d}_{k}} = \text{I}\{\overline{A}_{k} = \overline{d}_{k}(\overline{X}_{k})\} \] is the indicator that the treatment options received through Decision \(k\) are consistent with those dictated by the first \(k\) rules in \(d\) and

\[ \pi_{d,k}(h_{k}; \gamma_{k}) = P(\mathcal{C}_{d_{k}}=1| H_{k} = h_{k}) = \sum_{l=1}^{\ell_{k}}\text{I}\{s_{k}(h_{k}) = l\} \prod_{a_{k} = 1}^{m_{kl}} \omega_{k,l}(h_{k},a_{k};\gamma_{k})\text{I}\{d_{k}(h_{k}) = a_{k}\} \]

Let \[ \overline{\pi}_{d,k}(\overline{X_{k}}) = \left\{\prod_{j=2}^{k} \pi_{d,j}(\overline{X}_j)\right\}\pi_{d,1}(X_{1}), \quad k=2, \dots, K \] with \(\overline{\pi}_{d,0} \equiv 1\).

We must posit a model \(\omega_{k}(h_{k},a_{k}; \gamma_{k})\) for each \(\omega(h_{k},a_{k}) = P(A_{k} = a_{k}|H_{k} = h_{k}), k = 1, \dots, K\). Further, we must posit a model \(Q_{k}^{d}(\overline{X}_{k},d_{k-1}(\overline{X}_{k-1,i});\beta_{k})\), \(k = 1, \dots, K\) for each \(E\{V_{k+1}^{d}(H_{K+1})|H_{k} = h_{k}, A_{k} = a_{k})\), \(k = 1, \dots, K\).



We define function value_AIPW_md() to calculate the augmented inverse probability weighted estimator of the value of a fixed regime.



#------------------------------------------------------------------------------#
# AIPW estimator for value of a fixed tx regime with binary tx at each of the
#   K>1 decision points
#
# ASSUMPTIONS
#   each tx is 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
#  moQ    : list of modeling objects for Q-function regressions
#           kth element corresponds to kth decision point
#  data   : data.frame containing all covariates and txs
#  y      : outcome of interest
#  txName : vector of tx names (each tx must be coded as 0/1)
#  regime : matrix of 0/1 containing recommended txs
#           kth column corresponds to kth decision point
#
# RETURNS
#  a list containing
#     fitPS : list of value objects returned by propensity score regression
#      fitQ : list of value objects returned by Q-function regressions
#  valueHat : estimated value
#------------------------------------------------------------------------------#
value_AIPW_md <- function(moPS, moQ, data, y, txName, regime) {

  # number of decision points
  K <- length(x = moPS)

  # number of patients in data set
  n <- nrow(x = data)

  #### Propensity Score and Q-function Regressions

  qStep <- function(moQ, data, y, txName, regime) {
    # fit Q-function model
    fitObj <- modelObj::fit(object = moQ, data = data, response = y)

    # set tx received to recommended tx
    data[,txName] <- regime

    # predict Q-function for recommended tx
    qk <- modelObj::predict(object = fitObj, newdata = data)

    return (list("fit" = fitObj, "qk" = qk))
  }

  piStep <- function(moPS, data, y, regime) {
    # fit propensity score model
    fitObj <- modelObj::fit(object = moPS, data = data, response = y)

    # predict propensity score
    pk <- modelObj::predict(object = fitObj, newdata = data)
    if (is.matrix(x = pk)) pk <- pk[,ncol(x = pk),drop = TRUE]

    # pull propensity score corresponding to recommended treatment
    pidk <- pk*{regime == 1L} + {1.0 - pk}*{regime == 0L}

    return (list("fit" = fitObj, "pidk" = pidk))
  }

  fitQ <- list()
  res <- matrix(data = 0.0, nrow = n, ncol = K+1L)

  fitPS <- list()
  pid <- matrix(data = 1.0, nrow = n, ncol = K)


  res[,K+1L] <- y

  for (k in K:2L) {

    # use only individuals with more than 1 tx option
    one <- data[,txName[k-1L]] == 1L

    ## kth Q-function regression

    sk <- qStep(moQ = moQ[[ k ]],
                data = data[!one,,drop = FALSE],
                y = res[!one,k+1L],
                txName = txName[k],
                regime = regime[!one,k])

    fitQ[[ k ]] <- sk$fit

    res[,k] <- res[,k+1L]
    res[!one,k] <- sk$qk

    ## kth propensity score regression

    sk <- piStep(moPS = moPS[[ k ]], 
                 data = data[!one,,drop = FALSE],  
                 y = data[!one,txName[k]],
                 regime = regime[!one,k])

    fitPS[[ k ]] <- sk$fit
    pid[!one,k] <- sk$pidk

  }

  # k=1 Q-function regression
  s1 <- qStep(moQ = moQ[[ 1L ]],
              data = data,
              y = res[,2L],
              txName = txName[1L],
              regime = regime[,1L])

  fitQ[[ 1L ]] <- s1$fit
  res[,1L] <- s1$qk

  # k=1 propensity score regression
  s1 <- piStep(moPS = moPS[[ 1L ]], 
               data = data,  
               y = data[,txName[1L]], 
               regime = regime[,1L])

  fitPS[[ 1L ]] <- s1$fit
  pid[,1L] <- s1$pidk

  # indicators of 1:k txs received being those recommended
  cee <- t(x = apply(X = {data[,txName] == regime}, 
                     MARGIN = 1L,  
                     FUN = cumprod))

  # P(Ck == 1|Hk)
  pid <- t(x = apply(X = pid, MARGIN = 1L, FUN = cumprod))

  #### Value Estimate

  value <- cee[,K] * y / pid[,K]
  value <- value + {1.0 - cee[,1L] / pid[,1L]} * res[,1L]

  for (k in 2L:K) {
    value <- value + {cee[,k-1L] / pid[,k-1L] - cee[,k] / pid[,k]}*res[,k]
  }

  value <- mean(x = value)

  return( list("fitPS" = fitPS, "fitQ" = fitQ, "valueHat" = value) )
}


With value_AIPW_md() defined (see Implementation tab), it is straightforward to estimate the value. See \(Q_{k}(h_{k},a_{k};\beta_{k})\) and \(\pi_{k}(h_{k};\gamma_{k})\) in the sidebar for a review of the models and their basic diagnostics.


We compare the estimated value under the fixed regimes considered for the regression based estimator. Each tab below illustrates how to define the fixed regime under consideration and complete an analysis using the posited propensity score models.





The first fixed regime we consider is the static regime \(d = \{d_{1}(x_{1}) \equiv 1, d_{2}(\overline{x}_{2}) \equiv 1, d_{3}(\overline{x}_{3}) \equiv 1\}\). That is, all individuals are recommended to receive treatment \(A_{k} = 1\), \(k = 1, 2, 3\).

The estimated value of this fixed regime and its standard error are obtained as follows:

p1 <- modelObj::buildModelObj(model = ~ CD4_0,
                              solver.method = 'glm',
                              solver.args = list(family='binomial'),
                              predict.method = 'predict.glm',
                              predict.args = list(type='response'))
p2 <- modelObj::buildModelObj(model = ~ CD4_6,
                              solver.method = 'glm',
                              solver.args = list(family='binomial'),
                              predict.method = 'predict.glm',
                              predict.args = list(type='response'))
p3 <- modelObj::buildModelObj(model = ~ CD4_12,
                              solver.method = 'glm',
                              solver.args = list(family='binomial'),
                              predict.method = 'predict.glm',
                              predict.args = list(type='response'))
q1 <- modelObj::buildModelObj(model = ~ CD4_0*A1,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')
q2 <- modelObj::buildModelObj(model = ~ CD4_0 + CD4_6*A2,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')
q3 <- modelObj::buildModelObj(model = ~ CD4_0 + CD4_6 + CD4_12*A3,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')
regime <- matrix(data = 1L, nrow = nrow(x = dataMDPF), ncol = 3L)
value_AIPW_md(moPS = list(p1, p2, p3), 
              moQ = list(q1, q2, q3),
              data = dataMDPF,  
              y = dataMDPF$Y,  
              txName = c("A1","A2","A3"),
              regime = regime)
$fitPS
$fitPS[[1]]

Call:  glm(formula = YinternalY ~ CD4_0, family = "binomial", data = data)

Coefficients:
(Intercept)        CD4_0  
   2.385956    -0.006661  

Degrees of Freedom: 999 Total (i.e. Null);  998 Residual
Null Deviance:      1316 
Residual Deviance: 1227     AIC: 1231

$fitPS[[2]]

Call:  glm(formula = YinternalY ~ CD4_6, family = "binomial", data = data)

Coefficients:
(Intercept)        CD4_6  
   1.234491    -0.004781  

Degrees of Freedom: 631 Total (i.e. Null);  630 Residual
Null Deviance:      608.5 
Residual Deviance: 577.8    AIC: 581.8

$fitPS[[3]]

Call:  glm(formula = YinternalY ~ CD4_12, family = "binomial", data = data)

Coefficients:
(Intercept)       CD4_12  
   0.929409    -0.003816  

Degrees of Freedom: 513 Total (i.e. Null);  512 Residual
Null Deviance:      622.5 
Residual Deviance: 609.1    AIC: 613.1


$fitQ
$fitQ[[1]]

Call:
lm(formula = YinternalY ~ CD4_0 + A1 + CD4_0:A1, data = data)

Coefficients:
(Intercept)        CD4_0           A1     CD4_0:A1  
 846.727405    -0.305528    10.459837     0.006195  


$fitQ[[2]]

Call:
lm(formula = YinternalY ~ CD4_0 + CD4_6 + A2 + CD4_6:A2, data = data)

Coefficients:
(Intercept)        CD4_0        CD4_6           A2     CD4_6:A2  
   923.4224       2.0372      -1.8156     -78.5422      -0.0542  


$fitQ[[3]]

Call:
lm(formula = YinternalY ~ CD4_0 + CD4_6 + CD4_12 + A3 + CD4_12:A3, 
    data = data)

Coefficients:
(Intercept)        CD4_0        CD4_6       CD4_12           A3    CD4_12:A3  
   317.3743       2.0326       0.1371      -0.4478     603.5614      -1.9814  



$valueHat
[1] 722.3174

We see that the augmented inverse probability weighted estimator of the value under this regime is 722.32 cells/mm\(^3\).



The second fixed regime we consider is the static regime \(d = \{d_{1}(x_{1}) \equiv 0, d_{2}(\overline{x}_{2}) \equiv 0, d_{3}(\overline{x}_{3}) \equiv 0\}\). That is, all individuals are recommended to receive treatment \(A_{k} = 0\), \(k = 1, 2, 3\).

The estimated value of this fixed regime and its standard error are obtained as follows:

p1 <- modelObj::buildModelObj(model = ~ CD4_0,
                              solver.method = 'glm',
                              solver.args = list(family='binomial'),
                              predict.method = 'predict.glm',
                              predict.args = list(type='response'))
p2 <- modelObj::buildModelObj(model = ~ CD4_6,
                              solver.method = 'glm',
                              solver.args = list(family='binomial'),
                              predict.method = 'predict.glm',
                              predict.args = list(type='response'))
p3 <- modelObj::buildModelObj(model = ~ CD4_12,
                              solver.method = 'glm',
                              solver.args = list(family='binomial'),
                              predict.method = 'predict.glm',
                              predict.args = list(type='response'))
q1 <- modelObj::buildModelObj(model = ~ CD4_0*A1,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')
q2 <- modelObj::buildModelObj(model = ~ CD4_0 + CD4_6*A2,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')
q3 <- modelObj::buildModelObj(model = ~ CD4_0 + CD4_6 + CD4_12*A3,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')
regime <- matrix(data = 0L, nrow = nrow(x = dataMDPF), ncol = 3L)
value_AIPW_md(moPS = list(p1, p2, p3), 
              moQ = list(q1, q2, q3),
              data = dataMDPF,  
              y = dataMDPF$Y,  
              txName = c("A1","A2","A3"),
              regime = regime)
$fitPS
$fitPS[[1]]

Call:  glm(formula = YinternalY ~ CD4_0, family = "binomial", data = data)

Coefficients:
(Intercept)        CD4_0  
   2.385956    -0.006661  

Degrees of Freedom: 999 Total (i.e. Null);  998 Residual
Null Deviance:      1316 
Residual Deviance: 1227     AIC: 1231

$fitPS[[2]]

Call:  glm(formula = YinternalY ~ CD4_6, family = "binomial", data = data)

Coefficients:
(Intercept)        CD4_6  
   1.234491    -0.004781  

Degrees of Freedom: 631 Total (i.e. Null);  630 Residual
Null Deviance:      608.5 
Residual Deviance: 577.8    AIC: 581.8

$fitPS[[3]]

Call:  glm(formula = YinternalY ~ CD4_12, family = "binomial", data = data)

Coefficients:
(Intercept)       CD4_12  
   0.929409    -0.003816  

Degrees of Freedom: 513 Total (i.e. Null);  512 Residual
Null Deviance:      622.5 
Residual Deviance: 609.1    AIC: 613.1


$fitQ
$fitQ[[1]]

Call:
lm(formula = YinternalY ~ CD4_0 + A1 + CD4_0:A1, data = data)

Coefficients:
(Intercept)        CD4_0           A1     CD4_0:A1  
    318.206        1.754      538.981       -2.053  


$fitQ[[2]]

Call:
lm(formula = YinternalY ~ CD4_0 + CD4_6 + A2 + CD4_6:A2, data = data)

Coefficients:
(Intercept)        CD4_0        CD4_6           A2     CD4_6:A2  
   318.0672       1.9303      -0.1408     527.2873      -1.6443  


$fitQ[[3]]

Call:
lm(formula = YinternalY ~ CD4_0 + CD4_6 + CD4_12 + A3 + CD4_12:A3, 
    data = data)

Coefficients:
(Intercept)        CD4_0        CD4_6       CD4_12           A3    CD4_12:A3  
   317.3743       2.0326       0.1371      -0.4478     603.5614      -1.9814  



$valueHat
[1] 1085.224

We see that the augmented inverse probability weighted estimator of the value under this regime is 1085.22 cells/mm\(^3\).



The final fixed regime we consider is a covariate dependent regime

\[ d = \left\{ \begin{array}{l} d_{1}(x_{1}) = \text{I}(\text{CD4_0} < 250 ~\text{cells}/\text{mm}^{3}) \\ d_{2}(\overline{x}_{2}) = d_{1}(x_{1}) + \{1-d_{1}(x_{1})\} ~ \text{I}(\text{CD4_6} < 360 ~\text{cells}/\text{mm}^{3}) \\ d_{3}(\overline{x}_{2}) = d_{2}(\overline{x}_{2}) + \{1-d_{2}(\overline{x}_{2})\} ~ \text{I}(\text{CD4_12} < 300 ~\text{cells}/\text{mm}^{3}) \end{array} ~~~. \right. \]

The estimated value of this fixed regime and its standard error are obtained as follows:

p1 <- modelObj::buildModelObj(model = ~ CD4_0,
                              solver.method = 'glm',
                              solver.args = list(family='binomial'),
                              predict.method = 'predict.glm',
                              predict.args = list(type='response'))
p2 <- modelObj::buildModelObj(model = ~ CD4_6,
                              solver.method = 'glm',
                              solver.args = list(family='binomial'),
                              predict.method = 'predict.glm',
                              predict.args = list(type='response'))
p3 <- modelObj::buildModelObj(model = ~ CD4_12,
                              solver.method = 'glm',
                              solver.args = list(family='binomial'),
                              predict.method = 'predict.glm',
                              predict.args = list(type='response'))
q1 <- modelObj::buildModelObj(model = ~ CD4_0*A1,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')
q2 <- modelObj::buildModelObj(model = ~ CD4_0 + CD4_6*A2,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')
q3 <- modelObj::buildModelObj(model = ~ CD4_0 + CD4_6 + CD4_12*A3,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')
d1 <- as.integer(x = {dataMDPF$CD4_0 < 250.0})
d2 <- d1 + {1L-d1}*as.integer(x = {dataMDPF$CD4_6 < 360.0})
d3 <- d2 + {1L-d2}*as.integer(x = {dataMDPF$CD4_12 < 300.0})
regime <- cbind(d1, d2, d3)
value_AIPW_md(moPS = list(p1, p2, p3), 
              moQ = list(q1, q2, q3),
              data = dataMDPF,  
              y = dataMDPF$Y,  
              txName = c("A1","A2","A3"),
              regime = regime)
$fitPS
$fitPS[[1]]

Call:  glm(formula = YinternalY ~ CD4_0, family = "binomial", data = data)

Coefficients:
(Intercept)        CD4_0  
   2.385956    -0.006661  

Degrees of Freedom: 999 Total (i.e. Null);  998 Residual
Null Deviance:      1316 
Residual Deviance: 1227     AIC: 1231

$fitPS[[2]]

Call:  glm(formula = YinternalY ~ CD4_6, family = "binomial", data = data)

Coefficients:
(Intercept)        CD4_6  
   1.234491    -0.004781  

Degrees of Freedom: 631 Total (i.e. Null);  630 Residual
Null Deviance:      608.5 
Residual Deviance: 577.8    AIC: 581.8

$fitPS[[3]]

Call:  glm(formula = YinternalY ~ CD4_12, family = "binomial", data = data)

Coefficients:
(Intercept)       CD4_12  
   0.929409    -0.003816  

Degrees of Freedom: 513 Total (i.e. Null);  512 Residual
Null Deviance:      622.5 
Residual Deviance: 609.1    AIC: 613.1


$fitQ
$fitQ[[1]]

Call:
lm(formula = YinternalY ~ CD4_0 + A1 + CD4_0:A1, data = data)

Coefficients:
(Intercept)        CD4_0           A1     CD4_0:A1  
    372.657        1.646      484.531       -1.945  


$fitQ[[2]]

Call:
lm(formula = YinternalY ~ CD4_0 + CD4_6 + A2 + CD4_6:A2, data = data)

Coefficients:
(Intercept)        CD4_0        CD4_6           A2     CD4_6:A2  
   344.8571       1.8542      -0.1218     500.8351      -1.6030  


$fitQ[[3]]

Call:
lm(formula = YinternalY ~ CD4_0 + CD4_6 + CD4_12 + A3 + CD4_12:A3, 
    data = data)

Coefficients:
(Intercept)        CD4_0        CD4_6       CD4_12           A3    CD4_12:A3  
   317.3743       2.0326       0.1371      -0.4478     603.5614      -1.9814  



$valueHat
[1] 1114.073

We see that the augmented inverse probability weighted estimator of the value under this regime is 1114.07 cells/mm^3.



Comparison of Estimators





\(d = \{d_{1}(x_{1}) \equiv 1, d_{2}(\overline{x}_{2}) \equiv 1, d_{3}(\overline{x}_{3}) \equiv 1\}\)

The table below compares the estimated value obtained for the toy data set using each method discussed in the chapter. All values are in units of \(\text{cells}/\text{mm}^3\).

\(\widehat{\mathcal{V}}_{Q}(d)\) \(\widehat{\mathcal{V}}_{IPW}(d)\) \(\widehat{\mathcal{V}}_{AIPW}(d)\)
723.28 716.59 722.32
\(d = \{d_{1}(x_{1}) \equiv 0, d_{2}(\overline{x}_{2}) \equiv 0, d_{3}(\overline{x}_{3}) \equiv 0\}\)

The table below compares the estimated value obtained for the toy data set using each method discussed in the chapter. All values are in units of \(\text{cells}/\text{mm}^3\).

\(\widehat{\mathcal{V}}_{Q}(d)\) \(\widehat{\mathcal{V}}_{IPW}(d)\) \(\widehat{\mathcal{V}}_{AIPW}(d)\)
1102.81 1094.32 1085.22

\[ d = \left\{ \begin{array}{l} d_{1}(x_{1}) = \text{I}(\text{CD4_0} < 250 ~\text{cells}/\text{mm}^{3}) \\ d_{2}(\overline{x}_{2}) = d_{1}(x_{1}) + \{1-d_{1}(x_{1})\} ~ \text{I}(\text{CD4_6} < 360 ~\text{cells}/\text{mm}^{3}) \\ d_{3}(\overline{x}_{2}) = d_{2}(\overline{x}_{2}) + \{1-d_{2}(\overline{x}_{2})\} ~ \text{I}(\text{CD4_12} < 300 ~\text{cells}/\text{mm}^{3}) \end{array} ~~~. \right. \]

The table below compares the estimated value obtained for the toy data set using each method discussed in the chapter. All values are in units of \(\text{cells}/\text{mm}^3\).

\(\widehat{\mathcal{V}}_{Q}(d)\) \(\widehat{\mathcal{V}}_{IPW}(d)\) \(\widehat{\mathcal{V}}_{AIPW}(d)\)
1111.32 1102.45 1114.07



We have carried out a simulation study to evaluate the performances of the presented methods. We generated 1000 Monte Carlo data sets, each with \(n=1000\).





\(d = \{d_{1}(x_{1}) \equiv 1, d_{2}(\overline{x}_{2}) \equiv 1, d_{3}(\overline{x}_{3}) \equiv 1\}\)

The table below compares the Monte Carlo average and standard deviation for the estimated value obtained using each method discussed in the chapter. All values are in units of \(\text{cells}/\text{mm}^3\).

\(\widehat{\mathcal{V}}_{Q}(d)\) \(\widehat{\mathcal{V}}_{IPW}(d)\) \(\widehat{\mathcal{V}}_{AIPW}(d)\)
719.88 (3.76) 718.37 (5.79) 718.47 (3.83)
\(d = \{d_{1}(x_{1}) \equiv 0, d_{2}(\overline{x}_{2}) \equiv 0, d_{3}(\overline{x}_{3}) \equiv 0\}\)

The table below compares the Monte Carlo average and standard deviation for the estimated value obtained using each method discussed in the chapter. All values are in units of \(\text{cells}/\text{mm}^3\).

\(\widehat{\mathcal{V}}_{Q}(d)\) \(\widehat{\mathcal{V}}_{IPW}(d)\) \(\widehat{\mathcal{V}}_{AIPW}(d)\)
1113.46 (7.10) 1107.06 (9.87) 1107.13 (11.36)

\[ d = \left\{ \begin{array}{l} d_{1}(x_{1}) = \text{I}(\text{CD4_0} < 250 ~\text{cells}/\text{mm}^{3}) \\ d_{2}(\overline{x}_{2}) = d_{1}(x_{1}) + \{1-d_{1}(x_{1})\} ~ \text{I}(\text{CD4_6} < 360 ~\text{cells}/\text{mm}^{3}) \\ d_{3}(\overline{x}_{2}) = d_{2}(\overline{x}_{2}) + \{1-d_{2}(\overline{x}_{2})\} ~ \text{I}(\text{CD4_12} < 300 ~\text{cells}/\text{mm}^{3}) \end{array} ~~~. \right. \]

The table below compares the Monte Carlo average and standard deviation for the estimated value obtained using each method discussed in the chapter. All values are in units of \(\text{cells}/\text{mm}^3\).

\(\widehat{\mathcal{V}}_{Q}(d)\) \(\widehat{\mathcal{V}}_{IPW}(d)\) \(\widehat{\mathcal{V}}_{AIPW}(d)\)
1119.68 (5.82) 1119.98 (20.48) 1119.64 (5.86)


\(Q_{k}(h_{k}, a_{k};\beta_{k})\)



For Chapters 6 and 7, we consider a single model for each of the three propensity scores. It is not our objective to demonstrate a definitive analysis that one might do in practice but to illustrate how to apply the methods. The posited models are intentionally kept simple and likely to be familiar to most readers. By this, we avoid adding any additional complexity to the discussion. In practice, analysist would perform model and variable selection techniques, etc., to arrive at a posited model.


Click on each tab below to see the model and basic model diagnostic steps. For all of the methods discussed here, the Q-function models are fit using a backward iterative approach, which we also take here.



The posited model for \(Q_{3}(h_{3},a_{3}) = E(Y|H_{3} = h_{3}, A_{3} = a_{3})\) is misspecified as

\[ Q_{3}(h_{3},a_{3};{\beta}_{3}) = {\beta}_{30} + {\beta}_{31} \text{CD4_0} + {\beta}_{32} \text{CD4_6} + {\beta}_{33} \text{CD4_12} + a_{3}~({\beta}_{34} + {\beta}_{35} \text{CD4_12}). \]


Modeling Object


The parameters of this model will be estimated using ordinary least squares via R’s stats::lm() function. Predictions will be made using stats::predict.lm(), which by default returns predictions on the scale of the response variable as is the requirement of package DynTxRegime.

The modeling object for this regression step is

q3 <- modelObj::buildModelObj(model = ~ CD4_0 + CD4_6 + CD4_12*A3,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')

Model Diagnostics


Recall that individuals for whom only one treatment option is available should not be included in the regression analysis; i.e., individuals that received treatment \(A_{2} = 1\) remained on treatment 1 at the third stage.

oneA3 <- dataMDPF$A2 == 1L

For \(Q_{3}(h_{3},a_{3};{\beta}_{3})\) the regression can be completed as follows

Q3 <- modelObj::fit(object = q3, 
                    data = dataMDPF[!oneA3,],  
                    response = dataMDPF$Y[!oneA3])
Q3 <- modelObj::fitObject(object = Q3)
Q3

Call:
lm(formula = YinternalY ~ CD4_0 + CD4_6 + CD4_12 + A3 + CD4_12:A3, 
    data = data)

Coefficients:
(Intercept)        CD4_0        CD4_6       CD4_12           A3    CD4_12:A3  
   317.3743       2.0326       0.1371      -0.4478     603.5614      -1.9814  

where for convenience we have made use of modelObj’s fitObject() function to strip away the modeling object framework making Q3 an “lm” object.

In examining the diagnostic plots defined for “lm” objects

par(mfrow = c(2,2))
graphics::plot(x = Q3)

There is some indication of non-normality seen in the Normal Q-Q plot and the presence of outliers in the Residuals vs Leverage plot are consistent with the fact that this model is misspecified.

The methods discussed in this chapter use the backward iterative algorithm to obtain parameter estimates for the Q-function models. Thus, the pseudo-outcome is used in the second stage regression that follows. Recall, the pseudo-outcome at this decision point is an individuals expected outcome had they received treatment according the optimal treatment regime at Decision 3. The optimal regime is that which maximizes the Q-function, and thus we calculate the pseudo-outcome as follows. Notice that for individuals that received \(A_{2} = 1\), the pseudo-outcome is taken to be the final response .

A3 <- dataMDPF$A3

dataMDPF$A3 <- 0L
Q30 <- stats::predict(object = Q3, newdata = dataMDPF[!oneA3,])
dataMDPF$A3 <- 1L
Q31 <- stats::predict(object = Q3, newdata = dataMDPF[!oneA3,])
dataMDPF$A3 <- A3

V3 <- dataMDPF$Y
V3[!oneA3] <- pmax(Q30, Q31)


The posited model for \(Q_{2}(h_{2},a_{2}) = E\{V_{3}(H_{3})|H_{2} = h_{2}, A_{2} = a_{2}\}\) is misspecified as

\[ Q_{2}(h_{2},a_{2};\beta_{2}) = \beta_{20} + \beta_{21} \text{CD4_0} + \beta_{22} \text{CD4_6} + a_{2}~(\beta_{23} + \beta_{24} \text{CD4_6}). \]


Modeling Object


The parameters of this model will be estimated using ordinary least squares via R’s stats::lm() function. Predictions will be made using stats::predict.lm(), which by default returns predictions on the scale of the response variable.

The modeling objects for this regression step is

q2 <- modelObj::buildModelObj(model = ~ CD4_0 + CD4_6*A2,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')

Model Diagnostics


Recall that individuals for whom only one treatment option is available should not be included in the regression analysis; i.e., individuals that received treatment \(A_{1} = 1\) remained on treatment 1 at the second stage.

oneA2 <- dataMDPF$A1 == 1L

The pseudo-outcome, \(V_{3}(H_{3})\), was calculated under the previous tab. For \(Q_{2}(h_{2},a_{2};{\beta}_{2})\) the regression can be completed as follows

Q2 <- modelObj::fit(object = q2, 
                    data = dataMDPF[!oneA2,],  
                    response = V3[!oneA2])
Q2 <- modelObj::fitObject(object = Q2)
Q2

Call:
lm(formula = YinternalY ~ CD4_0 + CD4_6 + A2 + CD4_6:A2, data = data)

Coefficients:
(Intercept)        CD4_0        CD4_6           A2     CD4_6:A2  
   344.9733       1.8566      -0.1238     500.7085      -1.6028  

where for convenience we have again made use of modelObj’s fitObject() function to strip away the modeling object framework making Q2 an “lm” object.

In examining the diagnostic plots defined for “lm” objects

par(mfrow = c(2,2))
graphics::plot(x = Q2)

The non-normality of the residuals is consistent with the fact that this model is misspecified.

The methods discussed in this chapter use the backward iterative algorithm to obtain parameter estimates for the Q-function models. Thus, the pseudo-outcome is used in the first stage regression that follows. Recall, the pseudo-outcome at this decision point is an individuals expected outcome had they received treatment according the optimal treatment regime at Decisions 2 and 3. The optimal regime is that which maximizes the Q-function, and thus we calculate the pseudo-outcome as follows. Note that for individuals that received \(A_{1} = 1\) and were thus not included in the model fitting procedure, the pseudo-outcome is equal to the Decision 3 pseudo-outcome.

A2 <- dataMDPF$A2

dataMDPF$A2 <- 0L
Q20 <- stats::predict(object = Q2, newdata = dataMDPF[!oneA2,])
dataMDPF$A2 <- 1L
Q21 <- stats::predict(object = Q2, newdata = dataMDPF[!oneA2,])
dataMDPF$A2 <- A2

V2 <- V3
V2[!oneA2] <- pmax(Q20, Q21)


The posited model for \(Q_{1}(h_{1},a_{1}) = E\{V_{2}(H_{2})|H_{1} = h_{1}, A_{1} = a_{1}\}\) is taken to be

\[ Q_{1}(h_{1},a_{1};\beta_{1}) = \beta_{10} + \beta_{11} \text{CD4_0} + a_{1}~(\beta_{12} + \beta_{13} \text{CD4_0}). \]


Modeling Object


The parameters of this model will be estimated using ordinary least squares via R’s stats::lm() function. Predictions will be made using stats::predict.lm(), which by default returns predictions on the scale of the response variable as is the requirment of the DynTxRegime package.

The modeling objects for this regression step is

q1 <- modelObj::buildModelObj(model = ~ CD4_0*A1,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')

Model Diagnostics


Recall that \(Q_{1}(h_{1},a_{1})\) is the expectation of the value, \(V_{2}(H_{2})\); we use the results discussed in the previous regression analysis. For \(Q_{1}(h_{1},a_{1};\beta_{1})\) the regression can be completed as follows

Q1 <- modelObj::fit(object = q1, data = dataMDPF, response = V2)
Q1 <- modelObj::fitObject(object = Q1)
Q1

Call:
lm(formula = YinternalY ~ CD4_0 + A1 + CD4_0:A1, data = data)

Coefficients:
(Intercept)        CD4_0           A1     CD4_0:A1  
    379.564        1.632      477.623       -1.932  

where for convenience we have again made use of modelObj’s fitObject() function to strip away the modeling object framework making Q1 an “lm” object.

In examining the diagnostic plots defined for “lm” objects

par(mfrow = c(2,2))
graphics::plot(x = Q1)

The non-normality seen in the Normal Q-Q plot and the presence of outliers in the Residuals vs Leverage plot are consistent with the fact that this model is misspecified.

\(\omega_{k}(h_{k},a_{k};\gamma_{k })\)



For Chapters 6 and 7, we consider a single model for each of the three propensity scores. It is not our objective to demonstrate a definitive analysis that one might do in practice but to illustrate how to apply the methods. The posited models are intentionally kept simple and likely to be familiar to most readers. By this, we avoid adding any additional complexity to the discussion. In practice, analysist would perform model and variable selection techniques, etc., to arrive at a posited model.



Click on each tab below to see the model and basic model diagnostic steps.



The posited model for the first decision point is the model used to generate the data

\[ \omega_{1}(h_{1},1;\gamma_{1}) = \frac{\exp(\gamma_{10} + \gamma_{11}~\text{CD4_0})}{1+\exp(\gamma_{10} + \gamma_{11}~\text{CD4_0})}. \]


Modeling Object


The parameters of this model will be estimated using maximum likelihood via R’s stats::glm() function. Predictions will be made using stats::predict.glm(), which by default returns predictions on the scale of the linear predictors. For the implementations presented in this chapter, prediction for the propensity score must always be on the scale of the response variable, i.e., the probabilities. The modeling object for this model is specified as

p1 <- modelObj::buildModelObj(model = ~ CD4_0,
                              solver.method = 'glm',
                              solver.args = list(family='binomial'),
                              predict.method = 'predict.glm',
                              predict.args = list(type='response'))

Model Diagnostics


For \(\omega_{1}(h_{1},1;\gamma_{1})\) the regression can be completed as follows:

PS1 <- modelObj::fit(object = p1, data = dataMDPF, response = dataMDPF$A1)
PS1 <- modelObj::fitObject(object = PS1)
PS1

Call:  glm(formula = YinternalY ~ CD4_0, family = "binomial", data = data)

Coefficients:
(Intercept)        CD4_0  
   2.385956    -0.006661  

Degrees of Freedom: 999 Total (i.e. Null);  998 Residual
Null Deviance:      1316 
Residual Deviance: 1227     AIC: 1231

where for convenience we have made use of modelObj’s fitObject() function to strip away the modeling object framework making PS1 an object of class “glm.”

Though we know this model to be correctly specified, let’s examine the regression results.

summary(object = PS1)

Call:
glm(formula = YinternalY ~ CD4_0, family = "binomial", data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9111  -0.9416  -0.7097   1.2143   2.1092  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.3859559  0.3359581   7.102 1.23e-12 ***
CD4_0       -0.0066608  0.0007601  -8.763  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1315.8  on 999  degrees of freedom
Residual deviance: 1227.3  on 998  degrees of freedom
AIC: 1231.3

Number of Fisher Scoring iterations: 4

In comparing the null deviance (1315.8) and the residual deviance (1227.3), we see that including the independent variable, \(\text{CD4_0}\) does reduce the deviance. This result is consistent with the fact that it is the model used to generate the data.



The posited model for the second decision point is the model used to generate the data for individuals in subset \(s_{2}(h_{2}) = 2\), for whom more than one treatment option was available

\[ \omega_{2,2}(h_{2},1;\gamma_{2}) = \frac{\exp(\gamma_{20} + \gamma_{21}~\text{CD4_6})}{1+\exp(\gamma_{20} + \gamma_{21}~\text{CD4_6})}. \]


Modeling Object


The parameters of this model will be estimated using maximum likelihood via R’s stats::glm() function. Predictions will be made using stats::predict.glm(), which by default returns predictions on the scale of the linear predictors. For the implementations presented in this chapter, prediction for the propensity score must always be on the scale of the response variable, i.e., the probabilities. The modeling object for this model is specified as

p2 <- modelObj::buildModelObj(model = ~ CD4_6,
                              solver.method = 'glm',
                              solver.args = list(family='binomial'),
                              predict.method = 'predict.glm',
                              predict.args = list(type='response'))

Model Diagnostics


Recall that individuals for whom only one treatment option is available should not be included in the regression analysis; i.e., individuals that received treatment \(A_{1} = 1\) received treatment \(A_{2} = 1\) with probability 1.

oneA2 <- dataMDPF$A1 == 1L

For \(\omega_{2,2}(h_{2},1;\gamma_{2})\) the regression can be completed as follows:

PS2 <- modelObj::fit(object = p2, 
                     data = dataMDPF[!oneA2,],  
                     response = dataMDPF$A2[!oneA2])
PS2 <- modelObj::fitObject(object = PS2)
PS2

Call:  glm(formula = YinternalY ~ CD4_6, family = "binomial", data = data)

Coefficients:
(Intercept)        CD4_6  
   1.234491    -0.004781  

Degrees of Freedom: 631 Total (i.e. Null);  630 Residual
Null Deviance:      608.5 
Residual Deviance: 577.8    AIC: 581.8

where for convenience we have made use of modelObj’s fitObject() function to strip away the modeling object framework making PS2 an object of class “glm.”

Though we know this model to be correctly specified, let’s examine the regression results.

summary(object = PS2)

Call:
glm(formula = YinternalY ~ CD4_6, family = "binomial", data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.2815  -0.6751  -0.5572  -0.4062   2.4387  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.2344914  0.5042630   2.448   0.0144 *  
CD4_6       -0.0047811  0.0009013  -5.304 1.13e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 608.51  on 631  degrees of freedom
Residual deviance: 577.83  on 630  degrees of freedom
AIC: 581.83

Number of Fisher Scoring iterations: 4

In comparing the null deviance (608.5) and the residual deviance (577.8), we see that including the independent variable, \(\text{CD4_6}\) reduces the deviance. This result is consistent with the fact that it is the model used to generate the data.



The posited model for the final decision point is the model used to generate the data for individuals in subset \(s_{3}(h_{3}) = 2\), for whom more than one treatment option was available

\[ \omega_{3,2}(h_{3},1;\gamma_{3}) = \frac{\exp(\gamma_{30} + \gamma_{31}~\text{CD4_12})}{1+\exp(\gamma_{30} + \gamma_{31}~\text{CD4_12})}. \]


Modeling Object


The parameters of this model will be estimated using maximum likelihood via R’s stats::glm() function. Predictions will be made using stats::predict.glm(), which by default returns predictions on the scale of the linear predictors. For the implementations presented in this chapter, prediction for the propensity score must always be on the scale of the response variable, i.e., the probabilities. The modeling object for this model is specified as

p3 <- modelObj::buildModelObj(model = ~ CD4_12,
                              solver.method = 'glm',
                              solver.args = list(family='binomial'),
                              predict.method = 'predict.glm',
                              predict.args = list(type='response'))

Model Diagnostics


Recall that individuals for whom only one treatment option is available should not be included in the regression analysis; i.e., individuals that received treatment \(A_{2} = 1\) were assigned treatment \(A_{3} = 1\) with probability 1.

oneA3 <- dataMDPF$A2 == 1L

For \(\omega_{3,2}(h_{3},1;\gamma_{3})\) the regression can be completed as follows:

PS3 <- modelObj::fit(object = p3, 
                     data = dataMDPF[!oneA3,],  
                     response = dataMDPF$A3[!oneA3])
PS3 <- modelObj::fitObject(object = PS3)
PS3

Call:  glm(formula = YinternalY ~ CD4_12, family = "binomial", data = data)

Coefficients:
(Intercept)       CD4_12  
   0.929409    -0.003816  

Degrees of Freedom: 513 Total (i.e. Null);  512 Residual
Null Deviance:      622.5 
Residual Deviance: 609.1    AIC: 613.1

where for convenience we have made use of modelObj’s fitObject() function to strip away the modeling object framework making PS3 an object of class “glm.”

Though we know this model to be correctly specified, let’s examine the regression results.

summary(object = PS3)

Call:
glm(formula = YinternalY ~ CD4_12, family = "binomial", data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.2985  -0.8629  -0.7505   1.3664   1.8875  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.929409   0.508608   1.827 0.067646 .  
CD4_12      -0.003816   0.001069  -3.569 0.000359 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 622.45  on 513  degrees of freedom
Residual deviance: 609.14  on 512  degrees of freedom
AIC: 613.14

Number of Fisher Scoring iterations: 4

In comparing the null deviance (622.5) and the residual deviance (609.1), we see that including the independent variable, \(\text{CD4_12}\) reduces the deviance. This result is consistent with the fact that it is the model used to generate the data.