# Suite of programs for case-cohort analysis

# Main program

cch_function(formula, data=sys.parent(), subcoh, id, cohort.size, 
         method="Prentice"){
call <- match.call()
# Check id, subcoh and cohort.size variables
if(length(id)!=length(unique(id))) stop("Multiple records per id not allowed")
tt_table(subcoh)
if(min(charmatch(names(tt),c("0","1"),0))==0)
   stop("Permissible values for sucohort indicator are 0,1")
if(length(id)>cohort.size) stop("Number of records greater than cohort size")
nn_cohort.size	
# Evaluate model formula
m_match.call(expand.dots=F)
m$method_m$cohort.size_m$id_m$subcoh_NULL
m[[1]]_as.name("model.frame")
m_eval(m,sys.parent())
Terms_attr(m,"terms")
Y_model.extract(m, response)
if(!inherits(Y, "Surv")) stop("Response must be a survival object")
type <- attr(Y, "type")
itype<-charmatch(type,c("right","counting"),nomatch=0)
cens<-switch(itype+1,stop(paste("Cox model doesn't support \"", type, "\" 
      survival data", sep = "")),Y[,2],Y[,3])
cc<-cens+1-subcoh
texit<-switch(itype+1, stop(), Y[,1], Y[,2])
tenter<-switch(itype+1, stop(), rep(0,length(texit)), Y[,1])
X_model.matrix(Terms, m)
X_X[,2:ncol(X)]
imeth_charmatch(method,c("Prentice","SelfPren","LinYing"), nomatch=0) 
method.name_switch(imeth+1, stop("Method doesn't exist"), "Prentice", 
"SelfPren", "LinYing")
z_call(method.name, tenter=tenter, texit=texit, cc=cc, id=id, X=X, ntot=nn)
out_eval(z, local=sys.parent())
out$method_method
names(out$coef) <- dimnames(X)[[2]]
        if(!is.null(out$var))
           dimnames(out$var) <- list(dimnames(X)[[2]], dimnames(X)[[2]])
        if(!is.null(out$naive.var))
           dimnames(out$naive.var) <- list(dimnames(X)[[2]], dimnames(X)[[2]])
out$call_call
out$cohort.size_nn
out$subcohort.size_tt[2]
class(out)_"cch"
out
}

# Subprograms

Prentice_function(tenter, texit, cc,  id, X, ntot){
eps_0.00000001
cens_as.numeric(cc>0) # Censorship indicators
subcoh_as.numeric(cc<2) # Subcohort indicators

# Calculate Prentice estimate
ent2_tenter
ent2[cc==2]_texit[cc==2]-eps
fit1_coxph(Surv(ent2,texit,cens)~X,eps=eps,x=T)

# Calculate Prentice estimate and variance
nd_sum(cens) # Number of failures
nc_sum(subcoh) # Number in subcohort
ncd_sum(as.numeric(cc==1)) #Number of failures in subcohort
X_as.matrix(X)
aent_c(tenter[cc>0],tenter[cc<2])
aexit_c(texit[cc>0],texit[cc<2])
aX_rbind(as.matrix(X[cc>0,]),as.matrix(X[cc<2,]))
aid_c(id[cc>0],id[cc<2])
dum_rep(-100,nd)
dum_c(dum,rep(0,nc))
gp_rep(1,nd)
gp_c(gp,rep(0,nc))
fit_coxph(Surv(aent,aexit,gp)~aX+offset(dum)+cluster(aid),eps=eps,x=T,
iter.max=35,init=fit1$coefficients)
db_resid(fit,type="dfbeta")
db_as.matrix(db)
db_db[gp==0,]
fit$naive.var_fit$naive.var+(1-(nc/ntot))*t(db)%*%(db)
fit$coef_fit1$coefficients
fit
}

SelfPren_function(tenter, texit, cc,  id, X, ntot){
eps_0.00000001
cens_as.numeric(cc>0) # Censorship indicators
subcoh_as.numeric(cc<2) # Subcohort indicators

# Calculate Self-Prentice estimate and variance
nd_sum(cens) # Number of failures
nc_sum(subcoh) # Number in subcohort
ncd_sum(as.numeric(cc==1)) #Number of failures in subcohort
X_as.matrix(X)
aent_c(tenter[cc>0],tenter[cc<2])
aexit_c(texit[cc>0],texit[cc<2])
aX_rbind(as.matrix(X[cc>0,]),as.matrix(X[cc<2,]))
aid_c(id[cc>0],id[cc<2])
dum_rep(-100,nd)
dum_c(dum,rep(0,nc))
gp_rep(1,nd)
gp_c(gp,rep(0,nc))
fit_coxph(Surv(aent,aexit,gp)~aX+offset(dum)+cluster(aid),eps=eps,x=T)
db_resid(fit,type="dfbeta")
db_as.matrix(db)
db_db[gp==0,]
fit$naive.var_fit$naive.var+(1-(nc/ntot))*t(db)%*%(db)
fit
}

LinYing_function(tenter, texit, cc,  id, X, ntot){
eps_0.00000001
cens_as.numeric(cc>0) # Censorship indicators
subcoh_as.numeric(cc<2) # Subcohort indicators
nd_sum(cens) # Number of failures
nc_sum(subcoh) # Number in subcohort
ncd_sum(as.numeric(cc==1)) #Number of failures in subcohort

# Calculate Lin-Ying estimate and variance
offs_rep((ntot-nd)/(nc-ncd),length(texit))
offs[cc>0]_1
loffs_log(offs)
fit_coxph(Surv(tenter, texit, cens)~X+offset(loffs)+cluster(id),
eps=eps,x=T)
db_resid(fit,type="dfbeta")
db_as.matrix(db)
db_db[cens==0,,drop=F]
dbm_apply(db,2,mean)
db_sweep(db,2,dbm)
fit$naive.var_fit$naive.var+(1-(nc-ncd)/(ntot-nd))*t(db)%*%db
fit
}

"summary.cch"<-
function(object)
{
# produces summary from an object of the class "cch"
      call<-object$call
	coef <- object$coef
	method <- object$method
	se <- sqrt(diag(object$naive.var))
      Z<- coef/se
      p<- 2*(1-pnorm(abs(Z)))
      cohort.size<-object$cohort.size
      subcohort.size<-object$subcohort.size
coefficients <- matrix(0, nrow = length(coef), ncol = 4)
dimnames(coefficients) <- list(names(coef), c("Value", 
			"SE", "Z", "p"))
coefficients[, 1] <- coef
coefficients[, 2] <- se
coefficients[, 3] <- Z
coefficients[, 4] <- 2*(1-p)
structure(list(call=call, method=method, cohort.size=cohort.size,
subcohort.size=subcohort.size, coefficients = coefficients), 
class = "summary.cch")
}