# R/WRSEestimate.R In DTR: Estimation and Comparison of Dynamic Treatment Regimes

#### Documented in sub.WRSEestimateWRSEestimate

```###################################################
### Reference:
### Guo X, Tsiatis AA: A weighted risk set estimator for survival distributions in two-stage
### randomization designs with censored survival data. Int. J. Biostatistics 1:1-15, 2005
###################################################

###################################################
### code chunk number 1: chunklibraries
###################################################

require(survival)

###################################################
### code chunk number 2: chunkWRSE
###################################################

#Subfunction to calculate among each arm (A1 or A2)
sub.WRSEestimate <- function(pdata, t) {

#Retrieve data
n <- nrow(pdata)
TR <- pdata\$TR
R <- pdata\$R
Z <- pdata\$Z # Z=0 for B1, 1 for B2
U <- pdata\$U
delta <- pdata\$delta

#Calculate probability of Z given R
pi.z <- sum(R*Z) / sum(R)

#Calculate time-varying weight function for each observed Ui
#W1 for B1; W2 for B2
#Time as rows, individual as columns
w1 <- w2 <- matrix(0, nrow=n, ncol=n)
#Calculate summation of Wj(Ui)Yj(Ui)
s1 <- s2 <- rep(0, n)
#Calculate Yi(Uk) for second part of variance sigma
#Each individual time as one row
kind <- matrix(0, nrow=n, ncol=n)

for(i in 1:n) {

#Calculate I(TR <= Ui)
rind <- as.numeric(TR <= U[i])

#Calculate individual weights for each time Ui
w1[i,] <- 1 - R*rind + R*rind*(1-Z)/(1-pi.z) # weighting for A1B1
w2[i,] <- 1 - R*rind + R*rind*Z/pi.z # weighting for A1B2

#Calculate Yj(ui)
yind <- as.numeric(U >= U[i])

#Caculate summation Wj(Ui)Yj(Ui) for each time Ui
s1[i] <- sum(w1[i,]*yind)
s2[i] <- sum(w2[i,]*yind)

#Calculate Yi(Uk) for second part of variance sigma
kind[i,] <- as.numeric(U[i] >= U)

}

#Calculate the estimate and standard error for each time point
SURV1 <- SURV2 <- rep(NA, length(t))
SE1 <- SE2 <- COV12 <- rep(NA, length(t))

for(j in 1:length(t)) {

ind <- as.numeric(U <= t[j])

#Calculate the survival estimates
SURV1[j] <- exp(-sum(diag(w1)[which(s1!=0)]*delta[which(s1!=0)]*ind[which(s1!=0)]/s1[which(s1!=0)]))
SURV2[j] <- exp(-sum(diag(w2)[which(s2!=0)]*delta[which(s2!=0)]*ind[which(s2!=0)]/s2[which(s2!=0)]))

#Calculate first part (pI) of the variance sigma
#Calculate second part (pII) of the variance sigma
pI1 <- pI2 <- rep(0, n)
pII1 <- pII2 <- rep(0, n)
sdp1 <- sdp2 <- sdp12 <- rep(0, n)

for(k in 1:n) {

#Calculate first part for each individual
if(s1[k] != 0) pI1[k] <- w1[k,k]*delta[k]*ind[k]/s1[k]
if(s2[k] != 0) pI2[k] <- w2[k,k]*delta[k]*ind[k]/s2[k]

#Calculate second part for each individual
pII1[k] <- sum(w1[which(s1!=0),k]*kind[k,][which(s1!=0)]*diag(w1)[which(s1!=0)]*delta[which(s1!=0)]*ind[which(s1!=0)]/s1[which(s1!=0)]^2)
pII2[k] <- sum(w2[which(s2!=0),k]*kind[k,][which(s2!=0)]*diag(w2)[which(s2!=0)]*delta[which(s2!=0)]*ind[which(s2!=0)]/s2[which(s2!=0)]^2)

#Calcualte the square of the difference between first part and second part
sdp1[k] <- (pI1[k] - pII1[k])^2
sdp2[k] <- (pI2[k] - pII2[k])^2
sdp12[k] <- (pI1[k] - pII1[k])*(pI2[k] - pII2[k])

}

#Calculate the variance and covariace sigma
sig1 <- n*sum(sdp1)
sig2 <- n*sum(sdp2)
sig12 <- n*sum(sdp12)

#Calculate the variance and covariace
v1 <- SURV1[j]^2*sig1/n; SE1[j] <- sqrt(v1)
v2 <- SURV2[j]^2*sig2/n; SE2[j] <- sqrt(v2)
COV12[j] <- SURV1[j]*SURV2[j]*sig12/n

}

#Return results
est <- data.frame(t, SURV1, SURV2, SE1, SE2, COV12)
return(est)

}

WRSEestimate <- function(data # A data frame representing the data from sequentially randomized designs
# data = data frame {X, TR, R, Z, U, delta}
) {

#Chek for errors
if (is.null(data\$X)) stop("X can not be empty")
if (is.null(data\$TR)) stop("TR can not be empty")
if (is.null(data\$R)) stop("R can not be empty")
if (is.null(data\$Z)) stop("Z can not be empty")
if (is.null(data\$U)) stop("U can not be empty")
if (is.null(data\$delta)) stop("delta can not be empty")

#Times to be assessed
t <- unique(data\$U[which(data\$delta==1)])
#Order times
t <- t[order(t)]

#Number at risk
n.risk <- apply(as.array(t), 1, function(x) sum(as.numeric(data\$U >= x)))

#Number event
n.event <- apply(as.array(t), 1, function(x) length(which(data\$U==x & data\$delta==1)))

#Run sub.LDTestimate to obtain the defined estimates A1 arm
cat ("Estimating for A1 arm... \n")
est1 <- sub.WRSEestimate(pdata=data[which(data\$X==0),],t)

#Run sub.LDTestimate to obtain the defined estimates A2 arm
cat ("Estimating for A2 arm... \n")
est2 <- sub.WRSEestimate(pdata=data[which(data\$X==1),],t)

#Return class
censorDTR <- censortime <- NULL
results <- list(Call=match.call(),
DTR=c("A1B1", "A1B2", "A2B1", "A2B2"),
records=c(length(which((data\$X==0 & data\$R==0) | (data\$X==0 & data\$R==1 & data\$Z==0))),
length(which((data\$X==0 & data\$R==0) | (data\$X==0 & data\$R==1 & data\$Z==1))),
length(which((data\$X==1 & data\$R==0) | (data\$X==1 & data\$R==1 & data\$Z==0))),
length(which((data\$X==1 & data\$R==0) | (data\$X==1 & data\$R==1 & data\$Z==1)))),
events=c(sum(data\$delta[which((data\$X==0 & data\$R==0) | (data\$X==0 & data\$R==1 & data\$Z==0))]),
sum(data\$delta[which((data\$X==0 & data\$R==0) | (data\$X==0 & data\$R==1 & data\$Z==1))]),
sum(data\$delta[which((data\$X==1 & data\$R==0) | (data\$X==1 & data\$R==1 & data\$Z==0))]),
sum(data\$delta[which((data\$X==1 & data\$R==0) | (data\$X==1 & data\$R==1 & data\$Z==1))])),
censorDTR=c(rep("A1B1", length(which((data\$X==0 & data\$R==1 & data\$Z==0 & data\$delta==0) | (data\$X==0 & data\$R==0 & data\$delta==0)))),
rep("A1B2", length(which((data\$X==0 & data\$R==1 & data\$Z==1 & data\$delta==0) | (data\$X==0 & data\$R==0 & data\$delta==0)))),
rep("A2B1", length(which((data\$X==1 & data\$R==1 & data\$Z==0 & data\$delta==0) | (data\$X==1 & data\$R==0 & data\$delta==0)))),
rep("A2B2", length(which((data\$X==1 & data\$R==1 & data\$Z==1 & data\$delta==0) | (data\$X==1 & data\$R==0 & data\$delta==0))))),
censortime=c(data\$U[which((data\$X==0 & data\$R==1 & data\$Z==0 & data\$delta==0) | (data\$X==0 & data\$R==0 & data\$delta==0))],
data\$U[which((data\$X==0 & data\$R==1 & data\$Z==1 & data\$delta==0) | (data\$X==0 & data\$R==0 & data\$delta==0))],
data\$U[which((data\$X==1 & data\$R==1 & data\$Z==0 & data\$delta==0) | (data\$X==1 & data\$R==0 & data\$delta==0))],
data\$U[which((data\$X==1 & data\$R==1 & data\$Z==1 & data\$delta==0) | (data\$X==1 & data\$R==0 & data\$delta==0))]),
censorsurv=c(apply(as.array(data\$U[which((data\$X==0 & data\$R==1 & data\$Z==0 & data\$delta==0) | (data\$X==0 & data\$R==0 & data\$delta==0))]), 1,
function(x) { if(x<min(est1\$t)) 1 else est1\$SURV1[which(abs(est1\$t-x)==min(abs(est1\$t-x)[which(est1\$t<=x)]))] }),
apply(as.array(data\$U[which((data\$X==0 & data\$R==1 & data\$Z==1 & data\$delta==0) | (data\$X==0 & data\$R==0 & data\$delta==0))]), 1,
function(x) { if(x<min(est1\$t)) 1 else est1\$SURV2[which(abs(est1\$t-x)==min(abs(est1\$t-x)[which(est1\$t<=x)]))] }),
apply(as.array(data\$U[which((data\$X==1 & data\$R==1 & data\$Z==0 & data\$delta==0) | (data\$X==1 & data\$R==0 & data\$delta==0))]), 1,
function(x) { if(x<min(est2\$t)) 1 else est2\$SURV1[which(abs(est2\$t-x)==min(abs(est2\$t-x)[which(est2\$t<=x)]))] }),
apply(as.array(data\$U[which((data\$X==1 & data\$R==1 & data\$Z==1 & data\$delta==0) | (data\$X==1 & data\$R==0 & data\$delta==0))]), 1,
function(x) { if(x<min(est2\$t)) 1 else est2\$SURV2[which(abs(est2\$t-x)==min(abs(est2\$t-x)[which(est2\$t<=x)]))] })),
time=t, n.risk=n.risk, n.event=n.event,
SURV11=est1\$SURV1, SURV12=est1\$SURV2, SURV21=est2\$SURV1, SURV22=est2\$SURV2,
SE11=est1\$SE1, SE12=est1\$SE2, COV1112=est1\$COV12,
SE21=est2\$SE1, SE22=est2\$SE2, COV2122=est2\$COV12)
class(results) <- "DTR"
return(results)

}
```

## Try the DTR package in your browser

Any scripts or data that you put into this service are public.

DTR documentation built on May 30, 2017, 6:38 a.m.