R/nSurvNPH.R

In keaven/nphsim: Non proportional hazards sample size and simulation

#' Sample Size Under Non-Propotional Hazards Assumptions
#'
#' Description text; note that parameter descriptions below as well as text for 
#' help file needs to be updated.
#'
#' Details text
#'
#' @param lambdaC scalar, vector or matrix of event hazard rates for the 
#' control group; rows represent time periods while columns represent strata; a vector implies a single stratum.
#' @param hr hazard ratio (experimental/control) under the alternate hypothesis (scalar)
#' @param hr0 hazard ratio (experimental/control) under the null hypothesis (scalar)
#' @param etaC scalar, vector or matrix of dropout hazard rates for the control group; rows represent time periods while columns represent strata; if entered as a scalar, rate is constant accross strata and time periods; 
#' if entered as a vector, rates are constant accross strata
#' @param etaE matrix dropout hazard rates for the experimental group specified 
#' in like form as eta; if NULL, this is set equal to eta
#' @param NONCMPL TBD
#' @param DROPIN TBD
#' @param gamma a scalar, vector or matrix of rates of entry by time period (rows) 
#' and strata (columns); if entered as a scalar, rate is constant accross strata 
#' and time periods; if entered as a vector, rates are constant accross strata
#' @param R a scalar or vector of durations of time periods for recruitment rates 
#' specified in rows of gamma. Length is the same as number of rows in gamma. 
#' Note that when variable enrollment duration is specified (input T=NULL), 
#' the final enrollment period is extended as long as needed
#' @param S a scalar or vector of durations of piecewise constant event rates 
#' specified in rows of lambda, eta and etaE; this is NULL if there is a single 
#' event rate per stratum (exponential failure) or length of the number of rows 
#' in lambda minus 1, otherwise
#' @param T study duration; if T is input as NULL, this will be computed on output; 
#' see details.
#' @param minfup follow-up of last patient enrolled; if minfup is input as NULL, 
#' this will be computed on output; see details
#' @param ratio randomization ratio of experimental treatment divided by control; 
#' normally a scalar, but may be a vector with length equal to number of strata.
#' @param alpha type I error rate. Default is 0.025 since 1-sided testing is default
#' @param beta type II error rate. Default is 0.10 (90% power); 
#' NULL if power is to be computed based on other input values
#' @param sided 1 for 1-sided testing, 2 for 2-sided testing
#' @param tol or cases when T or minfup values are derived through root finding (T or minfup input as NULL), 
#' tol provides the level of error input to the uniroot() root-finding function. The default is the same as for uniroot.
#' @param SBDV This is a tuning parameter for the Lakatos method specifying how many intervals to use for
#' the Markov approximation applied in the method
#' @param SIMULT TBD
#' @param TRT_LAG lag treatment
#' @param nstrata TBD
#' @examples 
#' #no staggered entry, no lag time;
#' #no staggered entry, lag time;
#' #stagger entry, no lag time
#' #stagged entry; lag time;
#' #strata;
#' 
#' @export
nSurvNPH<-function ( lambdaC=log(2)/6
                    ,hr = 0.6
                    ,hr0 = 1
                    ,etaC = 0
                    ,etaE = NULL
                    ,NONCMPL=0
                    ,DROPIN= 0
                    ,gamma = 1
                    ,R = 12
                    ,S = NULL
                    ,T = NULL
                    ,minfup = NULL
                    ,ratio = 1
                    ,alpha = 0.025
                    ,beta = 0.1
                    ,sided = 1, tol = .Machine$double.eps^0.25
                    ,SBDV=20
                    ,SIMULT=FALSE
                    ,TRT_LAG=0 #lag treatment#
                    ,nstrata=1)

{
  zalpha <- -qnorm(alpha/sided)
  zbeta<--qnorm(beta)

  Qe<-ratio/(1+ratio) #proprotion of Experimental group
  Qc<-1-Qe #proportion of Control group

  N_INTRVL <- SBDV
  if (TRT_LAG>0)
  { a<-SBDV*TRT_LAG
    if (abs(a-round(a))>1e-6) stop("SBDV*TRT_LAG should be integer")
    NACTV    <- round(a)
    NSTATES  <- (2 * NACTV) + 2;
  }

  NN<-round(T*N_INTRVL) #total number of time points in Markov process
  #adjust gamma and recrutiment period given R and T;
  if (!is.matrix(gamma)) gamma <- matrix(gamma)
  cumR_org<-cumsum(R)
  cumR_adj<-unique(c(0,cumR_org[cumR_org<T],min(T,sum(R))))
  R<-diff(cumR_adj)
  nper<-length(R) #no. of recruitment period
  gamma<-matrix(gamma[1:nper,],ncol=nstrata)
  mat0<-matrix(0,NN,nstrata)
  gamma_NN<-mat0
  temp<-apply(gamma,2,FUN=rep,times=R*N_INTRVL)/N_INTRVL
  gamma_NN[1:nrow(temp),]<-temp
  RCRT_CUM<-matrix(apply(gamma_NN,2,cumsum),ncol=nstrata)
  AD_CENS<-mat0
  AD_CENS[NN:1,]<-gamma_NN/RCRT_CUM; #administrative censoring;

  if (is.null(S)) {nper<-1;S<-T}
  else nper<-length(S)+1
  if (!is.matrix(lambdaC)) lambdaC <- matrix(lambdaC, nrow=nper, ncol=nstrata) #row represent time, col represent strata
  if (!is.matrix(NONCMPL)) NONCMPL <- matrix(NONCMPL, nrow=nper, ncol=nstrata)
  if (!is.matrix(DROPIN)) DROPIN <- matrix(DROPIN, nrow=nper, ncol=nstrata)
  if (!is.matrix(hr)) hr <- matrix(hr,nrow=nper,ncol=nstrata)
  if (is.null(etaE)) etaE<-etaC
  if (!is.matrix(etaC)) etaC <- matrix(etaC,nrow=nper,ncol=nstrata)
  if (!is.matrix(etaE)) etaE <- matrix(etaE,nrow=nper,ncol=nstrata)

  #adjust event rate, drop out rate given S and T;
  cumS_org<-cumsum(S)
  cumS_adj<-unique(c(0,cumS_org[cumS_org<T],T))
  S<-diff(cumS_adj)
  lambdaE<-hr*lambdaC

  SN<-S*N_INTRVL
  nper<-length(S)

  NONCMPL_NP<-1-(1-NONCMPL[1:nper,,drop=FALSE])^(1/(N_INTRVL))
  DROPIN_NP<-1-(1-DROPIN[1:nper,,drop=FALSE])^(1/(N_INTRVL))
  etaC_NP<-1-(1-etaC[1:nper,,drop=FALSE])^(1/(N_INTRVL))
  etaE_NP<-1-(1-etaE[1:nper,,drop=FALSE])^(1/(N_INTRVL))
  pC_NP<-1-exp(-lambdaC[1:nper,,drop=FALSE])^(1/N_INTRVL)
  pE_NP<-1-exp(-lambdaE[1:nper,,drop=FALSE])^(1/N_INTRVL)


  NONCMPL_NN<-apply(NONCMPL_NP,2,FUN=rep,times=SN)
  DROPIN_NN<- apply(DROPIN_NP,2,FUN=rep,times=SN)
  etaC_NN<-apply(etaC_NP,2,FUN=rep,times=SN)#drop out rate for control;
  etaE_NN<-apply(etaE_NP,2,FUN=rep,times=SN)#drop out rate for experiment;
  pC_NN<-  apply(pC_NP,2,FUN=rep,times=SN)  #event rate for control;
  pE_NN<-  apply(pE_NP,2,FUN=rep,times=SN)  #event rate for experiment;

#MARKOV;
#INITAILIZE MATRICES;
DSTR<-vector("list", nstrata)


if (!TRT_LAG)
{ P_C<-rep(0,nstrata)
  P_E<-rep(0,nstrata)
  P_ALL<-rep(0,nstrata)
  eD_LR<-rep(0,nstrata)
  for (strata in 1:nstrata)
  {
    #initial state
  DISTR_E<-matrix(c(0,0,1,0),ncol=1)
  DISTR_C<-matrix(c(0,0,0,1),ncol=1)
  TRANS<-matrix(0,4,4)
  DSTR_E<-matrix(0,4,NN)
  DSTR_C<-matrix(0,4,NN)
  rownames(DSTR_C)<-c('LOSSES', 'EVENTS', 'DROPIN', 'ACTV C')
  rownames(DSTR_E)<-c('LOSSES', 'EVENTS', 'ACTV_E', 'NONCMPL')

   #transition matrix
    for ( I in 1:NN)
    { TRANS<-diag(4);
      TRANS[,3]<-c(etaE_NN[I,strata],pE_NN[I,strata],0,NONCMPL_NN[I,strata]);
      TRANS[,4]<-c(etaC_NN[I,strata],pC_NN[I,strata],DROPIN_NN[I,strata],0);
      diag(TRANS)<-0
      diag(TRANS)<-1-colSums(TRANS);
      DISTR_E<-TRANS%*%DISTR_E;
      DISTR_C<-TRANS%*%DISTR_C;

      if (!SIMULT) #start addressing stagged entry
        {
        TEMP_E<-DISTR_E[c(3,4)]*(1-AD_CENS[I,strata]);
        DISTR_E[1]<-DISTR_E[1]+sum(DISTR_E[c(3, 4)]-TEMP_E);
        DISTR_E[c(3,4)]<-TEMP_E;
        TEMP_C<-DISTR_C[c(3,4)]*(1-AD_CENS[I,strata]);
        DISTR_C[1]<-DISTR_C[1]+sum(DISTR_C[c(3, 4)]-TEMP_C);
        DISTR_C[c(3,4)]<-TEMP_C;
        }#end of addressing stagged entry


      DSTR_E[,I]<-DISTR_E;
      DSTR_C[,I]<-DISTR_C;
      }

  EVENT_C<-diff(c(0,DSTR_C[2,]))
  EVENT_E<-diff(c(0,DSTR_E[2,]))
  EVENT_ALL<-Qc*EVENT_C+Qe*EVENT_E
  LOSS_C<-diff(c(0,DSTR_C[1,]))
  LOSS_E<-diff(c(0,DSTR_C[1,]))
  ATRISK_C<-1-DSTR_C[1,]-DSTR_C[2,]+EVENT_C+LOSS_C
  ATRISK_E<-1-DSTR_E[1,]-DSTR_E[2,]+EVENT_E+LOSS_E
  PHI<-Qc/Qe*ATRISK_C/ATRISK_E #population ratio
  THETA<-log(1-EVENT_C/ATRISK_C)/log(1-EVENT_E/ATRISK_E) #risk ratio
  RHO<-EVENT_ALL/sum(EVENT_ALL)  # di/d
  GAMMA<-PHI*THETA/(1+PHI*THETA)-PHI/(1+PHI);
  ETA<-PHI/(1+PHI)^2;
  P_E[strata]<-DSTR_E[2,ncol(DSTR_E)]
  P_C[strata]<-DSTR_C[2,ncol(DSTR_C)]
  P_ALL[strata]<-Qc*P_C[strata]+Qe*P_E[strata]
  eD_LR[strata]<-sum(RHO*GAMMA)/sqrt(sum(RHO*ETA))
  DSTR[[strata]]<-list(DSTR_E=DSTR_E, DSTR_C=DSTR_C)

  }

} #end of no treatment lag

else if (TRT_LAG>0)
  {
  lambdaC<-lambdaC[1,1,drop=FALSE]
  lambdaE<-lambdaE[1,1,drop=FALSE]
  P_C<-rep(0,nstrata)
  P_E<-rep(0,nstrata)
  P_ALL<-rep(0,nstrata)
  eD_LR<-rep(0,nstrata)
  for (strata in 1:nstrata)
  {
    A<-matrix(0,NACTV,NACTV)
    B<-matrix(0,NACTV,NACTV)
    C<-matrix(0,NACTV,NACTV)
    D<-matrix(0,NACTV,NACTV)

    #initial state
    DISTR_E<-matrix(0,nrow=NSTATES,ncol=1)
    DISTR_C<-matrix(0,nrow=NSTATES,ncol=1)
    DISTR_E[3]<-1
    DISTR_C[NACTV+3]<-1
    TRANS<-matrix(0,NSTATES,NSTATES)
    DSTR_E<-matrix(0,4,NN)
    DSTR_C<-matrix(0,4,NN)
    rownames(DSTR_C)<-c('LOSSES', 'EVENTS', 'DROPIN', 'ACTV_C')
    rownames(DSTR_E)<-c('LOSSES', 'EVENTS', 'ACTV_E', 'NONCMPL')


    #transition matrix
    for ( I in 1:NN)
    {TRANS<-matrix(0,NSTATES,NSTATES);
    A<-matrix(0,NACTV,NACTV)
    B<-matrix(0,NACTV,NACTV)
    C<-matrix(0,NACTV,NACTV)
    D<-matrix(0,NACTV,NACTV)
    AA<-1-exp(-lambdaC-(0:NACTV)*(lambdaE-lambdaC)/NACTV)^(1/N_INTRVL)
    TRANS[1,]<-c(1,0,rep(etaE_NN[I,strata],NACTV),rep(etaC_NN[I,strata],NACTV))
    TRANS[2,]<-c(0,1,AA[2:(NACTV+1)],AA[1:(NACTV)])

    C<-diag(rep(DROPIN_NN[I,strata],NACTV))
    B<-diag(rep(NONCMPL_NN[I,strata],NACTV))

    if (NACTV>2) {
    diag(A[2:NACTV,1:(NACTV-1)])<-1-colSums(TRANS[,3:(NACTV+1),drop=FALSE])-NONCMPL_NN[I,strata]
    diag(D[1:(NACTV-1),2:(NACTV)])<-1-colSums(TRANS[,(NACTV+4):(2*NACTV+2),drop=FALSE])-DROPIN_NN[I,strata]
    }
    else
    {A[2,1]<-1-sum(TRANS[,3])-NONCMPL_NN[I,strata]
     D[1,2]<-1-sum(TRANS[,6])-DROPIN_NN[I,strata]
     }
    A[NACTV,NACTV]<-1-colSums(TRANS[,(NACTV+2),drop=FALSE])-NONCMPL_NN[I,strata]


    D[1,1]<-1-colSums(TRANS[,(NACTV+3),drop=FALSE])-DROPIN_NN[I,strata]
    TRANS[(3:NSTATES),(3:NSTATES)]<-cbind(rbind(A,B),rbind(C,D))
    DISTR_E<-TRANS%*%DISTR_E;
    DISTR_C<-TRANS%*%DISTR_C;

    if (!SIMULT) #start addressing stagged entry
    {
      TEMP_E<-DISTR_E[c(3:NSTATES)]*(1-AD_CENS[I,strata]);
      DISTR_E[1]<-DISTR_E[1]+sum(DISTR_E[c(3:NSTATES)]-TEMP_E);
      DISTR_E[3:NSTATES]<-TEMP_E;
      TEMP_C<-DISTR_C[3:NSTATES]*(1-AD_CENS[I,strata]);
      DISTR_C[1]<-DISTR_C[1]+sum(DISTR_C[3:NSTATES]-TEMP_C);
      DISTR_C[3:NSTATES]<-TEMP_C;
    }#end of addressing stagged entry


    DSTR_E[,I]<-c(DISTR_E[1:2],sum(DISTR_E[3:(NACTV+2)]),sum(DISTR_E[(NACTV+3):NSTATES]))
    DSTR_C[,I]<-c(DISTR_C[1:2],sum(DISTR_C[3:(NACTV+2)]),sum(DISTR_C[(NACTV+3):NSTATES]))
    }

    EVENT_C<-diff(c(0,DSTR_C[2,]))
    EVENT_E<-diff(c(0,DSTR_E[2,]))
    EVENT_ALL<-Qc*EVENT_C+Qe*EVENT_E
    LOSS_C<-diff(c(0,DSTR_C[1,]))
    LOSS_E<-diff(c(0,DSTR_C[1,]))
    ATRISK_C<-1-DSTR_C[1,]-DSTR_C[2,]+EVENT_C+LOSS_C
    ATRISK_E<-1-DSTR_E[1,]-DSTR_E[2,]+EVENT_E+LOSS_E
    PHI<-Qc/Qe*ATRISK_C/ATRISK_E #population ratio
    THETA<-log(1-EVENT_C/ATRISK_C)/log(1-EVENT_E/ATRISK_E) #risk ratio
    RHO<-EVENT_ALL/sum(EVENT_ALL)  # di/d
    GAMMA<-PHI*THETA/(1+PHI*THETA)-PHI/(1+PHI);
    ETA<-PHI/(1+PHI)^2;
    P_E[strata]<-DSTR_E[2,ncol(DSTR_E)]
    P_C[strata]<-DSTR_C[2,ncol(DSTR_C)]
    P_ALL[strata]<-Qc*P_C[strata]+Qe*P_E[strata]
    eD_LR[strata]<-sum(RHO*GAMMA)/sqrt(sum(RHO*ETA))
    DSTR[[strata]]<-list(DSTR_E=DSTR_E, DSTR_C=DSTR_C)

  }




  }  #end of treatment lag
q_strata<-colSums(gamma_NN)/sum(gamma_NN)
w_strata<-q_strata*P_E*P_C/P_ALL
w_strata<-w_strata/sum(w_strata)
dj<-q_strata*P_ALL/sum(q_strata*P_ALL)#proportion death for each stratum
C1<-sum(w_strata*eD_LR*dj^0.5)
d<-(zalpha+zbeta)^2/C1^2*sum(w_strata^2)
n<-d/sum(q_strata*P_ALL)

eNC<-n*q_strata*Qc
eNE<-n*q_strata*Qe
eN<-n*q_strata
eDC<-eNC*P_C
eDE<-eNE*P_E
eD<-eN*P_ALL

output<-list(DSTR=DSTR,alpha=alpha, beta=beta
             , eDC=eDC,eDE=eDE,d=eD
             , eNC=eNC,eNE=eNE,n=eN
             , pE=P_E, pC=P_C, p=P_ALL
             , d=round(d),n=round(n))
return(output)

}
#'
#' @importFrom(gsDesign)


#no staggered entry, no lag time;

#no staggered entry, lag time;
#stagger entry, no lag time
#stagged entry; lag time;
#strata;
#undebug("nSurvNPH")
#source("nSurvNPHDebug.R")