## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

run <- requireNamespace("dplyr", quietly = TRUE) 
knitr::opts_chunk$set(eval = run)

## ----trial design, echo=FALSE, fig.cap="SMART where only non-responders are re-randomized", fig.align="center", out.width="90%"----
knitr::include_graphics("four-arm-resp-g.png")

## ----load packages, message=FALSE, warning=FALSE------------------------------
library(rsmart)

## ----data---------------------------------------------------------------------
# This code block is used to generate the dataset used in the vignette. 
set.seed(1)
dat <- gen_no_trt_resp(n=300, s2=200, block_rep=2, r2p = 0.5)
regimes <- regime_list_no_trt_resp(emb_regimes = list(c(0,0), c(0,1), c(1,0), c(1,1)), 
                                  dat = dat, 
                                  resp_trt = list("r2" = list(0, 0, 0, 0)))

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

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

pi_list <- list(p1, 
                p2)

## ----regimes------------------------------------------------------------------
regime_all <- regime_list_no_trt_resp(emb_regimes = list(c(0,0), c(0,1), c(1,0), c(1,1)), 
                       dat, 
                       resp_trt=list("r2" = list(0, 0, 0, 0)))

## ----ipwe---------------------------------------------------------------------
ipweres <- iaipwe(df=dat, 
                  pi_list=pi_list, 
                  q_list=NULL, 
                  regime_all=regime_all, 
                  feasible_sets_indicator=TRUE, 
                  t_s=max(dat$t3)) 

## ----ipwe value---------------------------------------------------------------
ipweres$values
ipweres$se

## ----ipwe diff----------------------------------------------------------------
cont.mat <- matrix(data = c(1, 0, 0, -1, 
                0, 1, 0, -1, 
                0, 0, 1, -1), 
       nrow = 3, byrow=TRUE)
endi <- dim(ipweres$covariance)[2]
starti <- endi - length(ipweres$values) + 1
covvhat <- ipweres$covariance[starti:endi, starti:endi] / ipweres$nus$ns

chisq <- t(cont.mat %*% ipweres$values) %*% 
  (solve(cont.mat %*% covvhat %*% t(cont.mat))) %*% 
  (cont.mat %*% ipweres$values)
pchi <- 1-pchisq(q = chisq, df = 3)


## ----ipwe ci------------------------------------------------------------------
alpha <- 0.05
ci <- data.frame(
  regime = seq(1, 4),
  value = round(ipweres$values, 3),
  lower = round(ipweres$values - qnorm(1-alpha/2) * ipweres$se, 3),
  upper = round(ipweres$values + qnorm(1-alpha/2) * ipweres$se, 3)
)
ci

## ----q models-----------------------------------------------------------------
q2 <- modelObj::buildModelObj(model = ~ x11 + x12 + x21 +
                                a1 + a2 + a1:a2 + r2,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')
q1 <- modelObj::buildModelObj(model = ~ x11 + x12 + 
                                a1,
                              solver.method = 'lm',
                              predict.method = 'predict.lm')
q_list <- list(q1, q2)

## ----aipwe--------------------------------------------------------------------
aipweres <- iaipwe(df=dat, 
                  pi_list=pi_list, 
                  q_list=q_list, 
                  regime_all=regime_all, 
                  feasible_sets_indicator=TRUE, 
                  t_s=max(dat$t3)) 

## ----aipwe value--------------------------------------------------------------
aipweres$values
aipweres$se

## ----aipwe diff---------------------------------------------------------------
endi <- dim(aipweres$covariance)[2]
starti <- endi - length(aipweres$values) + 1
covvhat <- aipweres$covariance[starti:endi, starti:endi] / aipweres$nus$ns

chisq <- t(cont.mat %*% aipweres$values) %*% 
  (solve(cont.mat %*% covvhat %*% t(cont.mat))) %*% 
  (cont.mat %*% aipweres$values)
pchi <- 1-pchisq(q = chisq, df = 3)


## ----aipwe ci-----------------------------------------------------------------
ci <- data.frame(
  regime = seq(1, 4),
  value = round(aipweres$values, 3),
  lower = round(aipweres$values - qnorm(1-alpha/2) * aipweres$se, 3),
  upper = round(aipweres$values + qnorm(1-alpha/2) * aipweres$se, 3)
)
ci

## ----iaipwe-------------------------------------------------------------------
iaipweres <- iaipwe(df=dat, 
                  pi_list=pi_list, 
                  q_list=q_list, 
                  regime_all=regime_all, 
                  feasible_sets_indicator=TRUE, 
                  t_s=median(dat$t3))

## ----aipwe results------------------------------------------------------------
iaipweres$values
iaipweres$se

endi <- dim(iaipweres$covariance)[2]
starti <- endi - length(iaipweres$values) + 1
covvhat <- iaipweres$covariance[starti:endi, starti:endi] / iaipweres$nus$ns

chisq <- t(cont.mat %*% iaipweres$values) %*% 
  (solve(cont.mat %*% covvhat %*% t(cont.mat))) %*% 
  (cont.mat %*% iaipweres$values)
pchi <- 1-pchisq(q = chisq, df = 3)

ci <- data.frame(
  regime = seq(1, 4),
  value = round(iaipweres$values, 3),
  lower = round(iaipweres$values - qnorm(1-alpha/2) * iaipweres$se, 3),
  upper = round(iaipweres$values + qnorm(1-alpha/2) * iaipweres$se, 3)
)
ci