snapshot-adjusted.Rmd
library(XSRecency)
In this vignette, we show how to use the snapshot and adjusted estimators for cross-sectional incidence estimation. First, we will simulate cross-sectional data using functions from within this package. Then we illustrate how to use the estimation functions.
# Create an test recent function that approximates the following window and shadow period
# getTestRecentFunc() takes in window and shadow period in days as arguments
phi.func <- getTestRecentFunc(window=200, shadow=191)
ts <- seq(0, 10, by=0.01)
plot(phi.func(ts) ~ ts, type='l')
# Calculate the MDRI at 2 years for this function (we know that the window is 200/365.25)
window <- 200/365.25
mdri <- integratePhi(phi.func, maxT=2)
mdri
#> [1] 0.5363007
sim <- simCrossSect(phi.func=phi.func,
incidence_type="constant", prevalence=0.29, baseline_incidence=0.032)
set.seed(10)
data <- sim(5000)
head(data)
#> di ui ri
#> 1 0 NA NA
#> 2 0 NA NA
#> 3 0 NA NA
#> 4 0 NA NA
#> 5 0 NA NA
#> 6 0 NA NA
tail(data)
#> di ui ri
#> 4995 1 5.4014360 0
#> 4996 1 10.3831010 0
#> 4997 1 12.6781044 0
#> 4998 1 0.5060816 1
#> 4999 1 1.0612372 0
#> 5000 1 7.8467520 0
Records with are HIV-negative individuals, and are positive individuals. The variables and are infection times (based on epidemiological parameters supplied) and recent infection indicators (based on the phi function), respectively.
estSnapshot(n_r=sum(data$ri, na.rm=T), n_n=sum(!data$di), n_p=sum(data$di), n=5000,
mu=window, mu_var=0)
#> $est
#> [1] 0.03369465
#>
#> $var
#> [1] 1.778887e-05
In some scenarios, not all individuals who test positive for HIV will
receive a test. Here we show how to apply the estimator in these cases,
using the argument n_p_test
. The same argument can be used
in the estAdjusted
function.
# Remove 50% of the recency tests to simulate what would happen if only 50% of
# HIV positive individuals received a recency test
data$mask <- rbinom(n=5000, size=1, prob=0.5)
data$ri_mask <- data$ri
data[data$mask == 1, "ri_mask"] <- NA
estSnapshot(n_r=sum(data$ri_mask, na.rm=T), n_n=sum(!data$di), n_p=sum(data$di), n=5000,
mu=window, mu_var=0, n_p_test=sum(data$di & (!is.na(data$ri_mask))))
#> $est
#> [1] 0.03823064
#>
#> $var
#> [1] 3.89342e-05
# Suppose the FRR was known to be 0.1%
estAdjusted(n_r=sum(data$ri, na.rm=T), n_n=sum(!data$di), n_p=sum(data$di), n=100,
omega=mdri, omega_var=0, beta=0.001, beta_var=0, big_T=2)
#> $est
#> [1] 0.0337468
#>
#> $var
#> [1] 1.755951e-05