context("Tau variance estimator")
Tau_est <- function(yA, yB) {
U <- sapply(yB, function(j) (j > yA) - (j < yA))
m <- length(yA)
n <- length(yB)
Tau <- sum(U) / (m * n)
# Hanley-McNeil variance estimator
Q1 <- sum(rowSums(U)^2) / (m * n^2)
Q2 <- sum(colSums(U)^2) / (m^2 * n)
V_HM <- (1 - Tau^2 + (n - 1) * (Q1 - Tau^2) + (m - 1) * (Q2 - Tau^2)) / (m * n)
# Unbiased variance estimator
X <- sum(U^2) / (m * n)
V_U <- (1 + n * Q1 + m * Q2 - (m + n - 1) * Tau^2 - 2 * X) / ((m - 1) * (n - 1))
# Jack-knife variance estimator
V_JK <- Q1 / (m - 1) + Q2 / (n - 1) - Tau^2 * (1 / (m - 1) + 1 / (n - 1)) + (1 - X) * (1 / (m * (n - 1)) + 1 / (n * (m - 1)))
data.frame(Tau = Tau, HM = sqrt(V_HM), Unb = sqrt(V_U), JK = sqrt(V_JK), trunc = 1 / (m * n))
}
A_data <- rpois(10, lambda = 10)
B_data <- rpois(12, lambda = 12)
test_that("Tau SEs agree with above.", {
Tau_check <- Tau_est(yA = A_data, yB = B_data)
unbiased <- Tau(A_data, B_data, SE = "unbiased", confidence = NULL, trunc_const = TRUE)
Hanley <- Tau(A_data, B_data, SE = "Hanley", confidence = NULL)$SE
expect_equal(Tau_check$Unb, unbiased$SE)
expect_equal(Tau_check$HM, Hanley)
expect_equal(Tau_check$trunc, unbiased$trunc)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.