R/tsml.cara.rct.R
In achambaz/tsml.cara.rct: Targeted Sequential Minimum Loss CARA RCT Design and Inference
setMethodS3("update", "TSMLCARA", function(#Updates a TSMLCARA Object
### Updates a TSMLCARA object.
this,
### A \code{TSMLCARA} object, as created by \code{TSMLCARA}.
flavor=c("parametric", "lasso"),
### A \code{character} indicating the 'flavor' of the procedure.
...,
### Additional parameters.
verbose=FALSE
### A \code{logical} or an \code{integer} indicating the level of verbosity
### (defaults to 'FALSE').
) {
##alias<< update
##references<< Chambaz, van der Laan, Zheng, Chapter 16, Modern Adaptive Randomized Clinical Trials: Statistical, Operational, and Regulatory Aspects, by A. Sverdlov (CRC Press, 2015).
##seealso<<tsml.cara.rct, targetPsi
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## Validate arguments
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## Arguments 'flavor'
if (missing(flavor)) {
flavor <- getFlavor(this)
}
## Argument 'verbose'
verbose <- Arguments$getVerbose(verbose)
verbose <- less(verbose, 0)
## retrieving 'what'
what <- getWhat(this)
## retrieving 'obs', preparing 'A'
obs <- getObs(this)
A <- obs[, "A"]
## retrieving 'Qmin'
Qmin <- getQmin(this)
## retrieving 'tm.ref'
tm.ref <- getTm.ref(this)
## retrieving 'learnQ'
learnQ <- getLearnQ(this)
if (what=="ATE") {
## retrieving 'tm.model' and 'tm.control'
tm.model <- getTm.model(this)
tm.control <- getTm.control(this)
## retrieving 'targetLink'
targetLink <- getTargetLink(this)
}
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## learning
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
verbose && ruler(verbose)
verbose && enter(verbose, "UPDATING STEP")
## estimating 'Q'
verbose && enter(verbose, "Estimating Q")
W <- extractW(obs)
GA <- blik(A, obs[, "G"])
weights <- blik(A, tm.ref(W))
weights <- weights/GA
Q <- estimateQ(obs, learnQ=learnQ, weights=weights, Qmin=Qmin,
flavor=flavor, ..., verbose=verbose)
Qtab <- matrix(NA, nrow=nrow(obs), ncol=3)
colnames(Qtab) <- c("A", "A=0", "A=1")
Qtab[, "A=1"] <- Q(1, W)
Qtab[, "A=0"] <- Q(0, W)
Aone <- (A==1)
Qtab[Aone, "A"] <- Qtab[Aone, "A=1"]
Qtab[!Aone, "A"] <- Qtab[!Aone, "A=0"]
setQtab(this, Qtab)
setLearnedQ(this, attr(Q, "model"))
verbose && str(verbose, Qtab)
verbose && exit(verbose)
## targeting 'Gstar'
verbose && enter(verbose, "Targeting Gstar")
if (what=="ATE") {
weights <- (obs[, "Y"]-Qtab[, "A"])^2/GA
Gstar <- targetGstar(obs=obs, weights=weights, tm.model=tm.model,
targetLink=targetLink, tm.control=tm.control,
...,
verbose=verbose)
} else if (what=="MOR") {
Gstar <- targetOptRule(this, Q, ..., verbose=verbose)
}
setGstar(this, Gstar)
verbose && str(verbose, Gstar)
verbose && exit(verbose)
## updating 'Gtab'
verbose && enter(verbose, "Updating Gtab")
setGtab(this, Gstar(W))
verbose && str(verbose, Gstar)
verbose && exit(verbose)
## updating 'history'
verbose && enter(verbose, "Updating history")
step <- getStep(this)
step["update"] <- step["update"]+1
this$.step <- step
updateHistory(this, "update")
verbose && str(verbose, tsml.cara.rct::getHistory.TSMLCARA(this))
verbose && exit(verbose)
verbose && exit(verbose)
verbose && ruler(verbose)
})
setMethodS3("targetPsi", "TSMLCARA", function(#Targets a TSMLCARA Object Toward the Parameter Psi
### Targets a TSMLCARA object toward the parameter Psi.
this,
### A \code{TSMLCARA} object, as created by \code{TSMLCARA}.
...,
### Additional parameters.
verbose=FALSE
### A \code{logical} or an \code{integer} indicating the level of verbosity
### (defaults to 'FALSE').
) {
##alias<< targetPsi
##references<< Chambaz, van der Laan, Zheng, Chapter 16, Modern Adaptive Randomized Clinical Trials: Statistical, Operational, and Regulatory Aspects, by A. Sverdlov (CRC Press, 2015).
##seealso<< tsml.cara.rct, update
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## Validate arguments
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## Argument 'verbose'
verbose <- Arguments$getVerbose(verbose)
verbose <- less(verbose, 0)
## retrieving 'what'
what <- getWhat(this)
if (what=="ATE") {
targetLink <- NULL
}
else if (what=="MOR") {
targetLink <- getTargetLink(this)
}
## retrieving 'obs', 'Gtab', 'Qtab', preparing 'Y', 'A', 'G' and 'GA'
obs <- getObs(this)
Y <- obs[, "Y"]
A <- obs[, "A"]
G <- obs[, "G"]
GA <- blik(A, G)
GstarW <- getGtab(this)
Qtab <- getQtab(this)
## retrieving 'Qmin'
Qmin <- getQmin(this)
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## targeting Psi
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
verbose && ruler(verbose)
verbose && enter(verbose, "TARGETING PSI")
W <- extractW(obs)
GAstarW <- blik(A, GstarW)
weights <- GAstarW/GA
verbose && enter(verbose, "Fluctuating")
eps <- fluctuateQ(obs=obs, Qtab=Qtab, GAstarW=GAstarW,
what=what, weights=weights, Qmin=Qmin,
targetLink=targetLink,
..., verbose)
verbose && str(verbose, eps)
## QEpstab
epsHtab <- matrix(NA, nrow=nrow(obs), ncol=3)
colnames(epsHtab) <- colnames(Qtab)
if (what=="ATE") {
epsHtab[, "A=1"] <- eps/GstarW
epsHtab[, "A=0"] <- -eps/(1-GstarW)
} else if (what=="MOR") {
rA <- ((Qtab[, "A=1"]-Qtab[, "A=0"])>0)
epsHtab[, "A=1"] <- eps/GstarW
epsHtab[!rA, "A=1"] <- 0
epsHtab[, "A=0"] <- eps/(1-GstarW)
epsHtab[rA, "A=0"] <- 0
}
Aone <- (A==1)
epsHtab[Aone, "A"] <- epsHtab[Aone, "A=1"]
epsHtab[!Aone, "A"] <- epsHtab[!Aone, "A=0"]
scaledQtab <- scaleQmat(Qtab, wrt=Y, thr=Qmin)
if (what=="ATE") {
scaledQEpstab <- binomial()$linkinv(binomial()$linkfun(scaledQtab) + epsHtab) ##ie, expit(logit(...))
} else if (what=="MOR") {
scaledQEpstab <- targetLink$linkinv(targetLink$linkfun(scaledQtab) + epsHtab)
}
QEpstab <- scaleQmat(scaledQEpstab, reverse=TRUE)
attr(QEpstab, "eps") <- eps
setQEpstab(this, QEpstab)
verbose && str(verbose, QEpstab)
verbose && exit(verbose)
## estimating Psi
verbose && enter(verbose, "estimating psi")
## CAUTION: same form in both cases! no need to fork
if (what=="ATE") {
## epsHtab <- epsHtab/abs(eps) * matrix(weights, nrow=length(weights), ncol=3)
epsHtab <- epsHtab/eps * matrix(weights, nrow=length(weights), ncol=3)
psis <- estimatePsi(obs=obs, what=what, Qtab=NULL, QEpstab=QEpstab, epsHtab=epsHtab, verbose=verbose)
} else if (what=="MOR") {
epsHtab <- epsHtab/eps * matrix(weights, nrow=length(weights), ncol=3)
psis <- estimatePsi(obs=obs, what=what, Qtab=Qtab, QEpstab=QEpstab, epsHtab=epsHtab, verbose=verbose)
}
this$.psi <- psis["psi"]
this$.psi.sd <- psis["psi.sd"]
verbose && str(verbose, psis)
verbose && exit(verbose)
if (what=="MOR") {
## if targeting mean under optimal rule, then estimating the empirical
## regret as well
verbose && enter(verbose, "estimating the regret")
regret <- estimateRegret(obs=obs, what=what, Qtab=Qtab, QEpstab=QEpstab,
epsHtab=epsHtab, psi=psis["psi"], verbose=verbose)
this$.regret <- regret["regret"]
this$.regret.sd <- regret["regret.sd"]
verbose && str(verbose, regret)
verbose && exit(verbose)
}
## updating 'history'
verbose && enter(verbose, "Updating history")
step <- getStep(this)
step["target"] <- step["target"]+1
this$.step <- step
updateHistory(this, "target")
verbose && str(verbose, tsml.cara.rct::getHistory.TSMLCARA(this))
verbose && exit(verbose)
verbose && exit(verbose)
verbose && ruler(verbose)
})
tsml.cara.rct <- structure(function#Targeted Minimum Loss Covariate-Adjusted Response-Adaptive RCT Design and Statistical Analysis
### Simulates a targeted minimum loss covariate-adjusted response-adaptive RCT
### design and statistical analysis.
(what=c("ATE", "MOR"),
### A \code{character} indicating the parameter of interest to estimate.
### Either "ATE" for the Average Treatment Effect, the difference between the
### means under '\eqn{do(A=1)}' and '\eqn{do(A=0)}', or "MOR" for the Mean
### under the Optimal treatment Rule '\eqn{do(A=r(W))}'.
flavor=c("parametric", "lasso"),
### A \code{character} indicating the 'flavor' of the procedure.
ninit=50,
### An \code{integer}, number of subjects to sample at
### initialization. Defaults to 50.
by=25,
### An \code{integer}, number of subjects to sample at each step after
### initialization. Defaults to 25.
nmax=500,
### An \code{integer}, maximum number of subjects to sample during the
### trial. Must be larger than 'ninit+by'. Defaults to 500.
tm.init=oneOne,
### A \code{function} describing the initial treatment mechanism to be
### employed. Defaults to the balanced (1:1) treatment mechanism, ie,
### \code{function} \code{\link{oneOne}}.
tm.ref=oneOne,
### A \code{function} describing the reference treatment mechanism to be
### employed. Defaults to the balanced (1:1) treatment mechanism, ie,
### \code{function} \code{\link{oneOne}}.
learnQ,
### A model \eqn{{\cal Q}} of conditional expectations of \eqn{Y} given
### \eqn{(A,W)} for both flavors 'parametric' and 'lasso', given as a
### \code{formula} or a \code{function} outputing formulas. Defaults to
### \code{formula(Y~1)} for flavors 'parametric' and 'lasso'.
tm.model=formula(A~1),
### A parametric model \eqn{{\cal G}} of treatment mechanisms, used only when
### 'what' equals "ATE". The procedure targets the optimal treatment
### mechanism within this model. Defaults to \code{formula(A~1)}.
tm.control=glm.control(maxit=500),
### A \code{list} of options for the targeting of the treatment mechanism
### within the model defined by argument 'tm.model'. Used only when 'what'
### equals "ATE", it defaults to \code{glm.control(maxit=500)}.
Gmin=1e-2,
### A small positive \code{numeric}, with default value \code{1e-2}. When
### \code{what} equals 'ATE', it is the minimum value of elements of the
### parametric model \eqn{{\cal G}} of treatment mechanisms (see argument
### \code{tm.model}). The maximum value is \code{1-Gmin}. When \code{what}
### equals 'MOR', it is the minimum value of the conditional probability
### of \eqn{A=r_n(W)} given \eqn{W}.
Gexploit=Gmin,
### A small positive \code{numeric}, with default value that of \code{Gmin},
### or a function of sample size giving such small numbers, only used when
### \code{what} equals "MOR", in conjunction with \code{Gexplore}.
Gexplore=1e-2,
### Either a small positive \code{numeric}, with default value \code{1e-2}, or
### a function of sample size giving such small numbers, only used when
### \code{what} equals "MOR", in conjunction with \code{Gexploit}.
Qmin=1e-2,
### A small positive \code{numeric}, the minimum value of scaled outcomes
### \eqn{Y}. The maximum value is \code{1-Qmin}.
conf.level=0.95,
### A \code{numeric}, the confidence level of the resulting confidence
### interval.
verbose=FALSE,
### A \code{logical} or an \code{integer} indicating the level of verbosity
### (defaults to 'FALSE').
piV=c(1/2, 1/3, 1/6),
### Marginal distribution of \eqn{V}. Defaults to \code{c(1/2, 1/3, 1/6)}.
family=c("beta", "gamma"),
### A \code{character}, either "beta" (default) or "gamma", the nature of the
### law of outcome.
Qbar=Qbar1,
### A \code{function}, the conditional expectation of \eqn{Y} given
### \eqn{(A,W)}. Defaults to \code{Qbar1}.
Vbar=Vbar1,
### A \code{function}, the conditional variance of \eqn{Y} given
### \eqn{(A,W)}. Defaults to \code{Vbar1}.
Bn=1e5,
### An \code{integer}, the sample size used to estimate the true value of the
### data-adaptive parameter at each step of the procedure when 'what' equals
### 'MOR'. Defaults to 1e5.
slice.by=1e5
### An \code{integer}. If it is smaller than argument 'n' of 'getSample', then
### the simulation is decomposed into 'n%/%slice.by' smaller simulations of
### 'slice.by' observations and one of 'n%%slice.by' observations. Defaults to
### 1e5 (hence, no decomposition if 'n' smaller than 4e5). Mainly for internal
### use.
) {
##alias<< tsml.cara.rct
##references<< Chambaz, van der Laan, Zheng, Chapter 16, Modern Adaptive Randomized Clinical Trials: Statistical, Operational, and Regulatory Aspects, by A. Sverdlov (CRC Press, 2015).
##seealso<< update, targetPsi, getSample
##details<< Defines a lower-bound on the conditional probability of
## \eqn{do(A=1-r_n(W))} given \eqn{W}.
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## Validate arguments
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## Argument 'what'
what <- match.arg(what)
## Arguments 'flavor' and 'learnQ'
flavor <- match.arg(flavor)
if (missing(learnQ)) {
learnQ <- switch(flavor,
parametric=formula(Y~1),
lasso=formula(Y~1))
} else {
learnMode <- switch(flavor,
parametric=c("formula", "function"),
lasso=c("formula", "function"))
class <- class(learnQ)
if (!class%in%learnMode) {
throw("Argument 'learnQ' should be of class '", learnMode, "', not '", class, "' for flavor: ", flavor)
}
if (class=="formula") {
if (learnQ[[2]]!="Y") {
throw("Argument 'learnQ' should be a formula with response 'Y', not:", learnQ)
}
}
}
## Argument 'ninit'
ninit <- Arguments$getNumeric(ninit)
## Argument 'by'
by <- Arguments$getNumeric(by)
## Argument 'nmax'
nmax <- Arguments$getNumeric(nmax, c(ninit+by, Inf))
## Argument 'tm.init'
mode <- mode(tm.init)
if (mode != "function") {
throw("Argument 'tm.init' should be of mode 'function', not '", mode)
}
## Argument 'tm.ref'
mode <- mode(tm.ref)
if (mode != "function") {
throw("Argument 'tm.ref' should be of mode 'function', not '", mode)
}
## Argument 'tm.model'
if (!identical(class(tm.model), "formula") | tm.model[[2]]!="A") {
throw("Argument 'tm.model' should be a formula with response 'A', not:",
deparse(substitute(form, list(form=tm.model))))
}
## Argument 'tm.control'
## Argument 'Gmin'
Gmin <- Arguments$getNumeric(Gmin, c(0, 1/2))
## Argument 'Gexploit'
mode <- mode(Gexploit)
if (!(mode %in% c("numeric", "function"))) {
throw("Argument 'Gexploit' should be of mode either 'numeric' or 'function', not ", mode)
} else if (mode=="numeric") {
Gexploit <- Arguments$getNumeric(Gexploit, c(0, 1/2))
}
## Argument 'Gexplore'
mode <- mode(Gexplore)
if (!(mode %in% c("numeric", "function"))) {
throw("Argument 'Gexplore' should be of mode either 'numeric' or 'function', not ", mode)
} else if (mode=="numeric") {
Gexplore <- Arguments$getNumeric(Gexplore, c(0, 1/2))
}
## Argument 'verbose'
verbose <- Arguments$getVerbose(verbose)
verbose <- less(verbose, 10)
## Argument 'piV':
piV <- Arguments$getNumerics(piV, range=c(0,1))
if (sum(piV)!=1) {
throw("Argument 'piV' should consist of non-negative weights summing to one.")
}
## Argument 'family':
family <- match.arg(family)
## Argument 'Qbar':
mode <- mode(Qbar);
if (mode != "function") {
throw("Argument 'Qbar' should be of mode 'function', not '", mode)
}
## Argument 'Vbar':
mode <- mode(Vbar);
if (mode != "function") {
throw("Argument 'Vbar' should be of mode 'function', not '", mode)
}
## Argument 'Bn':
Bn <- Arguments$getInteger(Bn);
if (Bn <= 1e5) {
verbose && str(verbose, "isn't 'Bn' too small?")
}
## Argument 'slice.by'
slice.by <- Arguments$getInteger(slice.by, c(1, Inf));
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## Core
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## initial sampling (under 1:1 treatment mechanism)
## (a) sampling
obs <- getSample(ninit, tm=tm.ref,
piV=piV, Qbar=Qbar, Vbar=Vbar,
family=family)
nobs <- nrow(obs)
## (b) declaring the 'TSMLCARA' object
tsml.cara <- TSMLCARA(what=what,
flavor=flavor,
obs=obs,
tm.ref=tm.ref,
tm.model=tm.model,
tm.control=tm.control,
learnQ=learnQ,
Gmin=Gmin,
Gexploit=Gexploit,
Gexplore=Gexplore,
Qmin=Qmin,
conf.level=conf.level,
verbose=verbose)
## update of 'tsml.cara'
update(tsml.cara, verbose=verbose)
targetPsi(tsml.cara, verbose=verbose)
## if 'what=="MOR"', then compute and store the targeted data-adaptive parameters
if (what=="MOR") {
Qn <- getQ(tsml.cara)
ruleQn <- ruleQbar(Qn)
##
## first data-adaptive parameter:
## mean reward under the current best estimate of the optimal rule
##
truthn <- getSample(Bn, tm=oneOne, rule=ruleQn,
piV=piV, Qbar=Qbar, Vbar=Vbar, what,
family=family, slice.by=slice.by)[c("psi", "psi.sd")]
##
## second data-adaptive parameter:
## empirical cumulated regret
##
W <- extractW(obs)
rA <- as.numeric(ruleQn(W))
eregretn <- mean(obs[, "Y"] - Qbar(cbind(A=rA, W)))
##
## third data-adaptive parameter:
## counterfactual cumulated regret
##
col <- match(c("Y0", "Y1"), colnames(obs))
rA <- col[rA+1]
counterfactual <- as.numeric(obs[cbind(1:nrow(obs), rA)])
cregretn <- mean(obs[, "Y"] - counterfactual)
}
## successive updates of 'tsml.cara'
while (nobs < nmax) {
if (nobs+by>nmax) {
by <- (nmax-nobs)
}
newObs <- getSample(by, tm=getGstar(tsml.cara),
piV=piV, Qbar=Qbar, Vbar=Vbar,
family=family)
addNewSample(tsml.cara, newObs)
nobs <- nrow(getObs(tsml.cara))
update(tsml.cara, verbose=verbose)
targetPsi(tsml.cara, verbose=verbose)
if (what=="MOR") {
Qn <- getQ(tsml.cara)
ruleQn <- ruleQbar(Qn)
##
## first data-adaptive parameter:
## mean reward under the current best estimate of the optimal rule
##
truthn <- rbind(truthn,
getSample(Bn, tm=oneOne, rule=ruleQn,
piV=piV, Qbar=Qbar, Vbar=Vbar, what,
family=family, slice.by=slice.by)[c("psi", "psi.sd")])
##
## second data-adaptive parameter:
## empirical cumulated regret
##
obs <- getObs(tsml.cara)
W <- extractW(obs)
rA <- as.numeric(ruleQn(W))
eregretn <- c(eregretn,
mean(obs[, "Y"] - Qbar(cbind(A=rA, W))))
##
## third data-adaptive parameter:
## counterfactual cumulated regret
##
rA <- col[rA+1]
counterfactual <- as.numeric(obs[cbind(1:nrow(obs), rA)])
cregretn <- c(cregretn,
mean(obs[, "Y"] - counterfactual))
}
}
if (what=="MOR") {
truthn[, "psi.sd"] <- truthn[, "psi.sd"]/sqrt(Bn)
colnames(truthn) <- c("psin", "psin.sd")
rownames(truthn) <- NULL
attr(tsml.cara, "truthn") <- truthn
attr(tsml.cara, "eregretn") <- eregretn
attr(tsml.cara, "cregretn") <- cregretn
}
return(tsml.cara)
### Returns a \code{TSMLCARA} object which summarizes the TSMLCARA undertaken
### procedure.
}, ex=function() {
##
log <- Arguments$getVerbose(-1, timestamp=TRUE)
set.seed(12345)
## ########################
## AVERAGE TREATMENT EFFECT
## ########################
tm.model <- formula(A~.)
psi.sd <- sqrt(getOptVar(n=1e5,
tm.model=tm.model,
piV=c(1/2, 1/3, 1/6),
family="gamma",
Qbar=Qbar1,
Vbar=Vbar1))
truth <- c(psi=91/72, psi.sd=psi.sd)
## parametric example
learnQ <- formula(Y~I(as.integer(A)):(U+V)+I(as.integer(1-A)):(U+V))
ATE.param <- tsml.cara.rct(what="ATE",
flavor="parametric",
ninit=200,
by=100,
nmax=400,
tm.init=oneOne,
tm.ref=oneOne,
learnQ=learnQ,
tm.model=tm.model,
conf.level=0.95,
piV=c(1/2, 1/3, 1/6),
family="gamma",
Qbar=Qbar1,
Vbar=Vbar1)
ATE.param
## Not run:
plot(ATE.param, truth=truth)
## End(**Not run**)
## See the vignette for more examples...
})
############################################################################
## HISTORY:
## 2016-09-16
## o Created.
############################################################################
achambaz/tsml.cara.rct documentation built on May 10, 2019, 5:10 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
setMethodS3("update", "TSMLCARA", function(#Updates a TSMLCARA Object
### Updates a TSMLCARA object.
this,
### A \code{TSMLCARA} object, as created by \code{TSMLCARA}.
flavor=c("parametric", "lasso"),
### A \code{character} indicating the 'flavor' of the procedure.
...,
### Additional parameters.
verbose=FALSE
### A \code{logical} or an \code{integer} indicating the level of verbosity
### (defaults to 'FALSE').
) {
##alias<< update
##references<< Chambaz, van der Laan, Zheng, Chapter 16, Modern Adaptive Randomized Clinical Trials: Statistical, Operational, and Regulatory Aspects, by A. Sverdlov (CRC Press, 2015).
##seealso<<tsml.cara.rct, targetPsi
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## Validate arguments
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## Arguments 'flavor'
if (missing(flavor)) {
flavor <- getFlavor(this)
}
## Argument 'verbose'
verbose <- Arguments$getVerbose(verbose)
verbose <- less(verbose, 0)
## retrieving 'what'
what <- getWhat(this)
## retrieving 'obs', preparing 'A'
obs <- getObs(this)
A <- obs[, "A"]
## retrieving 'Qmin'
Qmin <- getQmin(this)
## retrieving 'tm.ref'
tm.ref <- getTm.ref(this)
## retrieving 'learnQ'
learnQ <- getLearnQ(this)
if (what=="ATE") {
## retrieving 'tm.model' and 'tm.control'
tm.model <- getTm.model(this)
tm.control <- getTm.control(this)
## retrieving 'targetLink'
targetLink <- getTargetLink(this)
}
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## learning
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
verbose && ruler(verbose)
verbose && enter(verbose, "UPDATING STEP")
## estimating 'Q'
verbose && enter(verbose, "Estimating Q")
W <- extractW(obs)
GA <- blik(A, obs[, "G"])
weights <- blik(A, tm.ref(W))
weights <- weights/GA
Q <- estimateQ(obs, learnQ=learnQ, weights=weights, Qmin=Qmin,
flavor=flavor, ..., verbose=verbose)
Qtab <- matrix(NA, nrow=nrow(obs), ncol=3)
colnames(Qtab) <- c("A", "A=0", "A=1")
Qtab[, "A=1"] <- Q(1, W)
Qtab[, "A=0"] <- Q(0, W)
Aone <- (A==1)
Qtab[Aone, "A"] <- Qtab[Aone, "A=1"]
Qtab[!Aone, "A"] <- Qtab[!Aone, "A=0"]
setQtab(this, Qtab)
setLearnedQ(this, attr(Q, "model"))
verbose && str(verbose, Qtab)
verbose && exit(verbose)
## targeting 'Gstar'
verbose && enter(verbose, "Targeting Gstar")
if (what=="ATE") {
weights <- (obs[, "Y"]-Qtab[, "A"])^2/GA
Gstar <- targetGstar(obs=obs, weights=weights, tm.model=tm.model,
targetLink=targetLink, tm.control=tm.control,
...,
verbose=verbose)
} else if (what=="MOR") {
Gstar <- targetOptRule(this, Q, ..., verbose=verbose)
}
setGstar(this, Gstar)
verbose && str(verbose, Gstar)
verbose && exit(verbose)
## updating 'Gtab'
verbose && enter(verbose, "Updating Gtab")
setGtab(this, Gstar(W))
verbose && str(verbose, Gstar)
verbose && exit(verbose)
## updating 'history'
verbose && enter(verbose, "Updating history")
step <- getStep(this)
step["update"] <- step["update"]+1
this$.step <- step
updateHistory(this, "update")
verbose && str(verbose, tsml.cara.rct::getHistory.TSMLCARA(this))
verbose && exit(verbose)
verbose && exit(verbose)
verbose && ruler(verbose)
})
setMethodS3("targetPsi", "TSMLCARA", function(#Targets a TSMLCARA Object Toward the Parameter Psi
### Targets a TSMLCARA object toward the parameter Psi.
this,
### A \code{TSMLCARA} object, as created by \code{TSMLCARA}.
...,
### Additional parameters.
verbose=FALSE
### A \code{logical} or an \code{integer} indicating the level of verbosity
### (defaults to 'FALSE').
) {
##alias<< targetPsi
##references<< Chambaz, van der Laan, Zheng, Chapter 16, Modern Adaptive Randomized Clinical Trials: Statistical, Operational, and Regulatory Aspects, by A. Sverdlov (CRC Press, 2015).
##seealso<< tsml.cara.rct, update
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## Validate arguments
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## Argument 'verbose'
verbose <- Arguments$getVerbose(verbose)
verbose <- less(verbose, 0)
## retrieving 'what'
what <- getWhat(this)
if (what=="ATE") {
targetLink <- NULL
}
else if (what=="MOR") {
targetLink <- getTargetLink(this)
}
## retrieving 'obs', 'Gtab', 'Qtab', preparing 'Y', 'A', 'G' and 'GA'
obs <- getObs(this)
Y <- obs[, "Y"]
A <- obs[, "A"]
G <- obs[, "G"]
GA <- blik(A, G)
GstarW <- getGtab(this)
Qtab <- getQtab(this)
## retrieving 'Qmin'
Qmin <- getQmin(this)
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## targeting Psi
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
verbose && ruler(verbose)
verbose && enter(verbose, "TARGETING PSI")
W <- extractW(obs)
GAstarW <- blik(A, GstarW)
weights <- GAstarW/GA
verbose && enter(verbose, "Fluctuating")
eps <- fluctuateQ(obs=obs, Qtab=Qtab, GAstarW=GAstarW,
what=what, weights=weights, Qmin=Qmin,
targetLink=targetLink,
..., verbose)
verbose && str(verbose, eps)
## QEpstab
epsHtab <- matrix(NA, nrow=nrow(obs), ncol=3)
colnames(epsHtab) <- colnames(Qtab)
if (what=="ATE") {
epsHtab[, "A=1"] <- eps/GstarW
epsHtab[, "A=0"] <- -eps/(1-GstarW)
} else if (what=="MOR") {
rA <- ((Qtab[, "A=1"]-Qtab[, "A=0"])>0)
epsHtab[, "A=1"] <- eps/GstarW
epsHtab[!rA, "A=1"] <- 0
epsHtab[, "A=0"] <- eps/(1-GstarW)
epsHtab[rA, "A=0"] <- 0
}
Aone <- (A==1)
epsHtab[Aone, "A"] <- epsHtab[Aone, "A=1"]
epsHtab[!Aone, "A"] <- epsHtab[!Aone, "A=0"]
scaledQtab <- scaleQmat(Qtab, wrt=Y, thr=Qmin)
if (what=="ATE") {
scaledQEpstab <- binomial()$linkinv(binomial()$linkfun(scaledQtab) + epsHtab) ##ie, expit(logit(...))
} else if (what=="MOR") {
scaledQEpstab <- targetLink$linkinv(targetLink$linkfun(scaledQtab) + epsHtab)
}
QEpstab <- scaleQmat(scaledQEpstab, reverse=TRUE)
attr(QEpstab, "eps") <- eps
setQEpstab(this, QEpstab)
verbose && str(verbose, QEpstab)
verbose && exit(verbose)
## estimating Psi
verbose && enter(verbose, "estimating psi")
## CAUTION: same form in both cases! no need to fork
if (what=="ATE") {
## epsHtab <- epsHtab/abs(eps) * matrix(weights, nrow=length(weights), ncol=3)
epsHtab <- epsHtab/eps * matrix(weights, nrow=length(weights), ncol=3)
psis <- estimatePsi(obs=obs, what=what, Qtab=NULL, QEpstab=QEpstab, epsHtab=epsHtab, verbose=verbose)
} else if (what=="MOR") {
epsHtab <- epsHtab/eps * matrix(weights, nrow=length(weights), ncol=3)
psis <- estimatePsi(obs=obs, what=what, Qtab=Qtab, QEpstab=QEpstab, epsHtab=epsHtab, verbose=verbose)
}
this$.psi <- psis["psi"]
this$.psi.sd <- psis["psi.sd"]
verbose && str(verbose, psis)
verbose && exit(verbose)
if (what=="MOR") {
## if targeting mean under optimal rule, then estimating the empirical
## regret as well
verbose && enter(verbose, "estimating the regret")
regret <- estimateRegret(obs=obs, what=what, Qtab=Qtab, QEpstab=QEpstab,
epsHtab=epsHtab, psi=psis["psi"], verbose=verbose)
this$.regret <- regret["regret"]
this$.regret.sd <- regret["regret.sd"]
verbose && str(verbose, regret)
verbose && exit(verbose)
}
## updating 'history'
verbose && enter(verbose, "Updating history")
step <- getStep(this)
step["target"] <- step["target"]+1
this$.step <- step
updateHistory(this, "target")
verbose && str(verbose, tsml.cara.rct::getHistory.TSMLCARA(this))
verbose && exit(verbose)
verbose && exit(verbose)
verbose && ruler(verbose)
})
tsml.cara.rct <- structure(function#Targeted Minimum Loss Covariate-Adjusted Response-Adaptive RCT Design and Statistical Analysis
### Simulates a targeted minimum loss covariate-adjusted response-adaptive RCT
### design and statistical analysis.
(what=c("ATE", "MOR"),
### A \code{character} indicating the parameter of interest to estimate.
### Either "ATE" for the Average Treatment Effect, the difference between the
### means under '\eqn{do(A=1)}' and '\eqn{do(A=0)}', or "MOR" for the Mean
### under the Optimal treatment Rule '\eqn{do(A=r(W))}'.
flavor=c("parametric", "lasso"),
### A \code{character} indicating the 'flavor' of the procedure.
ninit=50,
### An \code{integer}, number of subjects to sample at
### initialization. Defaults to 50.
by=25,
### An \code{integer}, number of subjects to sample at each step after
### initialization. Defaults to 25.
nmax=500,
### An \code{integer}, maximum number of subjects to sample during the
### trial. Must be larger than 'ninit+by'. Defaults to 500.
tm.init=oneOne,
### A \code{function} describing the initial treatment mechanism to be
### employed. Defaults to the balanced (1:1) treatment mechanism, ie,
### \code{function} \code{\link{oneOne}}.
tm.ref=oneOne,
### A \code{function} describing the reference treatment mechanism to be
### employed. Defaults to the balanced (1:1) treatment mechanism, ie,
### \code{function} \code{\link{oneOne}}.
learnQ,
### A model \eqn{{\cal Q}} of conditional expectations of \eqn{Y} given
### \eqn{(A,W)} for both flavors 'parametric' and 'lasso', given as a
### \code{formula} or a \code{function} outputing formulas. Defaults to
### \code{formula(Y~1)} for flavors 'parametric' and 'lasso'.
tm.model=formula(A~1),
### A parametric model \eqn{{\cal G}} of treatment mechanisms, used only when
### 'what' equals "ATE". The procedure targets the optimal treatment
### mechanism within this model. Defaults to \code{formula(A~1)}.
tm.control=glm.control(maxit=500),
### A \code{list} of options for the targeting of the treatment mechanism
### within the model defined by argument 'tm.model'. Used only when 'what'
### equals "ATE", it defaults to \code{glm.control(maxit=500)}.
Gmin=1e-2,
### A small positive \code{numeric}, with default value \code{1e-2}. When
### \code{what} equals 'ATE', it is the minimum value of elements of the
### parametric model \eqn{{\cal G}} of treatment mechanisms (see argument
### \code{tm.model}). The maximum value is \code{1-Gmin}. When \code{what}
### equals 'MOR', it is the minimum value of the conditional probability
### of \eqn{A=r_n(W)} given \eqn{W}.
Gexploit=Gmin,
### A small positive \code{numeric}, with default value that of \code{Gmin},
### or a function of sample size giving such small numbers, only used when
### \code{what} equals "MOR", in conjunction with \code{Gexplore}.
Gexplore=1e-2,
### Either a small positive \code{numeric}, with default value \code{1e-2}, or
### a function of sample size giving such small numbers, only used when
### \code{what} equals "MOR", in conjunction with \code{Gexploit}.
Qmin=1e-2,
### A small positive \code{numeric}, the minimum value of scaled outcomes
### \eqn{Y}. The maximum value is \code{1-Qmin}.
conf.level=0.95,
### A \code{numeric}, the confidence level of the resulting confidence
### interval.
verbose=FALSE,
### A \code{logical} or an \code{integer} indicating the level of verbosity
### (defaults to 'FALSE').
piV=c(1/2, 1/3, 1/6),
### Marginal distribution of \eqn{V}. Defaults to \code{c(1/2, 1/3, 1/6)}.
family=c("beta", "gamma"),
### A \code{character}, either "beta" (default) or "gamma", the nature of the
### law of outcome.
Qbar=Qbar1,
### A \code{function}, the conditional expectation of \eqn{Y} given
### \eqn{(A,W)}. Defaults to \code{Qbar1}.
Vbar=Vbar1,
### A \code{function}, the conditional variance of \eqn{Y} given
### \eqn{(A,W)}. Defaults to \code{Vbar1}.
Bn=1e5,
### An \code{integer}, the sample size used to estimate the true value of the
### data-adaptive parameter at each step of the procedure when 'what' equals
### 'MOR'. Defaults to 1e5.
slice.by=1e5
### An \code{integer}. If it is smaller than argument 'n' of 'getSample', then
### the simulation is decomposed into 'n%/%slice.by' smaller simulations of
### 'slice.by' observations and one of 'n%%slice.by' observations. Defaults to
### 1e5 (hence, no decomposition if 'n' smaller than 4e5). Mainly for internal
### use.
) {
##alias<< tsml.cara.rct
##references<< Chambaz, van der Laan, Zheng, Chapter 16, Modern Adaptive Randomized Clinical Trials: Statistical, Operational, and Regulatory Aspects, by A. Sverdlov (CRC Press, 2015).
##seealso<< update, targetPsi, getSample
##details<< Defines a lower-bound on the conditional probability of
## \eqn{do(A=1-r_n(W))} given \eqn{W}.
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## Validate arguments
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## Argument 'what'
what <- match.arg(what)
## Arguments 'flavor' and 'learnQ'
flavor <- match.arg(flavor)
if (missing(learnQ)) {
learnQ <- switch(flavor,
parametric=formula(Y~1),
lasso=formula(Y~1))
} else {
learnMode <- switch(flavor,
parametric=c("formula", "function"),
lasso=c("formula", "function"))
class <- class(learnQ)
if (!class%in%learnMode) {
throw("Argument 'learnQ' should be of class '", learnMode, "', not '", class, "' for flavor: ", flavor)
}
if (class=="formula") {
if (learnQ[[2]]!="Y") {
throw("Argument 'learnQ' should be a formula with response 'Y', not:", learnQ)
}
}
}
## Argument 'ninit'
ninit <- Arguments$getNumeric(ninit)
## Argument 'by'
by <- Arguments$getNumeric(by)
## Argument 'nmax'
nmax <- Arguments$getNumeric(nmax, c(ninit+by, Inf))
## Argument 'tm.init'
mode <- mode(tm.init)
if (mode != "function") {
throw("Argument 'tm.init' should be of mode 'function', not '", mode)
}
## Argument 'tm.ref'
mode <- mode(tm.ref)
if (mode != "function") {
throw("Argument 'tm.ref' should be of mode 'function', not '", mode)
}
## Argument 'tm.model'
if (!identical(class(tm.model), "formula") | tm.model[[2]]!="A") {
throw("Argument 'tm.model' should be a formula with response 'A', not:",
deparse(substitute(form, list(form=tm.model))))
}
## Argument 'tm.control'
## Argument 'Gmin'
Gmin <- Arguments$getNumeric(Gmin, c(0, 1/2))
## Argument 'Gexploit'
mode <- mode(Gexploit)
if (!(mode %in% c("numeric", "function"))) {
throw("Argument 'Gexploit' should be of mode either 'numeric' or 'function', not ", mode)
} else if (mode=="numeric") {
Gexploit <- Arguments$getNumeric(Gexploit, c(0, 1/2))
}
## Argument 'Gexplore'
mode <- mode(Gexplore)
if (!(mode %in% c("numeric", "function"))) {
throw("Argument 'Gexplore' should be of mode either 'numeric' or 'function', not ", mode)
} else if (mode=="numeric") {
Gexplore <- Arguments$getNumeric(Gexplore, c(0, 1/2))
}
## Argument 'verbose'
verbose <- Arguments$getVerbose(verbose)
verbose <- less(verbose, 10)
## Argument 'piV':
piV <- Arguments$getNumerics(piV, range=c(0,1))
if (sum(piV)!=1) {
throw("Argument 'piV' should consist of non-negative weights summing to one.")
}
## Argument 'family':
family <- match.arg(family)
## Argument 'Qbar':
mode <- mode(Qbar);
if (mode != "function") {
throw("Argument 'Qbar' should be of mode 'function', not '", mode)
}
## Argument 'Vbar':
mode <- mode(Vbar);
if (mode != "function") {
throw("Argument 'Vbar' should be of mode 'function', not '", mode)
}
## Argument 'Bn':
Bn <- Arguments$getInteger(Bn);
if (Bn <= 1e5) {
verbose && str(verbose, "isn't 'Bn' too small?")
}
## Argument 'slice.by'
slice.by <- Arguments$getInteger(slice.by, c(1, Inf));
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## Core
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## initial sampling (under 1:1 treatment mechanism)
## (a) sampling
obs <- getSample(ninit, tm=tm.ref,
piV=piV, Qbar=Qbar, Vbar=Vbar,
family=family)
nobs <- nrow(obs)
## (b) declaring the 'TSMLCARA' object
tsml.cara <- TSMLCARA(what=what,
flavor=flavor,
obs=obs,
tm.ref=tm.ref,
tm.model=tm.model,
tm.control=tm.control,
learnQ=learnQ,
Gmin=Gmin,
Gexploit=Gexploit,
Gexplore=Gexplore,
Qmin=Qmin,
conf.level=conf.level,
verbose=verbose)
## update of 'tsml.cara'
update(tsml.cara, verbose=verbose)
targetPsi(tsml.cara, verbose=verbose)
## if 'what=="MOR"', then compute and store the targeted data-adaptive parameters
if (what=="MOR") {
Qn <- getQ(tsml.cara)
ruleQn <- ruleQbar(Qn)
##
## first data-adaptive parameter:
## mean reward under the current best estimate of the optimal rule
##
truthn <- getSample(Bn, tm=oneOne, rule=ruleQn,
piV=piV, Qbar=Qbar, Vbar=Vbar, what,
family=family, slice.by=slice.by)[c("psi", "psi.sd")]
##
## second data-adaptive parameter:
## empirical cumulated regret
##
W <- extractW(obs)
rA <- as.numeric(ruleQn(W))
eregretn <- mean(obs[, "Y"] - Qbar(cbind(A=rA, W)))
##
## third data-adaptive parameter:
## counterfactual cumulated regret
##
col <- match(c("Y0", "Y1"), colnames(obs))
rA <- col[rA+1]
counterfactual <- as.numeric(obs[cbind(1:nrow(obs), rA)])
cregretn <- mean(obs[, "Y"] - counterfactual)
}
## successive updates of 'tsml.cara'
while (nobs < nmax) {
if (nobs+by>nmax) {
by <- (nmax-nobs)
}
newObs <- getSample(by, tm=getGstar(tsml.cara),
piV=piV, Qbar=Qbar, Vbar=Vbar,
family=family)
addNewSample(tsml.cara, newObs)
nobs <- nrow(getObs(tsml.cara))
update(tsml.cara, verbose=verbose)
targetPsi(tsml.cara, verbose=verbose)
if (what=="MOR") {
Qn <- getQ(tsml.cara)
ruleQn <- ruleQbar(Qn)
##
## first data-adaptive parameter:
## mean reward under the current best estimate of the optimal rule
##
truthn <- rbind(truthn,
getSample(Bn, tm=oneOne, rule=ruleQn,
piV=piV, Qbar=Qbar, Vbar=Vbar, what,
family=family, slice.by=slice.by)[c("psi", "psi.sd")])
##
## second data-adaptive parameter:
## empirical cumulated regret
##
obs <- getObs(tsml.cara)
W <- extractW(obs)
rA <- as.numeric(ruleQn(W))
eregretn <- c(eregretn,
mean(obs[, "Y"] - Qbar(cbind(A=rA, W))))
##
## third data-adaptive parameter:
## counterfactual cumulated regret
##
rA <- col[rA+1]
counterfactual <- as.numeric(obs[cbind(1:nrow(obs), rA)])
cregretn <- c(cregretn,
mean(obs[, "Y"] - counterfactual))
}
}
if (what=="MOR") {
truthn[, "psi.sd"] <- truthn[, "psi.sd"]/sqrt(Bn)
colnames(truthn) <- c("psin", "psin.sd")
rownames(truthn) <- NULL
attr(tsml.cara, "truthn") <- truthn
attr(tsml.cara, "eregretn") <- eregretn
attr(tsml.cara, "cregretn") <- cregretn
}
return(tsml.cara)
### Returns a \code{TSMLCARA} object which summarizes the TSMLCARA undertaken
### procedure.
}, ex=function() {
##
log <- Arguments$getVerbose(-1, timestamp=TRUE)
set.seed(12345)
## ########################
## AVERAGE TREATMENT EFFECT
## ########################
tm.model <- formula(A~.)
psi.sd <- sqrt(getOptVar(n=1e5,
tm.model=tm.model,
piV=c(1/2, 1/3, 1/6),
family="gamma",
Qbar=Qbar1,
Vbar=Vbar1))
truth <- c(psi=91/72, psi.sd=psi.sd)
## parametric example
learnQ <- formula(Y~I(as.integer(A)):(U+V)+I(as.integer(1-A)):(U+V))
ATE.param <- tsml.cara.rct(what="ATE",
flavor="parametric",
ninit=200,
by=100,
nmax=400,
tm.init=oneOne,
tm.ref=oneOne,
learnQ=learnQ,
tm.model=tm.model,
conf.level=0.95,
piV=c(1/2, 1/3, 1/6),
family="gamma",
Qbar=Qbar1,
Vbar=Vbar1)
ATE.param
## Not run:
plot(ATE.param, truth=truth)
## End(**Not run**)
## See the vignette for more examples...
})
############################################################################
## HISTORY:
## 2016-09-16
## o Created.
############################################################################
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.