R/FUN-cstat-ver003b.R

Defines functions cvC kmcens VEC2MAT Est.PH Inf.Cval.Delta Inf.Cval Est.Cval

Documented in cvC Est.Cval Est.PH Inf.Cval Inf.Cval.Delta kmcens VEC2MAT

#-- required package --#
# library(survival) ; 

#######################################################
# Function: Est.Cval: return point estimte
#           ver.002  --- 2011.2.15 --  can handle ties 
#           ver.003  --- 2011.9.21 --  FORTRAN (conc)
#           ver.003b --- 2013.2.13 --  add nofit option
#######################################################
Est.Cval<-function(mydata, tau, nofit=FALSE){

    ## =============== ##
    ## Weight          ##
    ## =============== ##
    cens <- kmcens(mydata[,1],mydata[,2], tau)
    GXi<-cens$surv[match(mydata[,1], cens$distinct, nomatch=1)]
    Wi<-1/GXi/GXi*mydata[,2]*as.numeric(mydata[,1]<tau) 

    ## ======================== ##
    ## Fit Cox if nofit==FALSE  ##
    ## ======================== ##

    if(nofit){
       rs = mydata[,3]
    }else{
       fit.cox <- Est.PH(mydata)
       rs=as.vector(fit.cox$rs)
    }
 
    ## =============== ##
    ## C-stat (D hat)  ##
    ## =============== ##
    cstat=conc(mydata[,1],mydata[,2],Wi,rs)
    
    ## =============== ##
    ## OUTPUT          ##
    ## =============== ##
    Z=list()
    Z$Dhat=cstat
    Z$rs=rs
    if(nofit==FALSE){
     Z$beta=fit.cox$beta
     Z$beta.var=fit.cox$var
     Z$rs=fit.cox$rs
     Z$Ui=fit.cox$Ui 
     Z$ft=fit.cox$ft
    }
    Z$cens.surv=cens$surv
    Z$cens.psii=cens$psii 
    Z$distinct=cens$distinct
    Z$wt=Wi
    Z$nofit=nofit

    return(Z)

}



#######################################################
# Function: Inf.Cval: give CI
#     ver.003 --- 2011.9.21 --  FORTRAN (conc, unoU2P)
#     ver.003b -- 2013.2.13 --  seed = null (default)
#######################################################
Inf.Cval<-function(mydata, tau, itr=1000, seed=NULL){

	if(!is.null(seed)) {set.seed(seed)}
	
	emp<-Est.Cval(mydata, tau)
	n<-nrow(mydata); p<-ncol(mydata)-2 ; 

	#-------------
	A<-as.matrix(emp$beta.var)
    distinct<-emp$distinct
    cens.surv<-emp$cens.surv

	#-------------
	time<-mydata[,1] ; status<-mydata[,2] ; covs<-mydata[,-1:-2] ;
    UJ<-emp$Ui
    cens.psii<-emp$cens.psii
    Wi<-emp$wt

	#-------------
	temp<-rep(0,itr)
#    Wa.stock<-Wb.stock<-WG.stock<-c()
#    b.star.stock<-c()
#    G.star.stock<-c()

	for (k in 1: itr){

 	    #----- random number ---#
        xi<-rexp(n)

        b.star = emp$beta + n*A%*%unoU2P(emp$Ui,xi) 
        rs.star<-as.matrix(covs)%*%as.vector(b.star)

        G.wk<-rep(0,length(distinct))
        G.wk=unoU2P(cens.psii,xi)
        G.star <- cens.surv - G.wk*cens.surv
        GXi<-G.star[match(mydata[,1], distinct, nomatch=1)]
        Wi.star<-1/GXi/GXi*mydata[,2]*as.numeric(mydata[,1]<tau)

        CW=unoCW(time, status, Wi, Wi.star, emp$rs, rs.star, xi, emp$Dhat)

        #----- result  ---#
	    temp[k] <-  (CW$Wa + CW$Wb + CW$Wg)*sqrt(n)

        #----- for check ----#
        #     b.star.stock<-rbind(b.star.stock, b.star)
        #     G.star.stock<-rbind(G.star.stock, G.star)
        #     Wa.stock<-c(Wa.stock, CW$Wa)
        #     Wb.stock<-c(Wb.stock, CW$Wb)
        #     WG.stock<-c(WG.stock, CW$Wg)
	}
	se=sqrt(var(temp)/n)
	low95<-emp$Dhat - 1.959*se
	upp95<-emp$Dhat + 1.959*se
	# return(list(Dhat=emp$Dhat, low95=low95, upp95=upp95, out=temp, b.star.stock=b.star.stock, G.star.stock=G.star.stock, Wa=Wa.stock, Wb=Wb.stock, WG=WG.stock, se=se))
	return(list(Dhat=emp$Dhat, low95=low95, upp95=upp95, se=se, ft=emp$ft))
}


#######################################################
# Function: Inf.Cval.Delta: give Delta C and its CI
#     ver.003b -- 2013.2.13 --  seed = null (default)
#######################################################
Inf.Cval.Delta<-function(mydata, covs0, covs1, tau, itr=1000, seed=NULL){
   
   mydata1=data.frame(mydata, covs1)
   mydata0=data.frame(mydata, covs0)


	if(!is.null(seed)) {set.seed(seed)}

	n<-nrow(mydata1); p<-ncol(mydata1)-2 ; 

	emp1<-Est.Cval(mydata1, tau)
	emp0<-Est.Cval(mydata0, tau)
	
	#-------------
    distinct<-emp1$distinct
    cens.surv<-emp1$cens.surv

	A1<-as.matrix(emp1$beta.var)
	A0<-as.matrix(emp0$beta.var)

	#-------------
	time<-mydata1[,1] ; status<-mydata1[,2] ; 
	covs1<-mydata1[,-1:-2] ;
	covs0<-mydata0[,-1:-2] ;

    cens.psii<-emp1$cens.psii
    Wi<-emp1$wt

    UJ1<-emp1$Ui
    UJ0<-emp0$Ui

	#-------------
	temp=temp1=temp0=rep(0,itr)

	for (k in 1: itr){


        xi<-rexp(n)
        
        b1.star = emp1$beta + n*A1%*%unoU2P(emp1$Ui,xi) 
        rs1.star<-as.matrix(covs1)%*%as.vector(b1.star)

        b0.star = emp0$beta + n*A0%*%unoU2P(emp0$Ui,xi) 
        rs0.star<-as.matrix(covs0)%*%as.vector(b0.star)

        G.wk<-rep(0,length(distinct))
        G.wk=unoU2P(cens.psii,xi)
        G.star <- cens.surv - G.wk*cens.surv
        GXi<-G.star[match(mydata1[,1], distinct, nomatch=1)]
        Wi.star<-1/GXi/GXi*mydata1[,2]*as.numeric(mydata1[,1]<tau)

        CW1=unoCW(time, status, Wi, Wi.star, emp1$rs, rs1.star, xi, emp1$Dhat)
        CW0=unoCW(time, status, Wi, Wi.star, emp0$rs, rs0.star, xi, emp0$Dhat)

        #----- result  ---#
	    temp1[k] <-  (CW1$Wa + CW1$Wb + CW1$Wg)*sqrt(n)
	    temp0[k] <-  (CW0$Wa + CW0$Wb + CW0$Wg)*sqrt(n)
	    temp[k] <-  temp1[k]-temp0[k]

 	}
	se=sqrt(var(temp)/n)
	Delta.hat<-emp1$Dhat - emp0$Dhat
	low95<-Delta.hat - 1.959*se
	upp95<-Delta.hat + 1.959*se

	se1=sqrt(var(temp1)/n)
	c1low95<-emp1$Dhat - 1.959*se1
	c1upp95<-emp1$Dhat + 1.959*se1

	se0=sqrt(var(temp0)/n)
	c0low95<-emp0$Dhat - 1.959*se0
	c0upp95<-emp0$Dhat + 1.959*se0

    C0=c(emp0$Dhat,se0, c0low95, c0upp95)
    C1=c(emp1$Dhat,se1, c1low95, c1upp95)
    DeltaC=c(Delta.hat, se, low95, upp95)
    out=as.matrix(rbind(C1,C0,DeltaC))
    rownames(out)=c("Model1","Model0","Delta")
    colnames(out)=c("Est","SE","Lower95","Upper95")
    return(out)   
#	return(list(C0=C0, C1=C1, DeltaC=DeltaC))
}



#######################################################
## Function: Est.PH: Fit Cox and return beta'x
## --fixed in version 2
#######################################################
Est.PH<-function(mydata){

	covs<-as.matrix(mydata[,-1:-2])
	ft<-coxph(Surv(mydata[,1],mydata[,2]) ~ covs)
	rs<-as.matrix(covs)%*%as.vector(ft$coefficient)
	
	#---- score ----#
	time <- mydata[,1]
	n<-nrow(covs); p<-ncol(covs)
	S0<-rep(0, n) ; S1<-matrix(0, n, p)
	for (i in 1:n){
		S0[i]=as.numeric(time>=time[i])%*%exp(rs)/n ;
		S1[i,]=apply(VEC2MAT(as.numeric(time>=time[i]),p)*VEC2MAT(exp(rs),p)*t(as.matrix(covs)), 1, mean) ;
		}
 	 Ui<-(covs - S1/S0)*mydata[,2]
	
	return(list(beta=ft$coefficient, var=ft$var, rs=rs, Ui=Ui, ft=ft))
	}


#============================================================
VEC2MAT<-function(vc, dm){
     matrix(vc, ncol=length(vc), nrow=dm, byrow=TRUE)
    }


#######################################################
## Function: kmcens: gives Keplan-Meier of Censoring ## 
#######################################################
## I: time      (nx1)
##    status    (nx1)  1=fail, 0=censor  
##    tau       (1x1)  truncate point P(C>tau)>0
##=====================================================
## O: surv	(tx1) G(t) <---------Kaplan-Meier
##    nelson	(tx1) Lambda_C(t) <--Nelson-Alan 
##    distinct	(tx1) t <------------Unique time point at observed event
##    pi_0(t)   (tx1) pr[1{X>=t}] 
##    pi_X(t)   (tx1) pi_X = pr[ 1{X>=t} * G(t) ]
##    pi_T(t)   (tx1) pi_T = pr[ 1{X>=t} / G(t) ]
##    Fn(t)     (tx1)   Fn = pr[ 1{X<=t} & stuatus==1]
##    Mi(t)     (nxt) Martingale for censoring
##    psii(t)	(nxt) iid term
#######################################################
kmcens<-function(time, status, tau){

#-- initial --#

 distinct<-unique(sort(time))
 #--- do not do this ---#
 # distinct<-unique(sort(time[status==1 & time < tau]))

 t<-length(distinct)
 n<-length(time)
 surv<-rep(0,t)
 nel.wk<-rep(0,t)
 nelson<-rep(0,t)

#- i = 1 -#
   yi<-sum(as.numeric(time>=distinct[1]))
   di<-sum(as.numeric(time==distinct[1] & status==0))
   surv[1]<- 1 * (1-di/yi)
   nel.wk[1]<-di/yi

#- i = 2 to t -#
 for (i in 2:t){
   yi<-sum(as.numeric(time>=distinct[i]))
   di<-sum(as.numeric(time==distinct[i] & status==0))
   surv[i]<-surv[i-1]*(1-di/yi)
   nel.wk[i]<-di/yi
   }

 #- shift -#
  surv[2:t]<-surv[1:(t-1)]
  surv[1]<-1
  nel.wk[2:t]<-nel.wk[1:(t-1)]
  nel.wk[1]<-0
 #---------#

  nelson<-cumsum(nel.wk)

 #--------------------------------------------------------#
 # pi_0(t) --- (tx1 vector)--- pr[ 1{X>=t}]              -#
 # pi_X(t) --- (tx1 vector)--- pr[ 1{X>=t} * G(t) ]      -#
 # pi_T(t) --- (tx1 vector)--- pr[ 1{X>=t} / G(t) ]      -#
 #   Fn(t) --- (tx1 vector)--- pr[ 1{X<=t} & stuatus==1] -#
 #--------------------------------------------------------#
 pi_0<-rep(0,t)
 pi_X<-rep(0,t)
 pi_T<-rep(0,t)
   Fn<-rep(0,t)
 for (i in 1:t){
  yi<- as.numeric(time >= distinct[i])
  ni<- as.numeric(time <= distinct[i] & status==1)
  pi_0[i]<-mean(yi)
  pi_X[i]<-mean(yi)*surv[i]
  pi_T[i]<-mean(yi)/surv[i]
    Fn[i]<-mean(ni)
 }

 #-------------------------------#
 #  Mi(t) --- (n x t) Martingale #
 #-------------------------------#
 Mi<-matrix(0, n, t)
 wk1<-matrix(0, n, t)
 wk2<-matrix(0, n, t)

 for (j in 1:t){ 
   wk1[,j]<- as.numeric(time <= distinct[j] & status==0)
   wk2[,j]<- as.numeric(time >= distinct[j])
 }
 for (k in 1:n){
   Mi[k,] <- wk1[k,] - wk2[k,]*nelson
 }
 

 #--------------------------------#
 # psii(t)	   (n x t matrix) #
 #--------------------------------#
 psii<-matrix(0, n, t) 
  wk1<-matrix(0, n, t)
  wk2<-matrix(0, n, t)

 for (j in 1:t){ 
   wk1[,j]<- as.numeric(time == distinct[j] & status==0)
   wk2[,j]<- as.numeric(time >= distinct[j])
 }

for (i in 1:n){ 
   psii[i,]<- (cumsum(wk1[i,]/pi_0) - cumsum(nel.wk*wk2[i,]/pi_0))
 }


 #----------------------------------------------------------#
 wk1<-NULL ; wk2<-NULL
 return(list(surv=surv, nelson=nelson, distinct=distinct, pi_0=pi_0, pi_X=pi_X, pi_T=pi_T, Mi=Mi, psii=psii, Fn=Fn))
}


#######################################################
# Function: cvC: return cross-validation estimates 
#######################################################
cvC=function(mydata, tau, cvK=10, Rep=10){

time=mydata[,1]
fail=mydata[,2]
X=as.matrix(mydata[,c(-1,-2)])

Dhat=rep(0,cvK)
Cest=rep(0, Rep)	

for(r in 1:Rep){
    group=sample(1:cvK,nrow(mydata),replace=TRUE)

    for(kk in 1:cvK){
     mod0=coxph(Surv(time,fail)~X,subset=group!=kk)
     XTest=as.matrix(X[group==kk,])
     scores=XTest%*%mod0$coef
     D=data.frame(cbind(time[group==kk],fail[group==kk],scores))
     Dhat[kk]=Est.Cval(D,tau)$Dhat
    }
    Cest[r]=weighted.mean(Dhat,table(group))
 }
  return(mean(Cest))
}

Try the survC1 package in your browser

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

survC1 documentation built on Feb. 10, 2021, 5:06 p.m.