stcch_function(formula, data=sys.parent(), subcoh, id, stratum, stratum.sizes, 
               method="I"){
call <- match.call()
# Check id, subcoh and stratum 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")
jj_max(stratum)
if(length(unique(stratum))!=jj ) stop("Problem with stratum indicators")
nn_stratum.sizes  # Cohort stratum sizes
if(any(nn==0)) stop("Zero stratum size in cohort")
if(length(nn) != jj) Stop("length(stratum.sizes) != number of strata")
n_table(stratum)  # Sample stratum sizes
m_table(stratum[subcoh==1])
if(length(m)<jj) Stop("At least one stratum not represented in subcohort")
if(max(as.vector(n)-nn)>0) stop("Cohort stratum smaller than sample stratum")
# Evaluate model formula
m_match.call(expand.dots=F)
m$method_m$stratum.sizes_m$stratum_m$vars_m$id_m$cc_m$tenter_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("I","II"), nomatch=0) 
method.name_switch(imeth+1,stop("Method doesn't exist"), "BorganI", "BorganII")
z_call(method.name, ent=tenter, exit=texit, cc=cc, id=id, X=X, stratum=stratum, stratum.sizes=nn)
out_eval(z, local=sys.parent())
out$call_call
out$method_method
out$stratum.sizes_nn
out$subcohort.stratum.sizes_n
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]])
        if(!is.null(out$opt))
           dimnames(out$opt)_list(NULL,dimnames(X)[[2]])
class(out)_"stcch"
out
}

BorganI_function(ent, exit, cc,  id, X, stratum, stratum.sizes){
eps_0.00000001
nobs_length(exit)
idx_1:length(nobs)
jj_max(stratum)
nn_stratum.sizes  # Cohort stratum sizes
n_table(stratum)  # Sample stratum sizes
d_table(stratum[cc>0]) # Failures in each stratum
tt_table(cc,stratum)
cens_as.numeric(cc>0) # Failure 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
m0_tt[1,] # Subcohort stratum sizes (noncases only)
if (ncd>0) m_m0+tt[2,] else m_m0 #Subcohort stratum sizes
X <- as.matrix(X)
kk_ncol(X) # Number of variables
wt <- as.vector(nn/m) # Weights for Estimator I
stratum_c(stratum[cc>0],stratum[cc<2])
w_wt[stratum]
ent <- c(ent[cc > 0], ent[cc < 2])
exit <- c(exit[cc > 0], exit[cc < 2])
X <- rbind(as.matrix(X[cc > 0,  ]), 
	as.matrix(X[cc < 2,  ]))
id <- 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))
w[gp==1]_1
fit_coxph(Surv(ent,exit,gp)~X+offset(dum)+cluster(id),
weights=w, eps=eps,x=T, iter.max=25)  
score <- resid(fit, type = "score", weighted=F)
sc_resid(fit, type="score", collapse=id, weighted=T)
score <- as.matrix(score)
score <- score[gp == 0,,drop=F]
st_stratum[gp==0]
sto_st %o% rep(1,kk)
Index_col(score)
tscore_tapply(score,list(sto,Index),mean)
pscore_tapply(score,list(sto,Index))
score_score-tscore[pscore]
delta_matrix(0,kk,kk)
opt_NULL
for (j in 1:jj) {
       temp_t(score[st==j,])%*%score[st==j,]/(m[j]-1)
       delta_as.matrix(delta+(wt[j]-1)*nn[j]*temp,kk,kk)
  if(is.null(opt)) 
       opt_nn[j]*sqrt(diag(fit$naive.var %*% temp %*% fit$naive.var)) 
   else
       opt_rbind(opt,nn[j]*sqrt(diag(fit$naive.var %*% temp %*% fit$naive.var))) 
                }
z_apply(opt,2,sum)
fit$opt_sweep(opt,2,z,FUN="/")
fit$var_fit$naive.var+fit$naive.var%*%delta%*%fit$naive.var
fit$delta_delta
fit$sc_sc
fit
}

BorganII_function(ent, exit, cc,  id, X, stratum, stratum.sizes){
eps_0.00000001
jj_max(stratum)
nn_stratum.sizes  # Cohort stratum sizes
n_table(stratum)  # Sample stratum sizes
d_table(stratum[cc>0]) # Failures in each stratum
tt_table(cc,stratum)
cens_as.numeric(cc>0) # Failure 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
m0_tt[1,] # Subcohort stratum sizes (controls only)
if (ncd>0) m_m0+tt[2,] else m_m0 #Subcohort stratum sizes
X <- as.matrix(X)
kk_ncol(X) # Number of variables
nn0_nn-as.vector(d) #Noncases in cohort
wt <- as.vector(nn0/m0)
w_wt[stratum]
w[cens==1]_1
fit_coxph(Surv(ent,exit,cens)~X+cluster(id),
weights=w,eps=eps,x=T, iter.max=25)  # Borgan Estimate II
score <- resid(fit, type = "score", weighted=F)
sc_resid(fit,type="score", collapse=id, weighted=T)
score <- as.matrix(score)
score <- score[cens == 0,,drop=F] # Scores for controls
st_stratum[cens==0] # Stratum indicators for controls
sto_st %o% rep(1,kk)
Index_col(score)
tscore_tapply(score,list(sto,Index),mean) # Within stratum control score means
pscore_tapply(score,list(sto,Index))
score_score-tscore[pscore] # Subtract off within stratum score means
delta_matrix(0,kk,kk)
opt_NULL
for (j in 1:jj) {
       temp_t(score[st==j,])%*%score[st==j,]/(m0[j]-1) # Borgan equation (19)
       delta_delta+(wt[j]-1)*nn0[j]*temp # Borgan equation (17)
   if(is.null(opt)) 
       opt_nn0[j]*sqrt(diag(fit$naive.var %*% temp %*% fit$naive.var)) 
   else
       opt_rbind(opt,nn0[j]*sqrt(diag(fit$naive.var %*% temp %*% fit$naive.var))) 
}
z_apply(opt,2,sum)
fit$opt_sweep(opt,2,z,FUN="/")
fit$var_fit$naive.var+fit$naive.var %*% delta %*% fit$naive.var
fit$delta_delta
fit$sc_sc
fit
}

"summary.stcch"<-
function(object)
{
# produces summary from an object of the class "stcch"
      call<-object$call
	coef <- object$coef
	method <- paste("Borgan",object$method)
	se <- sqrt(diag(object$var))
      Z<- coef/se
      p<- pnorm(Z)
      cohort.stratum.sizes<-object$stratum.sizes
      subcohort.stratum.sizes<-object$subcohort.stratum.sizes
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.stratum.sizes=cohort.stratum.sizes,
subcohort.stratum.sizes=subcohort.stratum.sizes, coefficients = coefficients), 
class = "summary.stcch")
}