R/delta2samp.R

In bpcp: Beta Product Confidence Procedure for Right Censored Data

# the following functions are used to calculate two-sample confidence intervals based on delta method modifications


parmtype_to_trans <- function(parmtype){
  # get one-sample trans from parmtype
  trans <- switch(parmtype,
                  difference="none",
                  oddsratio="logodds",
                  ratio="log",
                  efflogs="loglog",
                  effcdf="clog",
                cdfratio="clog",
                logsratio="loglog")
}

# the following function runs argument checks and calculates lower/upper limit and nullparm
# adapted from bpcp2samp code, so this could be used there too for less coding redundancy
twosampleChecks <- function(time, testtime, status, group, alternative, parmtype, method, nullparm, conf.level,changeGroupOrder, method_type=c("delta","bpcp") ){
  
  # in case group is a factor with more than 2 levels, but only 2 levels are selected
  # keep only the two levels with data
  #if (class(group)=="factor"){
  if (is(group,"factor")){
    # use order of levels, but only keep the levels with data 
    ug<- levels(group)
    ug<- ug[ug %in% group]
  } else ug<-sort(unique(group))
  
  ## argument checking
  
  alt<-match.arg(alternative, c("two.sided","less","greater"))

  if(method_type=="bpcp") ptype_opt <-  c("difference","oddsratio","ratio","efflogs","effcdf","cdfratio","logsratio","one.minus.ratio","one.minus.cdfratio") else ptype_opt <-  c("difference","oddsratio","ratio","efflogs","effcdf","one.minus.cdfratio")
  ptype<-match.arg(arg=parmtype, choices=ptype_opt, several.ok = T)[1]
  
  if(method_type=="delta") method_opt <-  c("standard","reg_hybrid","adj_hybrid","sh_adj_hybrid") else method_opt <-  c("melded","midp")
  meth<-match.arg(arg=method, choices=method_opt, several.ok = T)[1]
  
  
  if (length(ug)!=2) stop("group does not have 2 levels") 
  if ((length(group[group==ug[1]])<1) | (length(group[group==ug[2]])<1)){ 
    stop("should have at least one observation in each group") }
  
  if (changeGroupOrder){
    ug<-ug[2:1]
  }
  
  
  if (length(time)!=length(status) | length(time)!=length(group) ) stop("length of time, status and group should be the same")
  if(!is.null(testtime)){
    if (length(testtime)>1 | !is.numeric(testtime[1]) | testtime[1]<=0) stop("testtime must be positive and a vector of length 1")
    
    
  }
  
  # change "one.minus.cdfratio" to effcdf
  ptype <- ifelse(ptype=="one.minus.cdfratio","effcdf",ptype)
  
  trans <- parmtype_to_trans(ptype)
  h <- get_h(trans=trans)
  g <- get_g(ptype)
  
  lowerLimit<- g(h(0)-h(1))
  upperLimit<- g(h(1)-h(0))
  
  
  ptype.calc <- switch(ptype,
                       difference="difference",
                       oddsratio="oddsratio",
                       ratio="ratio",
                       efflogs="logsratio",
                       effcdf="cdfratio",
                       cdfratio="cdfratio",
                       logsratio="logsratio",
                       one.minus.ratio="ratio",
                       one.minus.cdfratio="cdfratio")
  
  checkNullparm<-function(nullparm){
    if (nullparm<lowerLimit | nullparm>upperLimit) stop("nullparm out of range of possible values")
  }
  
  if (!is.null(nullparm)) {
    checkNullparm(nullparm)
  Nullparm <- nullparm
  } else{
    g <- get_g(ptype)
  Nullparm <- g(0)
    
  }
  
  ## end of argument checking
  
 
  
  out <- list(ptype=ptype,ptype.calc=ptype.calc,lowerLimit=lowerLimit,upperLimit=upperLimit,
              Nullparm=Nullparm, alt=alt, ug=ug, alpha=1-conf.level, CL=conf.level, trans=trans, meth=meth)
  
  
}

create.htest <- function(testtime,alt,conf.level, lower, upper, p.value, est1, est2, beta,nullparm, ptype, ug,method,zero.one.adjustment=FALSE){
  
  
  ci<-c(lower,upper)
  attr(ci,"conf.level")<-conf.level
    
    full.method <- switch(method,
                          standard="Standard delta method",
                          reg_hybrid="Delta method using Borkowf's regular hybrid variance",
                          adj_hybrid="Delta method using Borkowf's adjusted hybrid variance",
                          sh_adj_hybrid="Delta method using Borkowf's adjusted hybrid variance and shrunken KM estimates",
                          melded="Two-Sample Melded BPCP Test",
                          midp="Two-Sample Melded BPCP Test (mid-p version)"
                          )

  if(zero.one.adjustment){
    full.method=paste0(full.method,", with zero-one adjustment")
  }
  
  stat<-est1
  parm<-est2
  
  # statistic= either S1(t) or F1(t)
  # parameter= either S2(t) or F2(t)
  if (ptype %in% c("cdfratio","one.minus.cdfratio","effcdf")){
    names(stat)<-paste("F(",testtime,";group=",ug[1],")",sep="")
    names(parm)<-paste("F(",testtime,";group=",ug[2],")",sep="")
  }  else {
    names(stat)<-paste("S(",testtime,";group=",ug[1],")",sep="")
    names(parm)<-paste("S(",testtime,";group=",ug[2],")",sep="")
  } 
  names(nullparm)<-paste(ptype)
  
  # get data.name based on parmtype
  data.name <- switch(ptype,
                      difference=paste("S(",testtime,";group=",ug[2],") - S(",testtime,";group=",ug[1],")",sep=""),
                      ratio=paste("S(",testtime,";group=",ug[2],")/ S(",testtime,";group=",ug[1],")",sep=""),
                      oddsratio=paste("odds[S(",testtime,";group=",ug[2],")] / odds[S(",testtime,";group=",ug[1],")]",sep=""),
                      cdfratio=paste("F(",testtime,";group=",ug[1],") / F(",testtime,";group=",ug[2],")",sep=""),
                      logsratio=paste("log[S(",testtime,";group=",ug[1],")] / log[S(",testtime,";group=",ug[2],")]",sep=""),
                      efflogs=paste("1-log[S(",testtime,";group=",ug[2],")] / log[S(",testtime,";group=",ug[1],")]",sep=""),
                      effcdf=paste("1-F(",testtime,";group=",ug[2],") / F(",testtime,";group=",ug[1],")",sep=""),
                      one.minus.ratio=paste("1-S(",testtime,";group=",ug[2],") / S(",testtime,";group=",ug[1],")",sep=""),
                      one.minus.cdfratio=testout$data.name<-paste("1-F(",testtime,";group=",ug[1],") / F(",testtime,";group=",ug[2],")",sep="")
                      )
  
  
  # create htest object for output
  structure(list(statistic = stat, parameter = parm, 
                 p.value = p.value, 
                 conf.int = ci, estimate = beta, 
                 null.value = nullparm, 
                 alternative = alt, method = full.method, 
                 data.name = data.name), class = "htest")
}

create.twosamp <- function(alt,conf.level, L, Lin, R, Rin, interval, g1, est_g1, lower_g1, upper_g1, g2,est_g2, lower_g2, upper_g2, beta, lower, upper,p.value,nullparm, ptype, method,zero.one.adjustment=FALSE){
  

    full.method <- switch(method,
                                standard="Standard delta method",
                                reg_hybrid="Delta method using Borkowf's regular hybrid variance",
                                adj_hybrid="Delta method using Borkowf's adjusted hybrid variance",
                                sh_adj_hybrid="Delta method using Borkowf's adjusted hybrid variance and shrunken KM estimates",
                          midp.m="Monte Carlo with mid-p",
                          melded="melded test using two beta CDs"
    ) 
  
    
    if(zero.one.adjustment){
      
      full.method=paste0(full.method,", with zero-one adjustment")
    }
  
  names(nullparm)<-paste(ptype)
  
  # create object for output
  structure(list(L=L, Lin=Lin, R=R, Rin=Rin, 
                interval=interval,
                 g1=g1,est_group1 = est_g1, lower_group1=lower_g1, upper_group1=upper_g1,
                 g2=g2,est_group2 = est_g2, lower_group2=lower_g2, upper_group2=upper_g2,
                 beta = beta, 
                 lower = lower,
                 upper = upper,
                 conf.level = conf.level, 
                 p.value = p.value, 
                 null.value = nullparm, 
                 alternative = alt, 
                 method = full.method), class = c("twosamp","list"))
}


# summary function for twosamp class
summary.twosamp <- function(object,...){
  x<-object
  out<-data.frame(x$interval,x$est_group1,x$est_group2,x$beta,x$lower,x$upper, x$p.value)
  dimnames(out)[[2]]<-c("time interval","group 1 surv", "group 2 surv", names(x$null.value),
                        paste("lower ",100*x$conf.level,"% CL",sep=""),
                        paste("upper ",100*x$conf.level,"% CL",sep=""),
                        "p.value")
  out
  
}


plot.twosamp <- function(x, ...) {
  # 1. Prepare the data.frame 
  x_df <- data.frame(
    time = c(x$L[which(x$Lin)], x$L[which(!x$Lin & !x$Rin)] + 0.01, x$R[which(x$Rin)]),
    group1_est = c(x$est_group1[which(x$Lin)], x$est_group1[which(!x$Lin & !x$Rin)], x$est_group1[which(x$Rin)]),
    group1_lower = c(x$lower_group1[which(x$Lin)], x$lower_group1[which(!x$Lin & !x$Rin)], x$lower_group1[which(x$Rin)]),
    group1_upper = c(x$upper_group1[which(x$Lin)], x$upper_group1[which(!x$Lin & !x$Rin)], x$upper_group1[which(x$Rin)]),
    group2_est = c(x$est_group2[which(x$Lin)], x$est_group2[which(!x$Lin & !x$Rin)], x$est_group2[which(x$Rin)]),
    group2_lower = c(x$lower_group2[which(x$Lin)], x$lower_group2[which(!x$Lin & !x$Rin)], x$lower_group2[which(x$Rin)]),
    group2_upper = c(x$upper_group2[which(x$Lin)], x$upper_group2[which(!x$Lin & !x$Rin)], x$upper_group2[which(x$Rin)])
  )
  
  # sort by time
  x_df <- x_df[order(x_df$time), ]
  
  # 2. Set up base R plotting: two panels, one for each group
  oldpar <- par(no.readonly=TRUE)
  on.exit(par(oldpar))
  par(mfrow=c(2,1), mar=c(4,4,2,1))
  
  # 3. Plot for group 1
  matplot(
    x_df$time, 
    cbind(x_df$group1_est, x_df$group1_lower, x_df$group1_upper), 
    type="s", lty=c(1,2,2), col=1, lwd=2,
    xlab="Time", ylab="Survival", 
    main=if(!is.null(x$g1)) x$g1 else "Group 1",
    ylim=range(c(x_df$group1_lower, x_df$group1_upper, x_df$group1_est), na.rm=TRUE)
  )
  
  # 4. Plot for group 2
  matplot(
    x_df$time, 
    cbind(x_df$group2_est, x_df$group2_lower, x_df$group2_upper), 
    type="s", lty=c(1,2,2), col=1, lwd=2,
    xlab="Time", ylab="Survival", 
    main=if(!is.null(x$g2)) x$g2 else "Group 2",
    ylim=range(c(x_df$group2_lower, x_df$group2_upper, x_df$group2_est), na.rm=TRUE)
  )
  }

# plot function (tidy version) for twosamp class
plot.twosamptidy <- function(x, ...){

  
  # convert output into data.frame
  x_df <- data.frame(
    time=rep(c(x$L[which(x$Lin)], x$L[which(!x$Lin & !x$Rin)] + 0.01 , x$R[which(x$Rin)]),2),
    group=rep(c("group1","group2"), each=length(x$L)),
    est=c(x$est_group1[which(x$Lin)], x$est_group1[which(!x$Lin & !x$Rin)], x$est_group1[which(x$Rin)],
          x$est_group2[which(x$Lin)], x$est_group2[which(!x$Lin & !x$Rin)], x$est_group2[which(x$Rin)]),
    
    lower=c(x$lower_group1[which(x$Lin)], x$lower_group1[which(!x$Lin & !x$Rin)], x$lower_group1[which(x$Rin)],
            x$lower_group2[which(x$Lin)], x$lower_group2[which(!x$Lin & !x$Rin)], x$lower_group2[which(x$Rin)]),
    
    upper=c(x$upper_group1[which(x$Lin)], x$upper_group1[which(!x$Lin & !x$Rin)], x$upper_group1[which(x$Rin)],
            x$upper_group2[which(x$Lin)], x$upper_group2[which(!x$Lin & !x$Rin)], x$upper_group2[which(x$Rin)])
    
  )
  
  ggplot(x_df, aes(x=x_df$time)) + geom_step(aes(y=x_df$est)) + 
    geom_step(aes(y=x_df$lower),linetype="dashed") + geom_step(aes(y=x_df$upper), linetype="dashed") +
    facet_wrap(vars(x_df$group), ncol=1, labeller=labeller(group=c(group1=x$g1,group2=x$g2))) + labs(y="Survival") +
    theme_bw()
  
  
}

# function to pull out beta estimate and CI for given test time(s) from twosamp object
BtCI<-function(x,tstar,...){
  nt<-length(tstar)
  I<-rep(NA,nt)
  index<-1:length(x$beta)
  ## picki gives TRUE/FALSE vector, TRUE where tval fits 
  ## into interval
  picki<-function(tval){
    (x$L<tval & x$R>tval) | (x$L==tval & x$Lin) | 
      (x$R==tval & x$Rin)
  }
  for (j in 1:nt){
    I[j]<-index[picki(tstar[j])]
  }
  out<-data.frame(time=tstar,beta=x$beta[I],
                  lower=x$lower[I],upper=x$upper[I], p.value=x$p.value[I])

  names(out)[2] <- names(x$null.value)
  attr(out,"conf.level")<- x$conf.level
  out     
}

# matches equation 1 in the supplement, except not transformed back with g yet
calc_eqS1 <- function(S1Est, S2Est, h, phi, Za, S1primeEst, V1Est, S2primeEst, V2Est, dh){
  
  Sigma1 <- calc_sigma(S1primeEst,V1Est, dh)
  Sigma2 <- calc_sigma(S2primeEst,V2Est, dh)
  
  DEst <- h(S2Est)-h(S1Est)
  lower <- DEst - phi*Za*sqrt(Sigma1 + Sigma2)
  upper <- DEst + phi*Za*sqrt(Sigma1 + Sigma2)

list(D=DEst,
     h.ci.lo=lower,
     h.ci.up=upper)
}


# matches equation 2 in the supplement
calc_eqS2 <- function(method, VEst, VarHs=NULL, ShatStar=NULL, n=NULL){
  
  # Find indices where VEst = 0
  ind <- which(VEst==0)
  
  # if length 0, just return VEst
  if(length(ind)==0){
    
    Vdot <- VEst
    
  } else{
  
  # new variance vector
  Vdot <- rep(NA, length(VEst))
  
  # any not 0 will stay the same
  Vdot[-ind] <- VEst[-ind]
  
  methodType <- ifelse(method=="standard","standard","Borkowf")
  
  if(methodType=="standard"){
    Vdot[ind] <- (ShatStar*(1-ShatStar)/n)[ind]
    }
  
  if(methodType=="Borkowf"){
    
    Vdot[ind] <- VarHs[ind]
    
    }
  
  }
  Vdot
}


# function to find replacements for zero-one adjustment based on survival estimates, parmtype, and group

zero.one.adjust <- function(Shat1, Shat2, method,Shatprime_g1, Shatprime_g2, phi, h, Za, VarEst_g1, VarHs_g1, VarEst_g2,VarHs_g2,ShatStar_g1,ShatStar_g2,n1, n2,dh){

  # recall estimates are beta=g(D), where D=h(S2)-h(S1)
  # make new vectors for adjusted D, lower, and upper CI limits
  D_adj <- rep(NA, length(Shat1))
  lower_adj <- rep(NA, length(Shat1))
  upper_adj <- rep(NA, length(Shat1))
  
  
  # if Shat1 and Shat2 are both not = to 0 or 1, no adjustment necessary
  eqS1_unadjusted <- calc_eqS1(S1Est=Shat1, S2Est=Shat2, h=h, phi=phi, Za=Za, S1primeEst=Shatprime_g1, V1Est=VarEst_g1, S2primeEst=Shatprime_g2, V2Est=VarEst_g2, dh=dh)
  D_adj[which(!Shat1 %in% c(0,1) & !Shat2 %in% c(0,1))] <- eqS1_unadjusted$D[which(!Shat1 %in% c(0,1) & !Shat2 %in% c(0,1))]
  lower_adj[which(!Shat1 %in% c(0,1) & !Shat2 %in% c(0,1))] <- eqS1_unadjusted$h.ci.lo[which(!Shat1 %in% c(0,1) & !Shat2 %in% c(0,1))]
  upper_adj[which(!Shat1 %in% c(0,1) & !Shat2 %in% c(0,1))] <- eqS1_unadjusted$h.ci.up[which(!Shat1 %in% c(0,1) & !Shat2 %in% c(0,1))]
  
  # equality value, lower limit, upper limit
  Deq <- 0
  lowerLimit_D<- h(0)-h(1)
  upperLimit_D<- h(1)-h(0)
  
  # First calculate equations S1-S3 with replacements
  
  ## Calculate Vdot whenever V = 0
  Vdot_g1 <- calc_eqS2(method=method, VEst=VarEst_g1, VarHs=VarHs_g1, ShatStar=ShatStar_g1, n=n1)
  Vdot_g2 <- calc_eqS2(method=method, VEst=VarEst_g2, VarHs=VarHs_g2, ShatStar=ShatStar_g2, n=n2)
  
  # Equation S1 with Vdot replacing V
  eqS1_vdot_replace <- calc_eqS1(S1Est=Shat1, S2Est=Shat2, h=h, phi=phi, Za=Za, S1primeEst=Shatprime_g1, V1Est=Vdot_g1, S2primeEst=Shatprime_g2, V2Est=Vdot_g2, dh=dh)
  
  ## replace Shat and Shatprime if h(S) is infinite
  Shat_replace_g1 <- ifelse(is.infinite(h(Shat1)), ShatStar_g1,Shat1)
  Shatprime_replace_g1 <- ifelse(is.infinite(h(Shat1)), ShatStar_g1,Shatprime_g1)
  
  Shat_replace_g2 <- ifelse(is.infinite(h(Shat2)), ShatStar_g2,Shat2)
  Shatprime_replace_g2 <- ifelse(is.infinite(h(Shat2)), ShatStar_g2,Shatprime_g2)
  
  ## equation S1 with Shat and vdot replacements
  eqS1_shat_replace <- calc_eqS1(S1Est=Shat_replace_g1, S2Est=Shat_replace_g2, h=h, phi=phi, Za=Za, S1primeEst=Shatprime_replace_g1, V1Est=Vdot_g1, S2primeEst=Shatprime_replace_g2, V2Est=Vdot_g2, dh=dh)
  
  # Now update D est, lower limit, and upper limit based on each type
  
  # Type 1: S1 = S2 = 0 or 1
  ## get type 1A indices
  ind_1A <- which(Shat1==Shat2 & (Shat1==0 | Shat1==1) & !is.infinite(h(Shat1)))
  
  D_adj[ind_1A] <- Deq
  lower_adj[ind_1A] <- eqS1_vdot_replace$h.ci.lo[ind_1A]
  upper_adj[ind_1A] <- eqS1_vdot_replace$h.ci.up[ind_1A]
  
  ## get type 1B indices
  ind_1B <- which(Shat1==Shat2 & (Shat1==0 | Shat1==1) & is.infinite(h(Shat1)))
  
  D_adj[ind_1B] <- Deq
  lower_adj[ind_1B] <- lowerLimit_D
  upper_adj[ind_1B] <- upperLimit_D
  
  # Type 2: (0,1) and (1,0)

  ## get type 2A indices (0, 1)
  ind_2A <- which(Shat1==0 & Shat2==1)
  
  ## get type 2B indices (1, 0)
  ind_2B <- which(Shat1==1 & Shat2==0)
  
  D_adj[ind_2A] <- upperLimit_D
  lower_adj[ind_2A] <- eqS1_shat_replace$h.ci.lo[ind_2A]
  upper_adj[ind_2A] <- upperLimit_D
  
  D_adj[ind_2B] <- lowerLimit_D
  lower_adj[ind_2B] <- lowerLimit_D  
  upper_adj[ind_2B] <- eqS1_shat_replace$h.ci.up[ind_2B]
  
  # Type 3: (0,s) or (s,0) or (1,s) or (s,1)
  ind_3 <- which( ((Shat1==0 | Shat1==1) & !(Shat2==0 | Shat2==1) ) | 
                    (Shat2==0 | Shat2==1) & !(Shat1==0 | Shat1==1))
  
  D_adj[ind_3] <- (h(Shat2)-h(Shat1))[ind_3]
  
  # calculate lower and upper limits
  lower_adj[ind_3] <- ifelse(D_adj[ind_3]==lowerLimit_D, lowerLimit_D, eqS1_shat_replace$h.ci.lo[ind_3])
  
  # calculate upper limit
  upper_adj[ind_3] <- ifelse(D_adj[ind_3]==upperLimit_D, upperLimit_D, eqS1_shat_replace$h.ci.up[ind_3])
  
  
  out <- list(D_adj=D_adj,
              lower_adj=lower_adj,
              upper_adj=upper_adj)
  
  out
  
  
}

get_h <- function(parmtype=NULL, trans=NULL){
  
  if(is.null(trans)){
    trans <- parmtype_to_trans(parmtype)
  }
  
  h <- switch(trans,
              none=function(S){S},
              log=function(S){ log(S)},
              logodds=function(S){ log(S/(1-S))},
              loglog=function(S){log(-log(S))},
              clog=function(S){log(1-S)},
              cloglog=function(S){ log(-log(1-S)) }
              )
  
  return(h)
}

get_dh <- function(parmtype=NULL, trans=NULL){
  
  if(is.null(trans)){
    trans <- parmtype_to_trans(parmtype)
  }
  # some transformations have ifelse statements because the 
  # the 0  or 1 survival cases should equal their limits 
  # from the correct direction, but using plug-in estimates 
  # does not equal those limits
  dh <- switch(trans,
              none= function(S){1},
              log=function(S){1/S},
              logodds=function(S){ 1/(S*(1-S)) },
              loglog=function(S){ifelse(S==0 | S==1, -Inf, 1/(S*log(S)))},
              clog=function(S){-1/(1-S)},
              cloglog=function(S){ifelse(S==0 | S==1,Inf, 1/((1-S)*(-log(1-S)))) }
  )
  
  return(dh)
}

get_hinv <- function(trans){
  
  hinv <- switch(trans,
                 none=function(y){y} ,
                 logodds=function(y){ ifelse(y>100,1,exp(y)/(1+exp(y))) }  ,
                 log=function(y){exp(y)},
                 loglog=function(y){exp(-exp(y))}  ,
                 clog=function(y){1-exp(y)} ,
                 cloglog=function(y){1-exp(-exp(y))} )
  
  return(hinv)
}

get_g <- function(parmtype){
  
  g <- switch(parmtype,
               difference=function(D){D},
               oddsratio=function(D){ exp(D) }  ,
               ratio=function(D){exp(D)},
               efflogs=function(D){1-exp(D)},
               effcdf=function(D){1-exp(D)},
               cdfratio=function(D){1/exp(D)},
              logsratio=function(D){ 
                out <- rep(NA, length(D))
                out[which(D==Inf)] <- 1/Inf
                out[which(D==-Inf)] <- Inf
                out[which(!is.infinite(D))] <- 1/exp(D[which(!is.infinite(D))])
                
                out
                })
  
  return(g)
}


get_ginv <- function(parmtype){
  
  ginv <- switch(parmtype,
              difference=function(beta){beta},
              oddsratio=function(beta){ log(beta) }  ,
              ratio=function(beta){log(beta)},
              efflogs=function(beta){log(1-beta)},
              effcdf=function(beta){log(1-beta)},
              cdfratio=function(beta){-log(beta)},
              logsratio=function(beta){-log(beta)
              })
  
  return(ginv)
}

get_all_times <- function(time, status, cens_symbol=TRUE){
  # get unique event/censoring times
  event_times <- unique(time[which(status==1)])
  cens_times <- unique(time[which(status==0)])
  sorting_i <- sort(c(event_times,cens_times), index.return=T)
  
  if(cens_symbol){
    all_times <- c(event_times,paste0(cens_times,"+"))[sorting_i$ix]
  } else{
    all_times <- c(event_times,cens_times)[sorting_i$ix]
  }
  
  return(all_times)
}

delta2samp <- function(time,status,group, testtime=NULL,conf.level=0.95,
                       zero.one.adjustment=FALSE,method=c("standard","reg_hybrid","adj_hybrid","sh_adj_hybrid"),
                       parmtype=c("difference","oddsratio","ratio","efflogs","effcdf"),
                       nullparm=NULL,
                       alternative=c("two.sided","less","greater"),changeGroupOrder=FALSE){
  
  
  # run two-sample argument checks and switches based on parmtype
  check_args <- twosampleChecks(time=time, testtime=testtime, conf.level=conf.level,status=status, group=group, alternative=alternative, parmtype=parmtype, method=method,nullparm=nullparm, changeGroupOrder=changeGroupOrder, method_type="delta")
  method <- check_args$meth
  
  all_times <- get_all_times(time=time, status=status)

  # next we will run delta.calc on each group separately
  time_g1 <- time[which(group==check_args$ug[1])]
  status_g1 <- status[which(group==check_args$ug[1])]
  delta_g1 <- delta.calc(all_times=all_times,
                         time_g=time_g1,
                         status_g=status_g1,
                         trans=check_args$trans, method=method, alpha=check_args$alpha,
                         zero.one.adjustment=zero.one.adjustment,two_samp_out=T)
  
  time_g2 <- time[which(group==check_args$ug[2])]
  status_g2 <- status[which(group==check_args$ug[2])]
  delta_g2 <- delta.calc(all_times=all_times,
                         time_g=time_g2,
                         status_g=status_g2,
                         trans=check_args$trans, method=method, alpha=check_args$alpha,
                         zero.one.adjustment=zero.one.adjustment,two_samp_out=T)
  
  
  # transformation functions based on parmtype
  h <-get_h(trans=check_args$trans)
  dh <-get_dh(trans=check_args$trans)
  g <-get_g(parmtype=check_args$ptype)
  ginv <- get_ginv(parmtype=check_args$ptype)
  
  # For each group, calculate h(S) and Sigma(S*, V) where S, S*, and V are based on the delta method
  if(method=="standard"){
    Surv_g1 <- delta_g1$KM
    Surv_g2 <- delta_g2$KM
    
    Shatprime_g1 <- delta_g1$KM
    VarEst_g1 <- delta_g1$gw
    
    Shatprime_g2 <- delta_g2$KM
    VarEst_g2 <- delta_g2$gw
  
    
  } else if (method=="reg_hybrid"){
    
    Surv_g1 <- delta_g1$KM
    Surv_g2 <- delta_g2$KM
    
    Shatprime_g1 <- delta_g1$KM
    VarEst_g1 <- delta_g1$VarH
    
    Shatprime_g2 <- delta_g2$KM
    VarEst_g2 <- delta_g2$VarH
   
  } else if (method=="adj_hybrid"){
    Surv_g1 <- delta_g1$KM
    Surv_g2 <- delta_g2$KM
    
    Shatprime_g1 <- delta_g1$KM
    VarEst_g1 <- delta_g1$VarHs
    
    Shatprime_g2 <- delta_g2$KM
    VarEst_g2 <- delta_g2$VarHs
    
  } else{
    
    Surv_g1 <- delta_g1$KMs
    Surv_g2 <- delta_g2$KMs
    
    Shatprime_g1 <- delta_g1$KM
    VarEst_g1 <- delta_g1$VarHs
    
    Shatprime_g2 <- delta_g2$KM
    VarEst_g2 <- delta_g2$VarHs
  }
  
  # calculate upper and lower CIs
  phi <- ifelse(check_args$ptype=="efflogs" | check_args$ptype=="effcdf", -1, 1)
  Za<-ifelse(check_args$alt=="two.sided",qnorm(1-check_args$alpha/2),
             qnorm(1-check_args$alpha))
  
  # need to calculate confidence interval differently if S1 or S2 = 0 or 1
  if(zero.one.adjustment){
    
   adjusted_res <- zero.one.adjust(
      Shat1=Surv_g1, Shat2=Surv_g2, 
      method=method, h=h, Za=Za, phi=phi,
      Shatprime_g1=Shatprime_g1, Shatprime_g2=Shatprime_g2,  
      VarEst_g1=VarEst_g1, VarEst_g2=VarEst_g2, 
      VarHs_g1=delta_g1$VarHs,VarHs_g2=delta_g2$VarHs,
      ShatStar_g1=delta_g1$KMs, ShatStar_g2=delta_g2$KMs, 
      n1=delta_g1$ni[1], n2=delta_g2$ni[1],
      dh=dh
    )
   
    
   D <- adjusted_res$D_adj
   h.ci.lo <- adjusted_res$lower_adj
   h.ci.up <- adjusted_res$upper_adj

   # calculate z_score and p-value
   # MISTAKE??!??   SE <- (h.ci.up - h.ci.lo)/Za
   SE<- (h.ci.up - h.ci.lo)/(2*phi*Za)
   #if (SE<=0) stop("fix code, the SE is not positive")
   # want a z-score that will give a p-value of 1 if SE blows up
   z_score <- ifelse(is.infinite(SE),
                     ifelse(check_args$alt=="two.sided",0,
                            ifelse(check_args$alt=="less", Inf,-Inf)),
                     phi*(D-ginv(check_args$Nullparm))/SE
                     
                     )
   
   
    
  } else{
    # no zero.one.adjustment
    
    D <- h(Surv_g2) - h(Surv_g1)
    
    Sigma_g1 <- calc_sigma(Shatprime_g1, VarEst_g1, dh)
    Sigma_g2 <- calc_sigma(Shatprime_g2, VarEst_g2, dh)
    
    h.ci.lo <- D - phi*Za*sqrt(Sigma_g1 + Sigma_g2)
    h.ci.up <- D + phi*Za*sqrt(Sigma_g1 + Sigma_g2)
    
    # calculate z-score and p-value
    z_score <- phi*(D - ginv(check_args$Nullparm) )/ sqrt(Sigma_g1 + Sigma_g2)
  
  }
  

  
  
  
  
  # calculate p-value
  if(check_args$alt=="two.sided"){
    p.value <- 2*(1-pnorm(abs(z_score)))
  } else if(check_args$alt=="greater"){
    p.value <- 1-pnorm(z_score)
  } else{
    p.value <- pnorm(z_score)
  }
  
  # transform to get D estimate and lower/upper conf int.
  beta <- g(D)
  g.h.ci.lo <- g(h.ci.lo)
  g.h.ci.up <- g(h.ci.up)
  # for one-sided tests, we only keep one side of the confidence limits
  if (check_args$alt=="less"){
    # if the alternative is: parm< Nullparm
    # then only the upper confidence limit matters 
    # set lower limit to its extreme
    g.h.ci.lo <- rep(check_args$lowerLimit,length(g.h.ci.up))
  } else if (check_args$alt=="greater"){
    # if the alternative is: parm> Nullparm
    # then only the lower confidence limit matters 
    # set upper limit to its extreme 
    g.h.ci.up <- rep(check_args$upperLimit,length(g.h.ci.lo) )
  }
  
  
  # create interval from L, Lin, Rz, Rin
  Interval<-intChar(delta_g1$L,delta_g1$R,delta_g1$Lin,delta_g1$Rin)
  
  # if test time is provided, then class htest for output
  if(!is.null(testtime)){
    # find index of interval that includes testtime
    # if testtime is beginning or end of an interval (L or R), select inclusive interval
    if(testtime %in% delta_g1$L | testtime %in% delta_g1$R){
      
      row_i <- which( (testtime == delta_g1$L & delta_g1$Lin) | (testtime == delta_g1$R & delta_g1$Rin))
      
    } else{
      # find interval that contains testtime
      row_i <- which(testtime > delta_g1$L & testtime < delta_g1$R)
    }
    
    out <- create.htest(testtime = testtime, alt=check_args$alt,conf.level=conf.level, lower=g.h.ci.lo[row_i], upper=g.h.ci.up[row_i], p.value=p.value[row_i], est1=Surv_g1[row_i], est2=Surv_g2[row_i], beta=beta[row_i],nullparm=check_args$Nullparm, ptype=check_args$ptype, ug=check_args$ug, method=method,zero.one.adjustment=zero.one.adjustment)
    
  } else{
  
  # else provide interval, beta est, CI, and p-value, as well as one sample
    

      
  
  out <- create.twosamp(alt=check_args$alt,
                        conf.level=conf.level,
                        L=delta_g1$L, Lin=delta_g1$Lin, R=delta_g1$R, Rin=delta_g1$Rin, 
                        interval=Interval,
                        g1=check_args$ug[1],
                        est_g1=Surv_g1,
                        lower_g1=delta_g1$lower,
                        upper_g1=delta_g1$upper,
                        g2=check_args$ug[2],
                        est_g2=Surv_g2,
                        lower_g2=delta_g2$lower,
                        upper_g2=delta_g2$upper,
                        beta=beta,
                        lower=g.h.ci.lo,
                        upper=g.h.ci.up,
                        p.value=p.value,
                        nullparm=check_args$Nullparm, 
                        ptype=check_args$ptype, 
                        method=method,
                        zero.one.adjustment=zero.one.adjustment
    )
  
  }
  
  return(out)
}