# R/funs-deriv.R In bcROCsurface: Bias-Corrected Methods for Estimating the ROC Surface of Continuous Diagnostic Tests

#### Defines functions est_func_spe_derivest_func_speest_func_ipw_derivest_func_ipwest_func_ie_derivest_func_iepi_est_funcmlog_est_funcpi_inv_derivrho_derivh_deriv

```####========================================================================####
## The R code for some specials functions, which are using in the computation ##
## for the asymptotic covariances of TCF's at a specified cut point.          ##
##                                                                            ##
####========================================================================####

## the derivative of function h(beta, theta, tau)
h_deriv <- function(thet, bet, n_nuis) {
res <- matrix(0, nrow = (6 + n_nuis), ncol = 3)
res[, 1] <- c(bet[1] / thet[1]^2, 0, -1 / thet[1], 0, 0, 0, rep(0, n_nuis))
res[, 2] <- c(0, -(bet[2] - bet[3]) / thet[2]^2, 0, 1 / thet[2],
-1 / thet[2], 0, rep(0, n_nuis))
res[, 3] <- c(bet[4] / (1 - thet[1] - thet[2])^2,
bet[4] / (1 - thet[1] - thet[2])^2,
0, 0, 0, 1 / (1 - thet[1] - thet[2]), rep(0, n_nuis))
return(t(res))
}

## the first derivatives of rho.
rho_deriv <- function(u_matrix, rho_hat, ref_level = c("1", "2")) {
# u_matrix: the design matrix is used to fit the disease model.
del_2 <- -u_matrix * rho_hat[, 1] * rho_hat[, 2]
ref_level <- match.arg(ref_level)
if (ref_level == "1") {
del_1 <- u_matrix * rho_hat[, 1] * (1 - rho_hat[, 1])
res <- cbind(del_1, del_2)
} else if (ref_level == "2") {
del_1 <- u_matrix * rho_hat[, 2] * (1 - rho_hat[, 2])
res <- cbind(del_2, del_1)
} else {
stop("Choose reference level in set {`1`, `2`}.")
}
return(res)
}

## the first derivatives of 1/pi.
pi_inv_deriv <- function(u_matrix, pi_hat, pi_cov,
model = c("logit", "probit", "threshold")) {
# u_matrix: the design matrix is used to fit the vefication model.
model <- match.arg(model)
res <- switch(EXPR = model,
logit = -u_matrix * (1 - pi_hat) / pi_hat,
probit = -u_matrix *
as.numeric(dnorm(u_matrix %*% pi_cov)) / pi_hat,
threshold = -u_matrix / pi_hat^2)
return(res)
}

## Estimating function for rho_k = Pr(D_k = 1 | T, A).
mlog_est_func <- function(dd, vv, u_matrix, rho_hat) {
temp <- u_matrix * vv
b1 <- (dd[, 1] - rho_hat[, 1]) * temp
b2 <- (dd[, 2] - rho_hat[, 2]) * temp
return(cbind(b1, b2))
}

## Estimating function for pi = Pr(V = 1| T, A).
pi_est_func <- function(vv, u_matrix, pi_hat, pi_cov,
model = c("logit", "probit", "threshold")) {
model <- match.arg(model)
res <- switch(EXPR = model,
logit = (vv - pi_hat) * u_matrix,
probit = u_matrix * as.numeric(dnorm(u_matrix %*% pi_cov)) *
(vv - pi_hat) / (pi_hat * (1 - pi_hat)),
threshold = (vv / pi_hat - 1) * u_matrix / (1 - pi_hat))
return(res)
}

## Estimating function of Imputation Estimator (FI and MSI).
est_func_ie <- function(dd, vv, tt, u_matrix, cc, cov_rho, rho_hat,
thet, bet, m) {
n_par <- 6 + length(cov_rho)
res <- matrix(0, nrow = nrow(u_matrix), ncol = n_par)
res[, 1] <- vv * (m * dd[, 1] - thet[1] + (1 - m) * rho_hat[, 1]) +
(1 - vv) * (rho_hat[, 1] - thet[1])
res[, 2] <- vv * (m * dd[, 2] - thet[2] + (1 - m) * rho_hat[, 2]) +
(1 - vv) * (rho_hat[, 2] - thet[2])
res[, 3] <- vv * ((tt >= cc[1]) * (m * dd[, 1] + (1 - m) * rho_hat[, 1]) -
bet[1]) + (1 - vv) * ((tt >= cc[1]) * rho_hat[, 1] -
bet[1])
res[, 4] <- vv * ((tt >= cc[1]) * (m * dd[, 2] + (1 - m) * rho_hat[, 2]) -
bet[2]) + (1 - vv) * ((tt >= cc[1]) * rho_hat[, 2] -
bet[2])
res[, 5] <- vv * ((tt >= cc[2]) * (m * dd[, 2] + (1 - m) * rho_hat[, 2]) -
bet[3]) + (1 - vv) * ((tt >= cc[2]) * rho_hat[, 2] -
bet[3])
res[, 6] <- vv * ((tt >= cc[2]) * (m * dd[, 3] + (1 - m) * rho_hat[, 3]) -
bet[4]) + (1 - vv) * ((tt >= cc[2]) * rho_hat[, 3] -
bet[4])
res[, 7:n_par] <- mlog_est_func(dd, vv, u_matrix, rho_hat)
return(res)
}

## Derivative of Estimating function of Imputation Estimator (FI and MSI).
est_func_ie_deriv <- function(vv, tt, u_matrix, cc, cov_rho, rho_hat, rho_hes,
thet, bet, m) {
der_rho_1 <- rho_deriv(u_matrix, rho_hat, ref_level = "1")
der_rho_2 <- rho_deriv(u_matrix, rho_hat, ref_level = "2")
der_rho_3 <- -(der_rho_1 + der_rho_2)
temp1 <- (1 - m * vv) * der_rho_1
temp2 <- (1 - m * vv) * der_rho_2
a <- rbind(colSums(temp1), colSums(temp2))
b1 <- colSums((tt >= cc[1]) * temp1)
b2 <- colSums((tt >= cc[1]) * temp2)
b3 <- colSums((tt >= cc[2]) * temp2)
b4 <- colSums((tt >= cc[2]) * (1 - m * vv) * der_rho_3)
b <- rbind(b1, b2, b3, b4)
res_left <- rbind(diag(-length(tt), 6, 6),
matrix(0, nrow = length(cov_rho), ncol = 6))
res_right <- rbind(a, b, -rho_hes)
res <- cbind(res_left, res_right)
return(res)
}

## Estimating function of IPW Estimator
est_func_ipw <- function(dd, vv, tt, u_matrix, cc, cov_pi, pi_hat, thet, bet,
model) {
n_par <- 6 + length(cov_pi)
res <- matrix(0, nrow = length(tt), ncol = n_par)
res[, 1] <- vv * (dd[, 1] - thet[1]) / pi_hat
res[, 2] <- vv * (dd[, 2] - thet[2]) / pi_hat
res[, 3] <- vv * ((tt >= cc[1]) * dd[, 1] - bet[1]) / pi_hat
res[, 4] <- vv * ((tt >= cc[1]) * dd[, 2] - bet[2]) / pi_hat
res[, 5] <- vv * ((tt >= cc[2]) * dd[, 2] - bet[3]) / pi_hat
res[, 6] <- vv * ((tt >= cc[2]) * dd[, 3] - bet[4]) / pi_hat
res[, 7:n_par] <- pi_est_func(vv, u_matrix, pi_hat, cov_pi, model)
return(res)
}

## Derivative of Estimating function of IPW Estimator
est_func_ipw_deriv <- function(dd, vv, tt, u_matrix, cc, cov_pi, pi_hat, pi_hes,
thet, bet, model) {
der_pi_inv <- pi_inv_deriv(u_matrix, pi_hat, cov_pi, model)
term <- vv * der_pi_inv
a1 <- colSums(term * (dd[, 1] - thet[1]))
a2 <- colSums(term * (dd[, 2] - thet[2]))
b11 <- colSums(term * ((tt >= cc[1]) * dd[, 1] - bet[1]))
b12 <- colSums(term * ((tt >= cc[1]) * dd[, 2] - bet[2]))
b22 <- colSums(term * ((tt >= cc[2]) * dd[, 2] - bet[3]))
b23 <- colSums(term * ((tt >= cc[2]) * dd[, 3] - bet[4]))
res_left <- rbind(diag(-sum(vv / pi_hat), 6, 6),
matrix(0, nrow = length(cov_pi), ncol = 6))
res_right <- rbind(a1, a2, b11, b12, b22, b23, -pi_hes)
res <- cbind(res_left, res_right)
return(res)
}

## Estimating function of SPE Estimator
est_func_spe <- function(dd, vv, tt, u_matrix_rho, u_matrix_pi, cc, cov_rho,
rho_hat, cov_pi, pi_hat, thet, bet, model) {
n_par <- 6 + length(cov_rho) + length(cov_pi)
res <- matrix(0, nrow = length(tt), ncol = n_par)
res[, 1] <- vv * (dd[, 1] - thet[1]) / pi_hat -
(rho_hat[, 1] - thet[1]) * (vv - pi_hat) / pi_hat
res[, 2] <- vv * (dd[, 2] - thet[2]) / pi_hat -
(rho_hat[, 2] - thet[2]) * (vv - pi_hat) / pi_hat
res[, 3] <- vv * ((tt >= cc[1]) * dd[, 1] - bet[1]) / pi_hat -
(vv - pi_hat) * ((tt >= cc[1]) * rho_hat[, 1] - bet[1]) / pi_hat
res[, 4] <- vv * ((tt >= cc[1]) * dd[, 2] - bet[2]) / pi_hat -
(vv - pi_hat) * ((tt >= cc[1]) * rho_hat[, 2] - bet[2]) / pi_hat
res[, 5] <- vv * ((tt >= cc[2]) * dd[, 2] - bet[3]) / pi_hat -
(vv - pi_hat) * ((tt >= cc[2]) * rho_hat[, 2] - bet[3]) / pi_hat
res[, 6] <- vv * ((tt >= cc[2]) * dd[, 3] - bet[4]) / pi_hat -
(vv - pi_hat) * ((tt >= cc[2]) * rho_hat[, 3] - bet[4]) / pi_hat
res[, 7:(6 + length(cov_rho))] <- mlog_est_func(dd, vv, u_matrix_rho, rho_hat)
res[, (7 + length(cov_rho)):n_par] <- pi_est_func(vv, u_matrix_pi, pi_hat,
cov_pi, model)
return(res)
}

## Derivative of Estimating function of SPE Estimator
est_func_spe_deriv <- function(dd, vv, tt, u_matrix_rho, u_matrix_pi, cc,
cov_rho, rho_hat, rho_hes, cov_pi, pi_hat,
pi_hes, thet, bet, model) {
der_rho_1 <- rho_deriv(u_matrix_rho, rho_hat, ref_level = "1")
der_rho_2 <- rho_deriv(u_matrix_rho, rho_hat, ref_level = "2")
der_rho_3 <- -(der_rho_1 + der_rho_2)
temp_v <- (1 - vv / pi_hat)
temp1 <- temp_v * der_rho_1
temp2 <- temp_v * der_rho_2
der_pi_inv <- pi_inv_deriv(u_matrix_pi, pi_hat, cov_pi, model)
temp3 <- vv * der_pi_inv
h <- rbind(colSums(temp1), colSums(temp2))
g1 <- colSums((tt >= cc[1]) * temp1)
g2 <- colSums((tt >= cc[1]) * temp2)
g3 <- colSums((tt >= cc[2]) * temp2)
g4 <- colSums((tt >= cc[2]) * temp_v * der_rho_3)
g <- rbind(g1, g2, g3, g4)
d1 <- colSums(temp3 * (rho_hat[, 1] - dd[, 1]))
d2 <- colSums(temp3 * (rho_hat[, 2] - dd[, 2]))
e11 <- colSums(temp3 * ((tt >= cc[1]) * (rho_hat[, 1] - dd[, 1])))
e12 <- colSums(temp3 * ((tt >= cc[1]) * (rho_hat[, 2] - dd[, 2])))
e22 <- colSums(temp3 * ((tt >= cc[2]) * (rho_hat[, 2] - dd[, 2])))
e23 <- colSums(temp3 * ((tt >= cc[2]) * (rho_hat[, 3] - dd[, 3])))
res_left <- rbind(diag(-length(tt), 6, 6),
matrix(0, nrow = length(cov_rho) + length(cov_pi),
ncol = 6))
res_right_rho <- rbind(h, g, -rho_hes,
diag(0, nrow = length(cov_pi), ncol = length(cov_rho)))
res_right_pi <- rbind(d1, d2, e11, e12, e22, e23,
diag(0, nrow = length(cov_rho), ncol = length(cov_pi)),
-pi_hes)
res <- cbind(res_left, res_right_rho, res_right_pi)
return(res)
}
```

## Try the bcROCsurface package in your browser

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

bcROCsurface documentation built on Sept. 9, 2023, 9:07 a.m.