```# sigmaone function
sigmaone <- function(fit)
{
N <- fit\$byprod\$N
score_process <- fit\$byprod\$score_process
#sum of squares of score_process
sigma_one <- t(score_process) %*% score_process
sigma_one <- sigma_one / N
fit\$byprod\$sigma_one <- sigma_one
sigma_one
}

# psihat function
psihat <- function(fit)
{
comp <- fit\$byprod\$comp
binary <- fit\$byprod\$binary
N <- fit\$byprod\$N
X <- fit\$byprod\$X
firstfit <- fit\$byprod\$firstfit
Z <- fit\$byprod\$Z
survtime <- fit\$byprod\$survtime
Z_int <- fit\$byprod\$Z_int
coef_est <- fit\$coef
#there are differences whether using binary, the derivative part
if(binary) sum_part <- (X * (exp(-fitted(firstfit)) / (1 + exp(-fitted(firstfit))) ** 2))
else sum_part <- X
#a big difference between whether competing or not
if(!comp) integral_part <- Z * survtime - Z_int
else{
cause <- fit\$byprod\$cause
censorsurv <- fit\$byprod\$censorsurv
G_int <- fit\$byprod\$G_int
GZ_int <- fit\$byprod\$GZ_int
integral_part <- Z * survtime - Z_int + (cause >= 2) / censorsurv * (Z * G_int - GZ_int)
}
#
psi_hat <- matrix(rowSums(sapply(1:N, function(i) outer(integral_part[i, ], sum_part[i, ]))), ncol = ncol(sum_part)) * coef_est[2]
psi_hat <- psi_hat / N
psi_hat
}

# sigmatwo function
sigmatwo <- function(fit)
{
N <- fit\$byprod\$N
psi_hat <- psihat(fit)
firstfit <- fit\$byprod\$firstfit
sigma_two <- N * psi_hat %*% vcov(firstfit) %*% t(psi_hat)
sigma_two
}

# pihatt function
pihatt <- function(fit)
{
survtime <- fit\$byprod\$survtime
pi_hatt <- seq(length(survtime), 1)
pi_hatt
}

# qhatt function
qhatt <- function(fit)
{
comp <- fit\$byprod\$comp
if(!comp) q_hatt <- 0
else{
#the following functions are a decompostion of q_hatt
#the main motivation of the following functions can be checked in the help file
{
}

{
}

firstpart <- function(s_zero, Z, cause, censorsurv)
{
sum_part <- apply(Z * (cause >= 2) / censorsurv, 2, cumsum)
sum_part <- rbind(rep(0, ncol(sum_part)), sum_part[1:(nrow(sum_part) - 1), ])
integral_part <- rev(cumsum(rev(censorsurv * (cause == 1) / s_zero)))
first_part <- sum_part * integral_part
first_part
}

secondpart <- function(s_zero, Z_bar, cause, censorsurv)
{
sum_part <- cumsum((cause >= 2) / censorsurv)
sum_part <- c(0, sum_part[1:(length(sum_part) - 1)])
integral_part <- apply(apply(apply(Z_bar * censorsurv * (cause == 1) / s_zero, 2, rev), 2, cumsum), 2, rev)
second_part <- sum_part * integral_part
second_part
}

thirdpart <- function(survtime, coef_est, Z, cause, censorsurv)
{
sum_part <- apply(Z * as.numeric(Z %*% coef_est) * (cause >= 2) / censorsurv, 2, cumsum)
sum_part <- rbind(rep(0, ncol(sum_part)), sum_part[1:(nrow(sum_part) - 1), ])
integral_part <- Gint(survtime, censorsurv)
third_part <- sum_part * integral_part
third_part
}

fourthpart <- function(survtime, coef_est, Z, Z_bar, cause, censorsurv)
{
sum_part <- apply(Z * (cause >= 2) / censorsurv, 2, cumsum)
sum_part <- rbind(rep(0, ncol(sum_part)), sum_part[1:(nrow(sum_part) - 1), ])
timediff <- c(survtime[1], diff(survtime))
integral_part <- cumsum(rev(censorsurv * as.numeric(Z_bar %*% coef_est) * timediff))
integral_part <- rev(c(0, integral_part[1:(length(integral_part) - 1)]))
fourth_part <- sum_part * integral_part
fourth_part
}

fifthpart <- function(survtime, coef_est, Z_bar, cause, censorsurv)
{
sum_part <- cumsum((cause >= 2) / censorsurv)
sum_part <- c(0, sum_part[1:(length(sum_part) - 1)])
timediff <- c(survtime[1], diff(survtime))
integral_part <- apply(apply(censorsurv * Z_bar * as.numeric(Z_bar %*% coef_est) * timediff, 2, rev), 2, cumsum)
integral_part <- rbind(rep(0, ncol(integral_part)), integral_part[1:(nrow(integral_part) - 1), ])
integral_part <- apply(integral_part, 2, rev)
fifth_part <- sum_part * integral_part
fifth_part
}
cause <- fit\$byprod\$cause
Z <- fit\$byprod\$Z
Z_bar <- fit\$byprod\$Z_bar
s_zero <- fit\$byprod\$s_zero
censorsurv <- fit\$byprod\$censorsurv
survtime <- fit\$byprod\$survtime
coef_est <- fit\$coef
first_part <- firstpart(s_zero, Z, cause, censorsurv)
second_part <- secondpart(s_zero, Z_bar, cause, censorsurv)
third_part <- thirdpart(survtime, coef_est, Z, cause, censorsurv)
fourth_part <- fourthpart(survtime, coef_est, Z, Z_bar, cause, censorsurv)
fifth_part <- fifthpart(survtime, coef_est, Z_bar, cause, censorsurv)
q_hatt <- leadingfirst_part - leadingsecond_part - (first_part - second_part + third_part - 2 * fourth_part + fifth_part)
}
q_hatt
}

# sigmathree function
sigmathree <- function(fit)
{
N <- fit\$byprod\$N
pZ <- fit\$byprod\$pZ
cause <- fit\$byprod\$cause
sigma_three <- qhatt(fit) / pihatt(fit)
sigma_three <- matrix(rowSums(apply(sigma_three, 1, function(x) outer(x, x)) * (cause == 0)), ncol = pZ)
sigma_three <- sigma_three / N
sigma_three
}

# betavarest function
betavarest <- function(fit)
{
N <- fit\$byprod\$N
useIV <- fit\$byprod\$useIV
comp <- fit\$byprod\$comp
omega_inv <- fit\$byprod\$omega_inv
if(!useIV && !comp){
sigma_one <- sigmaone(fit)
sigma_two <- 0
sigma_three <- 0
}
else if(useIV && !comp){
sigma_one <- sigmaone(fit)
sigma_two <- sigmatwo(fit)
sigma_three <- 0
}
else if(!useIV && comp){
sigma_one <- sigmaone(fit)
sigma_two <- 0
sigma_three <- sigmathree(fit)
}
else{
sigma_one <- sigmaone(fit)
sigma_two <- sigmatwo(fit)
sigma_three <- sigmathree(fit)
}
pararesult <- list(sigma_one = sigma_one,
sigma_two = sigma_two,
sigma_three = sigma_three)
fit\$byprod <- append(fit\$byprod, pararesult)
fit\$vcov <- omega_inv %*% (sigma_one + sigma_two + sigma_three) %*% omega_inv / N
fit
}

{
N <- fit\$byprod\$N
survtime <- fit\$byprod\$survtime
cause <- fit\$byprod\$cause
s_zero <- fit\$byprod\$s_zero
lead_part <- N * cumsum((cause == 1) / s_zero ^ 2)
}

# Dhatt function
Dhatt <- function(fit)
{
Z <- fit\$byprod\$Z
Z_bar <- fit\$byprod\$Z_bar
cause <- fit\$byprod\$cause
s_zero <- fit\$byprod\$s_zero
D_hatt <- apply((Z - Z_bar) * (cause == 1) / s_zero, 2, cumsum)
D_hatt
}

# YGint function
YGint <- function(fit)
{
#compute the integration of G(t) over s_zero
comp <- fit\$byprod\$comp
if(comp){
survtime <- fit\$byprod\$survtime
censorsurv <- fit\$byprod\$censorsurv
s_zero <- fit\$byprod\$s_zero
timediff <- c(survtime[1], diff(survtime))
temp <- cumsum(rev(timediff * censorsurv / s_zero))
temp <- c(0, temp[1:(length(temp) - 1)])
YG_int <- rev(temp)

}
else YG_int <- 0
YG_int
}

# szeroint function
szeroint <- function(fit)
{
survtime <- fit\$byprod\$survtime
s_zero <- fit\$byprod\$s_zero
timediff <- c(survtime[1], diff(survtime))
szero_int <- cumsum(timediff / s_zero)
szero_int
}

# Yint function
Yint <- function(fit)
{
#compute the integration of reciprocal of s_zero
szero_int <- szeroint(fit)
YG_int <- YGint(fit)
Y_int <- szeroint + YG_int
Y_int
}

# Ehatt function
Ehatt <- function(fit, i)
{
coef_est <- fit\$coef
firstfit <- fit\$byprod\$firstfit
X <- fit\$byprod\$X
N <- fit\$byprod\$N
binary <- fit\$byprod\$binary
censorsurv <- fit\$byprod\$censorsurv
cause <- fit\$byprod\$cause

if(binary) A <- (X * (exp(-fitted(firstfit)) / (1 + exp(-fitted(firstfit))) ** 2))
else A <- X
#the following functions are mainly motivated by the help file
#you may contact the author for the help file
szero_int <- szeroint(fit)
YG_int <- YGint(fit)
firstpart <- function(fit, i){
if(is.null(fit\$byprod\$Ehattfirst)){
sum_part <- coef_est[2] * A
integral_part <- szero_int
integral_part_diff <- c(integral_part[1], diff(integral_part))
first_part <- colSums(sum_part * integral_part)
fit\$byprod\$Ehattfirst <- list(sum_part = sum_part,
integral_part_diff = integral_part_diff,
first_part = first_part)
}
else{
sum_part <- fit\$byprod\$Ehattfirst\$sum_part
temp <- c(rep(0, i), rep(1, N - i))
sum_part <- sum_part * temp
integral_part_diff <- fit\$byprod\$Ehattfirst\$integral_part_diff
first_part <- fit\$byprod\$Ehattfirst\$first_part
first_part <- first_part - colSums(sum_part * integral_part_diff[i + 1])
fit\$byprod\$Ehattfirst\$first_part <- first_part
}
fit
}
secondpart <- function(fit, i){
comp <- fit\$byprod\$comp
if(comp){
if(is.null(fit\$byprod\$Ehattsecond)){
sum_part <- coef_est[2] * A * (cause >= 2) / censorsurv
integral_part <- YG_int
integral_part_diff <- c(rev(diff(rev(integral_part))), integral_part[N])
second_part <- colSums(sum_part * integral_part)
fit\$byprod\$Ehattsecond <- list(sum_part = sum_part,
integral_part_diff = integral_part_diff,
second_part = second_part)
}
else{
fit\$byprod\$Ehattsecond\$sum_part[i:N] <- 0
sum_part <- fit\$byprod\$Ehattsecond\$sum_part
integral_part_diff <- fit\$byprod\$Ehattsecond\$integral_part_diff
second_part <- fit\$byprod\$Ehattsecond\$second_part
second_part <- second_part - colSums(sum_part * integral_part_diff[i])
fit\$byprod\$Ehattsecond\$second_part <- second_part
}
}
else fit\$byprod\$Ehattsecond\$second_part <- 0
fit
}
fit <- firstpart(fit, i)
fit <- secondpart(fit, i)
first_part <- fit\$byprod\$Ehattfirst\$first_part
second_part <- fit\$byprod\$Ehattsecond\$second_part
fit\$byprod\$E_hatt <- first_part + second_part
fit
}

# EThetaEpart function
EThetaEpart <- function(fit)
{
N <- fit\$byprod\$N
useIV <- fit\$byprod\$useIV
comp <- fit\$byprod\$comp
if(!useIV) EThetaE_part <- 0
else{
if(is.null(fit\$byprod\$EThetaE_part)){
firstfit <- fit\$byprod\$firstfit
firstvcov <- vcov(firstfit)
EThetaE_part <- c()
for(i in N:1){
if(!is.null(fit\$byprod\$Ehattfirst)){
first <- !fit\$byprod\$Ehattfirst\$integral_part_diff[i + 1]
if(comp) second <- !fit\$byprod\$Ehattsecond\$integral_part_diff[i]
else second <- TRUE
if(first & second){
EThetaE_part <- c(EThetaE_part, temp)
next
}
}
fit <- Ehatt(fit, i)
E_hatt <- fit\$byprod\$E_hatt
temp <- t(E_hatt) %*% firstvcov %*% E_hatt
EThetaE_part <- c(EThetaE_part, temp)
}
EThetaE_part <- rev(EThetaE_part)
fit\$byprod\$EThetaE_part <- EThetaE_part
}
#else EThetaE_part <- fit\$byprod\$EThetaE_part
}
#fit <<- fit
#EThetaE_part
fit
}

# Ghatt function
Ghatt <- function(fit, newobsz)
{
survtime <- fit\$byprod\$survtime
Z_int <- fit\$byprod\$Z_int
G_hatt <- survtime %*% t(newobsz) - Z_int
}

# qprimehatt function
qprimehatt <- function(fit, i)
{
#this is the q_t(u) function in the paper for estimating the variance of hazard function
#i is the index of event time
#the following functions are mainly motivated by the help file
#you may contact the author for the help file
#the first part is automatically zero
secondpart <- function(fit, i)
{
N <- fit\$byprod\$N
cause <- fit\$byprod\$cause
s_zero <- fit\$byprod\$s_zero
censorsurv <- fit\$byprod\$censorsurv
#return a vector of v, the time t is fixed and decreasing each time it is called
if(is.null(fit\$byprod\$qprimehattsecond)){
sum_part <- cumsum((cause >= 2) / censorsurv)
sum_part <- c(0, sum_part[1:(length(sum_part) - 1)])
integral_part <- rev(cumsum(rev(censorsurv * (cause == 1) / s_zero ^ 2)))
integral_part_diff <- censorsurv * (cause == 1) / s_zero ^ 2
second_part <- sum_part * integral_part
fit\$byprod\$qprimehattsecond <- list(sum_part = sum_part,
integral_part_diff = integral_part_diff,
second_part = second_part)
}
else{
if(i < N){
if(i < N - 1){
fit\$byprod\$qprimehattsecond\$sum_part[(i + 2):N] <- 0
}
sum_part <- fit\$byprod\$qprimehattsecond\$sum_part
integral_part_diff <- fit\$byprod\$qprimehattsecond\$integral_part_diff
second_part <- fit\$byprod\$qprimehattsecond\$second_part
second_part <- second_part - sum_part * integral_part_diff[i + 1]
fit\$byprod\$qprimehattsecond\$sum_part[i + 1] <- 0
fit\$byprod\$qprimehattsecond\$second_part <- second_part
}
}
#fit <<- fit
fit
}
thirdpart <- function(fit, i)
{
N <- fit\$byprod\$N
coef_est <- fit\$coef
cause <- fit\$byprod\$cause
Z <- fit\$byprod\$Z
YG_int <- YGint(fit)
censorsurv <- fit\$byprod\$censorsurv
if(is.null(fit\$byprod\$qprimehattthird)){
sum_part <- cumsum(Z %*% coef_est * (cause == 2) / censorsurv)
sum_part <- c(0, sum_part[1:(length(sum_part) - 1)])
integral_part <- YG_int
integral_part_diff <- c(rev(diff(rev(integral_part))), YG_int[N])
third_part <- sum_part * integral_part
fit\$byprod\$qprimehattthird <- list(sum_part = sum_part,
integral_part_diff = integral_part_diff,
third_part = third_part)
}
else{
if(i < N){
fit\$byprod\$qprimehattthird\$sum_part[(i + 1):N] <- 0
sum_part <- fit\$byprod\$qprimehattthird\$sum_part
integral_part_diff <- fit\$byprod\$qprimehattthird\$integral_part_diff
third_part <- fit\$byprod\$qprimehattthird\$third_part
third_part <- third_part - sum_part * integral_part_diff[i + 1]
fit\$byprod\$qprimehattthird\$third_part <- third_part
}
}
#fit <<- fit
fit
}
YGZbetaint <- function(fit)
{
survtime <- fit\$byprod\$survtime
censorsurv <- fit\$byprod\$censorsurv
s_zero <- fit\$byprod\$s_zero
timediff <- c(survtime[1], diff(survtime))
coef_est <- fit\$coef
Z_bar <- fit\$byprod\$Z_bar
temp <- cumsum(rev(timediff * Z_bar %*% coef_est * censorsurv / s_zero))
temp <- c(0, temp[1:(length(temp) - 1)])
YGZbeta_int <- rev(temp)
YGZbeta_int
}

fourthpart <- function(fit, i)
{
N <- fit\$byprod\$N
cause <- fit\$byprod\$cause
censorsurv <- fit\$byprod\$censorsurv
if(is.null(fit\$byprod\$qprimehattfourth)){
sum_part <- cumsum((cause >= 2) / censorsurv)
sum_part <- c(0, sum_part[1:(length(sum_part) - 1)])
integral_part <- YGZbetaint(fit)
integral_part_diff <- c(rev(diff(rev(integral_part))), integral_part[N])
fourth_part <- sum_part * integral_part
fit\$byprod\$qprimehattfourth <- list(sum_part = sum_part,
integral_part_diff = integral_part_diff,
fourth_part = fourth_part)
}
else{
if(i < N){
fit\$byprod\$qprimehattfourth\$sum_part[(i + 1):N] <- 0
sum_part <- fit\$byprod\$qprimehattfourth\$sum_part
integral_part_diff <- fit\$byprod\$qprimehattfourth\$integral_part_diff
fourth_part <- fit\$byprod\$qprimehattfourth\$fourth_part
fourth_part <- fourth_part - sum_part * integral_part_diff[i + 1]
fit\$byprod\$qprimehattfourth\$fourth_part <- fourth_part
}
}
#fit <<- fit
fit
}
secondpart(fit, i)
thirdpart(fit, i)
fourthpart(fit, i)
second_part <- fit\$byprod\$qprimehattsecond\$second_part
third_part <- fit\$byprod\$qprimehattthird\$third_part
fourth_part <- fit\$byprod\$qprimehattfourth\$fourth_part
qprime_hatt <- 0 - second_part - third_part + fourth_part
fit\$byprod\$qprime_hatt <- qprime_hatt
#fit <<- fit
fit
}

# qprimepihatt function
qprimepihatt <- function(fit)
{
#compute the q_t(u) part in the paper
if(is.null(fit\$byprod\$qprimepi_hatt)){
comp = fit\$byprod\$comp
if(comp){
N <- fit\$byprod\$N
cause <- fit\$byprod\$cause
pi_hatt <- pihatt(fit)
qprimepi_hatt <- c()
for(i in N:1){
if(!is.null(fit\$byprod\$qprimehattsecond)){
second <- !fit\$byprod\$qprimehattsecond\$integral_part_diff[i + 1]
third <- !fit\$byprod\$qprimehattthird\$integral_part_diff[i + 1]
fourth <- !fit\$byprod\$qprimehattfourth\$integral_part_diff[i + 1]
if(second&third&fourth) {
qprimepi_hatt <- c(qprimepi_hatt, temp)
next
}
}
fit <- qprimehatt(fit, i)
temp <- fit\$byprod\$qprime_hatt / pi_hatt
temp <- temp ^ 2
temp <- sum(temp * (cause == 0))
qprimepi_hatt <- c(qprimepi_hatt, temp)
}
qprimepi_hatt <- rev(qprimepi_hatt)
qprimepi_hatt <- N * qprimepi_hatt
}else qprimepi_hatt <- 0
fit\$byprod\$qprimepi_hatt <- qprimepi_hatt
}#else qprimepi_hatt <- fit\$byprod\$qprimepi_hatt
#fit <<- fit
#qprimepi_hatt
fit
}

# hazardpredvarest function
hazardpredvarest <- function(newobsz, fit = NULL)
{
N <- fit\$byprod\$N
omega_inv <- fit\$byprod\$omega_inv
#qprimepi_hatt <- qprimepihatt(fit)
fit <- qprimepihatt(fit)
qprimepi_hatt <- fit\$byprod\$qprimepi_hatt
betavar_est <- fit\$vcov
D_hatt <- Dhatt(fit)
#EThetaE_part <- EThetaEpart(fit)
fit <- EThetaEpart(fit)
EThetaE_part <- fit\$byprod\$EThetaE_part
G_hatt <- Ghatt(fit, newobsz)
qprimepi_hatt = qprimepi_hatt,
D_hatt = D_hatt,
EThetaE_part = EThetaE_part,
G_hatt = G_hatt)
fit\$byprod <- append(fit\$byprod, hazard_aux)
#the variance estimate at each point for the hazard function
hazardpredvar_est <- (lead_part + qprimepi_hatt + N * rowSums(G_hatt %*% betavar_est * G_hatt)
+ EThetaE_part + 2 * rowSums(G_hatt %*% omega_inv * D_hatt)) / N
fit\$hazardpredvar_est <- hazardpredvar_est
fit
}
```