# R/stmodelKM.R In stepp: Subpopulation Treatment Effect Pattern Plot (STEPP)

#### Documented in stepp.KM

```#################################################################
#
# stmodelKM.R
#
#######################
# stepp model: KM     #
#######################
#
kmest1 <- function(y,m,n,ndf,t,s,v,ntpt,tpt,nrr,ndd) {
# initialization
f<-1.0
kr<-n
nrr[1]<-n
ndd[1]<-0

for (i in 2:ntpt){
nrr[i]<-0
ndd[i]<-0
}

ltp<-1
var<-0.0
l<-1
t[1]<-0
s[1]<-1
v[1]<-0
i<-1

# main loop
while (i<= n) {
k<-i+1
k2<-0

while (k2<=0) {
if (k > n) k2<-1
else {
if (y[k] != y[i]) k2<-1 else k<-k+1
}

} # end while

k<-k-1
nd<-0

for (j in i:k) nd<-nd+m[j]

while (ltp<=ntpt && y[i]>tpt[ltp+1]) {
ltp<-ltp+1
nrr[ltp]<-kr
}  # end while

ndd[ltp]<-ndd[ltp]+nd
if (nd>0) {
t1<- nd / kr
f<-f*(1-t1)
if (nd<kr) var<-var+t1/(kr-nd)
t[l+1]<-y[i]
s[l+1]<-s[l]
v[l+1]<-v[l]
l<-l+2
t[l]<-y[i]
s[l]<-f
v[l]<-var*f*f
} # end if

i<-k+1
kr<-n-k     # kr<-n-k-1
k<-i
} # end while - main loop

l<-l+1
t[l]<-y[n]
s[l]<-s[l-1]
v[l]<-v[l-1]

# return a list of all the arguments
kmtest1<-list(y=y,m=m,n=n,ndf=ndf,t=t,s=s,v=v,ntpt=ntpt,tpt=tpt,nrr=as.integer(nrr),ndd=as.integer(ndd))
kmtest1

}

kmest <- function (time,status,group,tpt,pv=TRUE,pv.strat,pv.sub,rho=0,subset,na.action=na.omit) {
d <- data.frame(time=time,status=status,
group=as.factor(if (missing(group)) rep(1,length(time)) else group),
pv.strat=as.factor(if (missing(pv.strat)) rep(1,length(time)) else pv.strat),
pv.sub = as.factor(if (missing(pv.sub)) rep(TRUE, length(time)) else pv.sub))
if (!missing(subset)) d <- d[subset,]
tmp <- nrow(d)
d <- na.action(d)
if (nrow(d) != tmp) cat(format(tmp-nrow(d)),'cases omitted due to missing values\n')
no <- nrow(d)
if (any(d\$status != 0 & d\$status != 1)) stop("invalid status values")
T1 <- table(d\$group)
subgl <- T1>0
T8 <- d\$group[d\$status==1]
if (length(T8) > 0) {
T8 <- cbind(T1,table(d\$group[d\$status==1]))
} else {
T8 <- cbind(T1,rep(0,length(T8)))
}
T8 <- T8[subgl,,drop=FALSE]
T1 <- names(T1)[subgl]
time <- ifelse(d\$time<=0,.00001,d\$time)
if (missing(tpt)) {
tpt <- pretty(time)}
else {
ym <- round(max(time),2)
tpt <- c(0,tpt,ym)
}
ntpt <- length(tpt)-1
lev <- vector("character",ntpt)
for (i in 1:ntpt) lev[i] <- paste(format(tpt[i]),format(tpt[i+1]),sep="-")
nrd <- rep(0, ntpt)
o <- order(time)
Tl <- length(T1)
z <- as.list(1:Tl)
for(i in 1:Tl) {
Ty <- (time[o])[d\$group[o]==T1[i]]
Tm <- (d\$status[o])[d\$group[o]==T1[i]]
ndf <- length(unique(Ty[Tm==1]))
t <- double(2*ndf+2)

a <- kmest1(Ty, Tm, length(Ty), ndf, t,t,t,ntpt, tpt, nrd, nrd)

tt <- paste(a[[11]],a[[10]],sep="/")
names(tt) <- lev
z[[i]] <- list(time=a[[5]],est=a[[6]],var=a[[7]],tint=tt,nnd=T8[i,])
}
names(z) <- T1
if (pv & Tl > 1) {
stop("internal error")
# pv <- logrank(time=d\$time,status=d\$status,group=d\$group,strata=d\$pv.strat,rho=rho,subset=(d\$pv.sub=='TRUE'))\$pv
# attr(z,"pv") <- pv
}

class(z) <- "kmest"
z
}

#
tpest1 <- function(x,n,ind,tp5,ntp) {
l <- ntp
for (i in rev(1:ntp)){
if (x[n] >= tp5[i]) break;
ind[l]<-0
l <- l - 1
}

if (l <=0) return (ind);

if (x[n] == tp5[l]){
ind[l] <- n
l <- l - 1
}

# assuming unique values sorted in ascending order
k <- n-1
loop <- TRUE
while(loop){
if (l <= 0) return (ind);

loop <- FALSE
for (i in rev(1:k)) {
if (x[k] <= tp5[l]){
ind[l] <- k+1
l <- l - 1
loop <- TRUE
break #out of the for loop
} else
k <- k-1
} # end for loop
} # end while loop

# error in the following loop corrected 9-28-04
for (i in 1:l) ind[i] <- 0
return (ind);
}

tpest <- function(w,times) {
if (!is.null(w\$Tests)) w <- w[names(w) != 'Tests']
ng <- length(w)
times <- sort(unique(times))
nt <- length(times)
storage.mode(times) <- "double"
storage.mode(nt) <- "integer"
ind <- matrix(0,ncol=nt,nrow=ng)
oute <- matrix(NA,ncol=nt,nrow=ng)
outv <- oute
storage.mode(ind) <- "integer"
slct <- rep(TRUE,ng)
for (i in 1:ng) {
if (is.null((w[[i]])\$est)) { slct[i] <- FALSE} else {
z1 <- as.integer(tpest1(w[[i]][[i]], length(w[[i]][[i]]), ind[i,], times, nt))
ind[i,] <- z1
oute[i, ind[i,]>0] <- w[[i]][[2]][z1]
if (length(w[[i]])>2) outv[i,ind[i,]>0] <- w[[i]][[3]][z1]
}
}

dimnames(oute) <- list(names(w)[1:ng],as.character(times))
dimnames(outv) <- dimnames(oute)
list(est=oute[slct,,drop=FALSE],var=outv[slct,,drop=FALSE])
}

setClass("stmodelKM",
representation(
coltrt    = "numeric",	# treatment
survTime  = "numeric",  # time to event
censor    = "numeric",	# time to censor
trts      = "numeric",	# trt encoding
timePoint = "numeric"	  # evaluated time
)
)

setMethod("estimate",
signature = "stmodelKM",
definition = function(.Object, sp, ...) {
modtype   <- sp@win@type
nsubpop   <- sp@nsubpop
subpop    <- sp@subpop
coltrt    <- .Object@coltrt
trts	    <- .Object@trts
timePoint <- .Object@timePoint
survTime  <- .Object@survTime
censor    <- .Object@censor

# set up internal assignment
ntrts	    <- length(trts)
txassign  <- rep(NA, length(coltrt))

for (j in 1:ntrts)
txassign[which(coltrt == trts[j])] <- j

# set up return structure
Obs	<- list(sObs = rep(0,nsubpop),
sSE  = rep(0,nsubpop),
oObs = 0,
oSE  = 0)
TrtEff <- array(list(Obs),ntrts)

HRs <- list(skmw     = 0,
logHR    = rep(0,nsubpop),
logHRSE  = rep(0,nsubpop),
ologHR   = 0,
ologHRSE = 0,
logHRw   = 0)
Ratios <- array(list(HRs),ntrts-1)

# Estimate treatment effect for each treatment j
for (j in 1:ntrts) {

skmObs <- rep(0,nsubpop)
skmSE  <- rep(0,nsubpop)
for (i in 1:nsubpop) {
trteff    <- tpest(kmest(survTime[txassign==j & subpop[,i]==1],
censor[txassign==j & subpop[,i]==1]),timePoint)
skmObs[i] <- max(trteff\$est,0)
skmSE[i]  <- sqrt(trteff\$var)
}
overalltrteff <- tpest(kmest(survTime[txassign==j],censor[txassign==j]),timePoint)
overallSkmObs <- max(overalltrteff\$est, 0)
overallSkmSE  <- sqrt(overalltrteff\$var)

# check to make sure that subpopulation estimates are OK.
if (sum(is.na(skmObs)) != 0 | is.na(overallSkmObs) != FALSE) {
cat("\n")
print(paste("Unable to estimate survival time at ", timePoint, " time-unit(s) for trt ",
j, "because there are too few events within one or more subpopulation(s)."))
print(paste("The problem may be avoided by constructing larger subpopulations and/or by selecting a different timepoint for estimation."))
stop()
}

TrtEff[[j]]=list(sObs=skmObs, sSE=skmSE, oObs=overallSkmObs, oSE=overallSkmSE)
}

# Estimate the relative treatment effect comparing each one with trt 1
for (j in 2:ntrts) {
txassign <- rep(-1, length(coltrt))
txassign[which(coltrt == trts[1])] <- 1
txassign[which(coltrt == trts[j])] <- 0

logHR    <- rep(0,nsubpop)
logHRSE  <- rep(0,nsubpop)
sel      <- txassign==1 | txassign==0 #j
for (i in 1:nsubpop) {
seli       <- sel & subpop[,i]==1
LogRank    <- survdiff(Surv(survTime[seli],censor[seli])~txassign[seli])
logHR[i]   <- -(LogRank\$obs[1]-LogRank\$exp[1])/LogRank\$var[1,1]
logHRSE[i] <- sqrt(1/LogRank\$var[1,1])
}
LogRank        <- survdiff(Surv(survTime[sel],censor[sel]) ~ txassign[sel])
overallLogHR   <- -(LogRank\$obs[1]-LogRank\$exp[1])/LogRank\$var[1,1]
overallLogHRSE <- sqrt(1/LogRank\$var[1,1])
logHRw <- sum(logHR/logHRSE)
skmw   <- sum((TrtEff[[1]]\$sObs-TrtEff[[j]]\$sObs)/sqrt(TrtEff[[1]]\$sSE^2 + TrtEff[[j]]\$sSE^2))

Ratios[[j-1]] = list(skmw     = skmw,
logHR    = logHR,
logHRSE  = logHRSE,
ologHR   = overallLogHR,
ologHRSE = overallLogHRSE,
logHRw   = logHRw)

}

# return the estimate as a list
estimate <- list( model	   = "KMe",
ntrts	   = ntrts,
TrtEff   = TrtEff,
Ratios   = Ratios
)
return(estimate)
}
)

setMethod("test",
signature = "stmodelKM",
definition = function(.Object, nperm=2500, sp, effect, showstatus=TRUE, Cox=FALSE, MM=NULL, ...) {

test <- NULL

# no permutation test is done if nperm is 0; return immediately
if (nperm > 0){
win		<- sp@win
nsubpop	<- sp@nsubpop
osubpop	<- sp@subpop
coltrt	<- .Object@coltrt
trts	<- .Object@trts
survTime  <- .Object@survTime
timePoint <- .Object@timePoint
censor    <- .Object@censor

# set up internal assignment
ntrts	<- length(trts)
txassign  <- rep(-1, length(coltrt))	# do not use NA here

for (j in 1:ntrts)
txassign[which(coltrt == trts[j])] <- j

# set up return structure
tsigma    <- matrix(rep(0, (nperm * nsubpop)), ncol = nsubpop)
PermTest	<- list(sigma    = tsigma,
HRsigma    = tsigma,
pvalue     = 0,
chi2pvalue = 0,
HRpvalue   = 0)
Res       <- array(list(PermTest),ntrts-1)

for (j in 2:ntrts){
#
#   do the permutations
#	compare trt1 with the rest
#
differences <- matrix(rep(0, (nperm * nsubpop)), ncol = nsubpop)
logHRs      <- matrix(rep(0, (nperm * nsubpop)), ncol = nsubpop)
no          <- 0
p           <- 0
terminate   <- 0
Ntemp       <- nrow(osubpop)
IndexSet1   <- (1:Ntemp)[txassign == 1]
IndexSet2   <- (1:Ntemp)[txassign == j]

if (Cox) {
if (is.null(MM)){
fmla <- as.formula(paste("Surv(survTime,censor) ~ txassign"))
} else {
xnam <- colnames(MM)[2:dim(MM)[2]]
fmla <- as.formula(paste("Surv(survTime,censor) ~ txassign + ", paste(xnam, collapse= "+")))
}
}

if (showstatus){
title <- paste("\nComputing the p-value with", nperm)
title <- paste(title, "number of permutations comparing trt", trts[j], "with trt", trts[1], "\n")
cat(title)
pb    <- txtProgressBar(min=0, max=nperm-1, style=3)
}

while (no < nperm) {
## PERMUTE THE VECTOR OF SUBPOPULATIONS WITHIN TREATMENTS ##
if (showstatus) setTxtProgressBar(pb, no)

Subpop        <- as.matrix(osubpop)
permuteSubpop <- matrix(0, nrow = Ntemp, ncol = nsubpop)
permuteSubpop[txassign == 1] <- Subpop[sample(IndexSet1),]
permuteSubpop[txassign == j] <- Subpop[sample(IndexSet2),]
subpop        <- permuteSubpop
skm1          <- rep(0, nsubpop)
skm2          <- rep(0, nsubpop)
skmSE1	  <- rep(0, nsubpop)
skmSE2	  <- rep(0, nsubpop)
slogHRs       <- rep(0, nsubpop)
slogHRSEs     <- rep(0, nsubpop)

for (i in 1:nsubpop) {
seff1	  <- tpest(kmest(survTime[txassign == 1 & subpop[, i] == 1], censor[txassign == 1 & subpop[, i] == 1]), timePoint)
seff2	  <- tpest(kmest(survTime[txassign == j & subpop[, i] == 1], censor[txassign == j & subpop[, i] == 1]), timePoint)
skm1[i]     <- max(seff1\$est, 0)
skm2[i]     <- max(seff2\$est, 0)
skmSE1[i]   <- sqrt(seff1\$var)
skmSE2[i]   <- sqrt(seff2\$var)
sel         <- (txassign == 1 | txassign == j) & subpop[,i]==1
if (Cox) {
CoxM       	<- coxph(fmla, subset=(sel))
if (is.null(MM)){
slogHRs[i]    <- CoxM\$coefficient
slogHRSEs[i]  <- sqrt(CoxM\$var)
} else {
slogHRs[i]    <- CoxM\$coefficient[1]
slogHRSEs[i]  <- sqrt(CoxM\$var[1,1])
}
} else {
LogRank         <- survdiff(Surv(survTime[sel],censor[sel]) ~ txassign[sel])
slogHRs[i]      <- -(LogRank\$obs[1]-LogRank\$exp[1])/LogRank\$var[1,1]
slogHRSEs[i]    <- sqrt(1/LogRank\$var[1,1])
}
}
selo         	  <- (txassign == 1 | txassign == j)

slogHRw		  <- sum(slogHRs/slogHRSEs)

overallSkm1         <- max(tpest(kmest(survTime[txassign == 1], censor[txassign == 1]), timePoint)\$est, 0)
overallSkm2         <- max(tpest(kmest(survTime[txassign == j], censor[txassign == j]), timePoint)\$est, 0)

if (Cox) {
CoxM		       <- coxph(fmla, subset=(selo))
if (is.null(MM)){
overallSLogHR    <- CoxM\$coefficient
overallLogHRSE   <- sqrt(CoxM\$var)
} else {
overallSLogHR    <- CoxM\$coefficient[1]
overallSLogHRSE  <- sqrt(CoxM\$var[1,1])
}
} else {
LogRank            <- survdiff(Surv(survTime[selo],censor[selo]) ~ txassign[selo])
overallSLogHR      <- -(LogRank\$obs[1]-LogRank\$exp[1])/LogRank\$var[1,1]
overallSLogHRSE    <- sqrt(1/LogRank\$var[1,1])
}
if (sum(is.na(skm1)) == 0 & sum(is.na(skm2)) == 0 & is.na(overallSkm1) == FALSE & is.na(overallSkm2) == FALSE) {
no <- no + 1
p <- p + 1
for (s in 1:nsubpop) {
differences[p, s] <- (skm1[s] - skm2[s]) - (overallSkm1 - overallSkm2)
logHRs[p,s]       <- slogHRs[s] - overallSLogHR
}
}
terminate <- terminate + 1
if (terminate >= nperm + 10000) {
print(paste("After permuting ", nperm, "plus 10000, or ",
nperm + 10000, " times, the program is unable to generate the permutation distribution based on ",
nperm, "permutations of the data"))
print(paste("Consider creating larger subpopulations or selecting a different timepoint for estimation"))
print(paste("when comparing trt ",trts[1], "vs trt ",trts[j]))
stop()
}
}
# generating the sigmas and p-values
sObs   <- effect\$TrtEff[[1]]\$sObs-effect\$TrtEff[[j]]\$sObs
oObs   <- effect\$TrtEff[[1]]\$oObs-effect\$TrtEff[[j]]\$oObs
logHR  <- effect\$Ratios[[j-1]]\$logHR
ologHR <- effect\$Ratios[[j-1]]\$ologHR

mname  <- paste0("SP",seq(1,dim(differences)[2]),"-Overall")
if (win@type == "tail-oriented"){
# remove the trivial case of the full cohort
rm <- length(win@r1)+1
differences <- differences[,-rm]
mname <- mname[-rm]
sObs <- sObs[-rm]
oObs <- oObs[-rm]
logHRs <- logHRs[,-rm]
logHR <- logHR[-rm]
ologHR <- ologHR[-rm]
}

sigmas        		<- ssigma(differences)
rownames(sigmas\$sigma)	<- mname
colnames(sigmas\$sigma)	<- mname
Res[[j-1]]\$sigma		<- sigmas\$sigma
Res[[j-1]]\$chi2pvalue     <- ppv (differences, sigmas\$sigmainv, sObs, oObs, nperm)
Res[[j-1]]\$pvalue     	<- ppv2(differences, sObs, oObs, nperm)
sigmas        		<- ssigma(logHRs)
rownames(sigmas\$sigma)	<- mname
colnames(sigmas\$sigma)	<- mname
Res[[j-1]]\$HRsigma    	<- sigmas\$sigma
Res[[j-1]]\$HRpvalue   	<- ppv2(logHRs, logHR, ologHR, nperm)

if (showstatus) close(pb)
}
test = list( model  = "KMt",
ntrts  = ntrts,
Res    = Res
)
}
return(test)
}
)

#
# printing support functions for KM model
#
print.estimate.KM <- function(x, timePoint, trts) {
for (j in 1:x@effect\$ntrts){
cat("\n")
write(paste("Survival estimates for treatment group", trts[j],
"at time point", timePoint), file = "")

temp <- matrix(c(1:x@subpop@nsubpop, round(x@effect\$TrtEff[[j]]\$sObs, digits = 4),
round(x@effect\$TrtEff[[j]]\$sSE, digits = 4)), ncol = 3)
write("                         Survival", file = "")
write("     Subpopulation     Probability      Std. Err.",
file = "")
for (i in 1:x@subpop@nsubpop) {
if (x@subpop@win@type == "tail-oriented" & i == (length(x@subpop@win@r1)+1)){
write(paste(format(temp[i, 1], width = 12), format(temp[i,2], width = 19, nsmall = 4),
format(temp[i, 3], width = 15, nsmall = 4), "Overall"), file = "")
} else {
write(paste(format(temp[i, 1], width = 12), format(temp[i,2], width = 19, nsmall = 4),
format(temp[i, 3], width = 15, nsmall = 4)), file = "")
}
}
if (x@subpop@win@type != "tail-oriented"){
write(paste("        Overall",
format(round(x@effect\$TrtEff[[j]]\$oObs, digits = 4), nsmall = 4, width = 16),
format(round(x@effect\$TrtEff[[j]]\$oSE, digits = 4), nsmall = 4, width = 15)),
file = "")
}
cat("\n")
}

cat("\n")
write("Survival differences at time point and hazard ratio estimates", file = "")

for (j in 2:x@effect\$ntrts){
cat("\n")
write(paste("trt ", trts[j], "vs. trt ", trts[1]), file = "")

temp <- matrix(c(1:x@subpop@nsubpop,
round(x@effect\$TrtEff[[1]]\$sObs - x@effect\$TrtEff[[j]]\$sObs, digits = 4),
round(sqrt(x@effect\$TrtEff[[1]]\$sSE^2 + x@effect\$TrtEff[[j]]\$sSE^2), digits = 4)), ncol = 3)
cat("\n")
write(paste("Survival differences at time point", timePoint), file = "")

write(paste("Comparing trt ", trts[j], " vs. trt ", trts[1]), file = "")
cat("\n")
write("                         Survival", file = "")
write("     Subpopulation      Difference      Std. Err.",
file = "")
for (i in 1:x@subpop@nsubpop) {
if (x@subpop@win@type == "tail-oriented" & i == (length(x@subpop@win@r1)+1)){
write(paste(format(temp[i, 1], width = 12), format(temp[i, 2], width = 19, nsmall = 4),
format(temp[i, 3], width = 15, nsmall = 4), "Overall"), file = "")
} else {
write(paste(format(temp[i, 1], width = 12), format(temp[i, 2], width = 19, nsmall = 4),
format(temp[i, 3], width = 15, nsmall = 4)), file = "")
}
}
if (x@subpop@win@type != "tail-oriented"){
write(paste("        Overall",
format(round(x@effect\$TrtEff[[1]]\$oObs - x@effect\$TrtEff[[j]]\$oObs, digits = 4), nsmall = 4, width = 16),
format(round(sqrt(x@effect\$TrtEff[[1]]\$oSE^2 + x@effect\$TrtEff[[j]]\$oSE^2), digits = 4), nsmall = 4, width = 15)),
file = "")
}
cat("\n")
write("Hazard ratio estimates", file = "")
temp <- matrix(c(1:x@subpop@nsubpop,
round(x@effect\$Ratios[[j-1]]\$logHR, digits = 6),
round(x@effect\$Ratios[[j-1]]\$logHRSE, digits = 6),
round(exp(x@effect\$Ratios[[j-1]]\$logHR), digits = 2)), ncol = 4)
write("     Subpopulation        Log HR       Std. Err.       Hazard Ratio",
file = "")
for (i in 1:x@subpop@nsubpop) {
if (x@subpop@win@type == "tail-oriented" & i == (length(x@subpop@win@r1)+1)){
write(paste(format(temp[i, 1], width = 12),
format(temp[i, 2], width = 19, nsmall = 6),
format(temp[i, 3], width = 14, nsmall = 6),
format(temp[i, 4], width = 15, nsmall = 2), "Overall"),
file = "")
}
else {
write(paste(format(temp[i, 1], width = 12),
format(temp[i, 2], width = 19, nsmall = 6),
format(temp[i, 3], width = 14, nsmall = 6),
format(temp[i, 4], width = 15, nsmall = 2)),
file = "")
}
}
if (x@subpop@win@type != "tail-oriented"){
write(paste("        Overall",
format(round(x@effect\$Ratios[[j-1]]\$ologHR,      digits = 6), nsmall = 6, width = 16),
format(round(x@effect\$Ratios[[j-1]]\$ologHRSE,    digits = 6), nsmall = 6, width = 14),
format(round(exp(x@effect\$Ratios[[j-1]]\$ologHR), digits = 2), nsmall = 2, width = 15)),
file = "")
}
cat("\n")
}
}

print.cov.KM <- function(stobj, timePoint, trts) {
if (!is.null(stobj@result)) {
for (j in 1:(stobj@result\$ntrts-1)){
ns <- stobj@subpop@nsubpop
if (stobj@subpop@win@type == "tail-oriented") ns <- ns-1
cat("\n")
write(paste("The covariance matrix of the Kaplan-Meier differences at",
timePoint, "time units for the", ns, "subpopulations is:"),
file = "")
write(paste("trt ", trts[j+1], "vs. trt ", trts[1]), file = "")
print(stobj@result\$Res[[j]]\$sigma)

cat("\n")
write(paste("The covariance matrix of the log hazard ratios for the",
ns, "subpopulations is:"), file = "")
print(stobj@result\$Res[[j]]\$HRsigma)
cat("\n")

#write(paste("The covariance matrix (based on homogeneous association) of the Kaplan-Meier differences at",
#     	timePoint, "time units for the", stobj@subpop@nsubpop, "subpopulations is:"),
#     	file = "")
#print(stobj@result\$hasigma)

#cat("\n")
#write(paste("The covariance matrix (based on homogeneous association) of the log hazard ratios for the",
#    	stobj@subpop@nsubpop, "subpopulations is:"), file = "")
#print(stobj@result\$haHRsigma)
#cat("\n")
}
}
}

print.stat.KM <- function(stobj, trts) {
if (!is.null(stobj@result)) {
for (j in 1:(stobj@result\$ntrts-1)){

t <- stobj@result\$Res[[j]]
cat("\n")
write(paste("Supremum test results"), file = "")
write(paste("trt", trts[j + 1], "vs. trt", trts[1]), file = "")
write(paste("Interaction p-value based on Kaplan-Meier estimates:", t\$pvalue), file = "")
write(paste("Interaction p-value based on hazard ratio estimates:", t\$HRpvalue), file = "")

cat("\n")
write(paste("Chi-square test results"), file = "")
write(paste("Interaction p-value based on Kaplan-Meier estimates:",
t\$chi2pvalue), file = "")

#cat("\n")
#write(paste("Homogeneous association test results"), file = "")
#write(paste("Interaction p-value based on Kaplan-Meier estimates:",
#	hapvalue), file = "")
#write(paste("Interaction p-value based on hazard ratio estimates:",
#	haHRpvalue), file = "")

cat("\n")
}
}
}

setMethod("print",
signature=  "stmodelKM",
definition=  function(x, stobj, estimate=TRUE, cov=TRUE, test=TRUE, ...) {

#
#  1. estimates
#
if (estimate){
print.estimate.KM(stobj, x@timePoint, x@trts)
}

#
#   2. covariance matrices
#
if (cov){
print.cov.KM(stobj, x@timePoint, x@trts)
}

#
#   3. Supremum test and Chi-square test results
#
if (test){
print.stat.KM(stobj, x@trts)
}
}
)

# constructor function for stmodelKM
stepp.KM <- function(coltrt, survTime, censor, trts, timePoint){
model <- new("stmodelKM", coltrt=coltrt, survTime=survTime, censor=censor,
trts=trts, timePoint=timePoint)
return(model)
}
```

## Try the stepp package in your browser

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

stepp documentation built on Jan. 13, 2021, 5:25 p.m.