anilbattalahalli / smartutils

Import the library

install.packages("devtools")
library(devtools)

install_github("AnilBattalahalli/SMARTutils")
library(SMARTutils)

Generate data from conditional parameters

data <- cSMART.dgen()

##  DTR      Treatment Mean      Treatment Variance      Treatment Correlation 
## ______    ______________      __________________      _____________________ 
##  1, 1         8.6           100.84          0.1769139 
##  1,-1         4.4           113.44          0.2374824 
## -1, 1         -0.8          105.76          0.1603631 
## -1,-1         3.2           100.16          0.08146965 
## ______    ______________      __________________      _____________________ 
## 
## True Beta:
## 
## Beta[0] - 3.85 
## Beta[1] - 2.65 
## Beta[2] - 0.05 
## Beta[3] - 2.05

head(data)

##   i j A1 R A2         Y
## 1 1 1 -1 0  1 -3.452152
## 2 1 2 -1 0  1 -2.735390
## 3 1 3 -1 0  1 -1.758195
## 4 1 4 -1 0  1  8.243450
## 5 1 5 -1 0  1  3.232144
## 6 2 1 -1 0  1 13.987650

Get the estimates for the regression parameters

report <- cSMART.mm(Y~a1+a2+I(a1*a2), data, verbose=T)

## Parameter     Estimate    Std.Err     Z Score     Pr(>|z|) 
## ___________   --------    -------     -------     ------- 
## (Intercept)   4.47717     0.36671     12.20888    2.786906e-34 
## a1            2.74937     0.36671     7.497308    6.514162e-14 
## a2            -0.22003    0.36671     -0.5999946      0.5485099 
## I(a1 * a2)    1.10251     0.36671     3.006454    0.00264314 
## ___________   --------    -------     -------     ------- 
## 
## Marginal Mean Model:  Y ~ a1 + a2 + I(a1 * a2) 
## 
## Working covariance structure: 'EXCH'  (Homogeneous-Exchangeable covariance structure)
## Variance      102.2465
## Correlation   0.1591008
## 
## Variance-Covariance matrix of the estimates
##              (Intercept)           a1           a2   I(a1 * a2)
## (Intercept)  0.134479319  0.003630063 -0.013093342 -0.004863585
## a1           0.003630063  0.134479319 -0.004863585 -0.013093342
## a2          -0.013093342 -0.004863585  0.134479319  0.003630063
## I(a1 * a2)  -0.004863585 -0.013093342  0.003630063  0.134479319

smartutils's Issues

Specification of outcome variable in cSMART.mm

I'm currently trying to specify a particular variable, y_adj, as my outcome variable in cSMART.mm() and I think there's something off with how the function handles the argument Y = "y_adj". I specify it as follows:

report <- cSMART.mm(Y ~ A1 + A2 + I(A1 * A2), df1, verbose = TRUE, covstr = 'IND', i = "school", Y = "y_adj")

but that provides the same output as when I use:

report <- cSMART.mm(Y ~ A1 + A2 + I(A1 * A2), df1, verbose = TRUE, covstr = 'IND', i = "school", Y = "Y")

(and note that it provides an error when I start with cSMART.mm(y_adj ~ ...). When I instead recode y_adj as simply Y, I obtain different (and presumably, correct) output.

Wrong probability for the Bernoulli random variable

    if (A1 == 1) {
      R <-  rbinom(1, 1, p_1)
      if (R == 1) {
        treat_code <- "1,1,."
        A2 <- 1
      } else if (R == 0) {
        A2 <-  (2 * rbinom(1, 1, p_A2)) - 1
        treat_code <- sprintf("1,0,%d", A2)
      }
    } else if (A1 == -1) {
      R <-  rbinom(1, 1 , p_1)
      if (R == 1) {
        treat_code <- "-1,1,."
        A2 <- 1
      }
      else if (R == 0) {
        A2 <-  (2 * rbinom(1, 1, p_A2)) - 1
        treat_code <- sprintf("-1,0,%d", A2)
      }
    }

Nullify A2 when R=1

      if (R == 1) {
        treat_code <- "-1,1,."
        A2 <- 1
      }

and

     if (R == 1) {
        treat_code <- "1,1,."
        A2 <- 1
      }

Specification of formula argument in the cSMART.mm function

How an end-user should approach specifying the formula argument in the cSMART.mm function might need to be revisited.

Coefficient of beta 2

Right now, it seems as if beta[2] is providing the effect estimate of the the effect of A_2 given non-response rather than just the effect of A_2. If you run 1000 simulations, then the mean effect estimate of beta[2] = true beta[2]/(1 - r) where r denotes the probability of response. I can't seem to figure out why this is occurring. The other three coefficients seem to be correct.
Code:

Results:

True Betas: