Nothing
# Incidence Estimation Tools Copyright (C) 2015-2018, DST-NRF Centre of
# Excellence in Epidemiological Modelling and Analysis (SACEMA), Stellenbosch
# University and individual contributors.
#
# This program is free software: you can redistribute it and/or modify it under
# the terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version. This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details. You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
DM_FirstOrderTerms <- function(prevH, prevR, mdri, frr, bigt) {
fot_prevH <- (prevR - frr)/(((1 - prevH)^2) * (mdri - frr * bigt)) #E.G. d(I)/d(P_H)
fot_prevR <- prevH/((1 - prevH) * (mdri - frr * bigt))
fot_mdri <- (frr * prevH - prevR * prevH)/((1 - prevH) * ((mdri - frr * bigt)^2))
fot_frr <- (prevH * (bigt * prevR - mdri))/((1 - prevH) * ((mdri - frr * bigt)^2))
return(c(fot_prevH, fot_prevR, fot_mdri, fot_frr))
}
#' Precision and Required Sample Size
#'
#' Calculate precision of incidence estimate and required sample size for a
#' given precision of the incidence estimate.
#'
#' @param I Expected Incidence.
#' @param RSE_I Relative Standard Error of Incidence Estimate. If this is the desired output, set to "out".
#' @param PrevH Prevalence of HIV.
#' @param CR Coverage rate: probability (0-1) of being tested for recency when positive for HIV.
#' @param MDRI mean duration of recent infection in days (vector/integer).
#' @param RSE_MDRI Relative standard error of MDRI (vector/integer).
#' @param FRR False recent rate (vector/integer).
#' @param RSE_FRR Relative standard error of FRR (vector/integer).
#' @param BigT post-infection time cut-off for true vs. false recency. Default is 730 days.
#' @param DE_H Design effect of HIV prevalence test (vector/integer).
#' @param DE_R Design effect of recency test (vector/integer).
#' @param n Sample Size: Set to a hypothetical value if the desired output is RSE_I, othewise set to "out" to obtain required sample size.
#' @param step number of steps between minimum I and maximum I in the calculation of a range of output.
#' @return Either sample size necessary for a given precision under a given set of testing characteristics and a hypothetical prevalence/incidence scenario, or precision under a particular sample size scenario, with a given hypothetical prevalence/incidence scenario.
#' @details The package contains long form documentation in the form of vignettes that cover the use of the main fucntions. Use browseVignettes(package="inctools") to access them.
#'
#' This function summarizes performance of a recent infection test into a standard error of the incidence estimate, given the
#' estimated test properties and hypothetical survey context or the sample size necessary for a given level of precision.
#'
#' Up to two arguments can be specified as ranges, with the input parameter `step` specifying the number of increments
#' between the endpoints of the two ranges supplied under the argument name. This yields output for each step.
#' See the second and third example below
#' for an illustration of this output.
#'
#' Either the argument RSE_I or the argument n must be set to "out".
#'
#' @examples
#' incprecision(I = 0.015, RSE_I = 0.25, PrevH = 0.2, CR = 1,
#' MDRI = 200, RSE_MDRI = 0.05, FRR = 0.01, RSE_FRR = 0.2,
#' BigT = 730, DE_H = 1.1, DE_R = 1, n = 'out')
#'
#' incprecision(I = c(0.015,0.02), RSE_I = 0.25, PrevH = c(0.10,0.20),
#' CR = 1, MDRI = 200, RSE_MDRI = 0.05, FRR = 0.01, RSE_FRR = 0.2,
#' BigT = 700, DE_H = 1, DE_R = 1, n = 'out', step = 5)
#'
#' incprecision(I = 0.017, RSE_I = 'out', PrevH = c(0.10,0.20),
#' CR = 1, MDRI = 211, RSE_MDRI = 0.05, FRR = 0.009, RSE_FRR = 0.2,
#' BigT = 720, n = 5000, step = 5)
#' @export
incprecision <- function(I, RSE_I, PrevH, CR, MDRI, RSE_MDRI, FRR, RSE_FRR, BigT = 730,
DE_H = 1, DE_R = 1, n = "out", step = 5) {
var_list <- list(I = I, RSE_I = RSE_I, PrevH = PrevH, CR = CR, MDRI = MDRI, RSE_MDRI = RSE_MDRI,
FRR = FRR, RSE_FRR = RSE_FRR, BigT = BigT, DE_H = DE_H, DE_R = DE_R, n = n,
step = step)
# CHECK TO MAKE SURE ONLY TWO VARIABLES ARE ALLOWED TO VARY
max.list <- 0
for (i in 1:length(var_list)) {
if (length(var_list[[i]]) > 1) {
max.list = max.list + 1
}
if (max.list > 2) {
stop("Only a maximum of 2 variables are allowed to vary")
}
}
if (sum(var_list == "out") > 1) {
stop("Only one of the variables RSE_I or n can be requested at a time")
}
if (length(var_list[1:8]) < 8) {
stop("Not enough variables have been specified")
}
for (i in 1:12) {
if (length(var_list[[i]]) > 2 | length(var_list[[i]]) < 1) {
stop(paste("specify (only) min & max values for ", names(var_list)[i]),
sep = "")
}
}
if (any(I < 0)) {
stop("I must be >= 0")
}
if (any(RSE_I != "out" & RSE_I < 0)) {
stop("RSE_I must be >= 0 or 'out'")
}
if (any(PrevH <= 0)) {
stop("PrevH must be > 0")
}
if (any(PrevH > 1)) {
stop("PrevH must be <= 1")
}
if (any(CR <= 0)) {
stop("CR must be > 0")
}
if (any(CR > 1)) {
stop("CR must be <= 1")
}
if (any(RSE_MDRI < 0)) {
stop("RSE_MDRI must be >= 0")
}
if (any(FRR < 0)) {
stop("FRR must be >= 0")
}
if (any(FRR > 1)) {
stop("FRR must be <= 1")
}
if (any(RSE_FRR < 0)) {
stop("RSE_FRR must be >= 0")
}
if (any(DE_H < 0)) {
stop("DE_H must be >= 0")
}
if (any(DE_H < 1) | any(DE_R < 1)) {
warning("Design effect smaller than 1, implying complex sampling produces smaller standard errors than simple random sampling")
}
# for (i in c(1, 3, 4, 6:8)) { #(i in c(1, 3, 4, 6:8)) {
# if (is.numeric(var_list[[1:length(var_list[i])]]) > 0 & (sum(var_list[[i]] <=
# 1) != length(var_list[[i]]) | sum(var_list[[i]] >= 0) != length(var_list[[i]]))) {
# stop("Some input values are less than 0 or greater than 1")
# }
# }
# THESE CHECKS NEED TO BE FIXED
# if (sum(RSE_I != "out") > 0) {
# if (is.numeric(var_list[[1:length(var_list[2])]]) > 0 & (sum(var_list[[2]] <=
# 1) != length(var_list[[2]]) | sum(var_list[[2]] >= 0) != length(var_list[[2]]))) {
# stop("Some input values are less than 0 or greater than 1")
# }
# }
# above code does what below code does, only in 1 line. Which should we keep? if
# (is.numeric(I)>0) {stopifnot (I<=1 & I>=0)} if (is.numeric(RSE_I)>0) {stopifnot
# (RSE_I<=1 & RSE_I>=0)} if (is.numeric(PrevH)>0) {stopifnot (PrevH<=1 &
# PrevH>=0)} if (is.numeric(CR)>0) {stopifnot (CR<=1 & CR>=0)} if
# (is.numeric(MDRI)>0) {stopifnot (MDRI>=0)} if (is.numeric(RSE_MDRI)>0)
# {stopifnot (RSE_MDRI<=1 & RSE_MDRI>=0)} if (is.numeric(FRR)>0) {stopifnot
# (FRR<=1 & FRR>=0)} if (is.numeric(RSE_FRR)>0) {stopifnot (RSE_FRR<=1 &
# RSE_FRR>=0)} if (is.numeric(DE_H)>0) {stopifnot (DE_H>=1)} if
# (is.numeric(DE_R)>0) {stopifnot (DE_R>=1)} if (is.numeric(n)>0) {stopifnot
# (n>100)}
if (sum(BigT <= 182)) {
warning("BigT is smaller than half a year")
}
if (sum(BigT < MDRI) > 0) {
stop("MDRI cannot be greater than BigT")
}
if (sum(RSE_MDRI < 0.01) > 0) {
warning("RSE of estimated MDRI is less than 1%")
}
if (sum(FRR == 0) > 0) {
warning("Zero estimated FRR")
}
if (sum(FRR > 0.1) > 0) {
warning("Estimated FRR is greater than 10%")
}
if (sum(RSE_FRR > 0.3) > 0) {
warning("RSE of estimated FRR is greater than 30%")
}
if (sum(RSE_FRR < 0.05) > 0) {
warning("RSE of estimated FRR is less than 5%")
}
if (sum(I > 0.2) > 0) {
warning(paste("Possible error in incidence input.", max(I), " seems exceptionally high",
sep = ""))
}
if (sum(n < 1000) > 0) {
warning("Sample size is smaller than 1000")
}
# Need error that reads: 'ERROR: Test properties not consistent with test for
# recent infection'
######### THIS WHOLE SECTION HERE IS FOR IF ONE OR MORE OF THE VARIABLES IS ALLOWED TO
######### VARY####### CREATES TWO NULL VALUES, I BELIEVE FOR WHICH VARIABLES ARE TO VARY
vary1 <- NULL
vary2 <- NULL
# NOW FOR EACH VARIABLE IN LIST, IF IT'S ALLOWED TO VARY, THEN MAKE A MATRIX OF
# VALUES FOR THE STEP, FROM BEGINNING TO END, STEP NUMBER OF COLUMNS. IF THE
# VARIABLE IS THE FIRST ONE IN THIS LIST TO VARY THEN WE'RE CALLING THE VARIABLE
# IT'S SAME NAME, AS A MATRIX, AND IF IT'S THE SECOND, WE'RE CALLING THAT
# VARIABLE AS IT'S SAME NAME, BUT AS A TRANSPOSED MATRIX, SO NOW THE ROWS ARE THE
# STEPPING VALUES there's got to be a way to make this more efficient, like if
# two have been met, STOP...
if (length(I) == 2) {
I <- matrix(rep(seq(from = min(I), to = max(I), length.out = step), times = step),
ncol = step, nrow = step)
if (length(vary1) == 0) {
vary1 <- I
vary_name1 <- "I"
} else {
vary2 = I = t(I)
vary_name2 = "I"
}
}
if (length(RSE_I) == 2) {
RSE_I <- matrix(rep(seq(from = min(RSE_I), to = max(RSE_I), length.out = step),
times = step), ncol = step, nrow = step)
if (length(vary1) == 0) {
vary1 <- RSE_I
vary_name1 <- "I"
} else {
vary2 = RSE_I = t(RSE_I)
vary_name2 = "I"
}
}
if (length(PrevH) == 2) {
PrevH <- matrix(rep(seq(from = min(PrevH), to = max(PrevH), length.out = step),
times = step), ncol = step, nrow = step)
if (length(vary1) == 0) {
vary1 <- PrevH
vary_name1 <- "PrevH"
} else {
vary2 = PervH = t(PrevH)
vary_name2 <- "PrevH"
}
}
if (length(CR) == 2) {
CR <- matrix(rep(seq(from = min(CR), to = max(CR), length.out = step), times = step),
ncol = step, nrow = step)
if (length(vary1) == 0) {
vary1 <- CR
vary_name1 <- "CR"
} else {
vary2 = CR = t(CR)
vary_name2 = "CR"
}
}
if (length(MDRI) == 2) {
MDRI <- matrix(rep(seq(from = min(MDRI), to = max(MDRI), length.out = step),
times = step), ncol = step, nrow = step)
# IF THIS VARIABLE IS TO VARY, MAKE A MATRIX OF DIM (STEP*STEP) WHERE EACH COLUMN
# IS THE SEQUENCE OF VALUES
if (length(vary1) == 0) {
vary1 <- MDRI
vary_name1 <- "MDRI"
} else {
vary2 = MDRI = t(MDRI)
vary_name2 = "MDRI"
}
}
if (length(RSE_MDRI) == 2) {
RSE_MDRI <- matrix(rep(seq(from = min(RSE_MDRI), to = max(RSE_MDRI), length.out = step),
times = step), ncol = step, nrow = step)
# IF THIS VARIABLE IS TO VARY, MAKE A MATRIX OF DIM (STEP*STEP) WHERE EACH COLUMN
# IS THE SEQUENCE OF VALUES
if (length(vary1) == 0) {
vary1 <- RSE_MDRI
vary_name1 <- "RSE_MDRI"
} else {
vary2 = RSE_MDRI = t(RSE_MDRI)
vary_name2 = "RSE_MDRI"
}
}
if (length(FRR) == 2) {
FRR <- matrix(rep(seq(from = min(FRR), to = max(FRR), length.out = step),
times = step), nrow = step, ncol = step)
if (length(vary1) == 0) {
vary1 <- FRR
vary_name1 <- "FRR"
} else {
vary2 = FRR = t(FRR)
vary_name2 = "FRR"
}
}
if (length(RSE_FRR) == 2) {
RSE_FRR <- matrix(rep(seq(from = min(RSE_FRR), to = max(RSE_FRR), length.out = step),
times = step), nrow = step, ncol = step)
if (length(vary1) == 0) {
vary1 <- RSE_FRR
vary_name1 <- "RSE_FRR"
} else {
vary2 = RSE_FRR = t(RSE_FRR)
vary_name2 = "RSE_FRR"
}
}
if (length(BigT) == 2) {
BigT <- matrix(rep(seq(from = min(BigT), to = max(BigT), length.out = step),
times = step), nrow = step, ncol = step)
if (length(vary1) == 0) {
vary1 <- BigT
vary_name1 <- "BigT"
} else {
vary2 = BigT = t(BigT)
vary_name2 = "BigT"
}
}
if (length(DE_H) == 2) {
DE_H <- matrix(rep(seq(from = min(DE_H), to = max(DE_H), length.out = step),
times = step), nrow = step, ncol = step)
if (length(vary1) == 0) {
vary1 <- DE_H
vary_name1 <- "DE_H"
} else {
vary2 = DE_H = t(DE_H)
vary_name2 = "DE_H"
}
}
if (length(DE_R) == 2) {
DE_R <- matrix(rep(seq(from = min(DE_R), to = max(DE_R), length.out = step),
times = step), nrow = step, ncol = step)
if (length(vary1) == 0) {
vary1 <- DE_R
vary_name1 <- "DE_R"
} else {
vary2 = DE_R = t(DE_R)
vary_name2 = "DE_R"
}
}
if (length(n) == 2) {
n <- matrix(rep(seq(from = min(n), to = max(n), length.out = step), times = step),
nrow = step, ncol = step)
if (length(vary1) == 0) {
vary1 <- n
vary_name1 <- "n"
} else {
vary2 = n = t(n)
vary_name2 = "n"
}
}
######### END SECTION FOR IF ONE OR MORE OF THE VARIABLES IS ALLOWED TO VARY#######
# NOW MAKE VARIABLES IN DAYS TO BE IN UNITS YEARS
if (is.numeric(MDRI)) {
MDRI <- MDRI/365.25
}
if (is.numeric(BigT)) {
BigT <- BigT/365.25
}
PrevR <- ((I * (1 - PrevH) * (MDRI - FRR * BigT))/PrevH + FRR)
out2 <- PrevHR <- round(PrevH * PrevR, digits = 5) #Prev.HIV&recent
out3 <- PrevHnR <- round(PrevH - PrevHR, digits = 5) #Prev.HIV&nonrecent
fot <- DM_FirstOrderTerms(PrevH, PrevR, MDRI, FRR, BigT)
# if the output of each term of DM_FirstOrderTerms is univariate, do one thing,
# otherwise, do another...
if (sum(lengths(var_list) > 1) == 0 | length(fot) == 4) {
fot_PrevH <- fot[1]
fot_PrevR <- fot[2]
fot_MDRI <- fot[3]
fot_FRR <- fot[4]
} else {
# this needs to reflect what happens to the output fot matrix which depends on
# which variables are allowed to vary.
if (length(I) > 1 & sum(lengths(var_list)) == 14) {
fot_PrevH <- matrix(fot[1:(step * step)], nrow = step, ncol = step, byrow = FALSE)
fot_PrevR <- fot[(step * step + 1)] #things change if and only if only Incidence is allowed to vary.
fot_MDRI <- matrix(fot[((step * step) + 2):(((step * step) * 2 + 1))],
nrow = step, ncol = step)
fot_FRR <- matrix(fot[(((step * step) * 2 + 2)):length(fot)], nrow = step,
ncol = step)
} else {
# here is the situation if only a variable besides incidence is allowed to vary,
# or any two parameters are allowed to vary.
fot_PrevH <- matrix(fot[1:(step * step)], nrow = step, ncol = step, byrow = FALSE)
fot_PrevR <- matrix(fot[((step * step) * 1 + 1):(((step * step) * 2))],
nrow = step, ncol = step)
fot_MDRI <- matrix(fot[((step * step) * 2 + 1):(((step * step) * 3))],
nrow = step, ncol = step)
fot_FRR <- matrix(fot[((step * step) * 3 + 1):(((step * step) * 4))],
nrow = step, ncol = step)
}
}
# IF SAMPLE SIZE n IS THE OUPUT VARIABLE (SO PRECISION/RSE_I IS FIXED)
if (sum(n == "out") > 0) {
out4 <- RSE_I_inf_ss <- round(sqrt((fot_MDRI * RSE_MDRI * MDRI)^2 + (fot_FRR *
RSE_FRR * FRR)^2)/I, digits = 5) #RSE.I.inf.sample
out1 <- n <- ceiling(((fot_PrevH^2) * PrevH * (1 - PrevH) * DE_H + (fot_PrevR^2) *
(PrevR * (1 - PrevR) * DE_R/(CR * PrevH)))/((RSE_I^2 - RSE_I_inf_ss^2) *
I^2))
if (sum(((PrevH * (1 - PrevH))/n) * DE_H <= 0) > 0) {
stop("no sample size will meet input constraints")
}
out5 <- RSE_PrevH <- round(sqrt(((PrevH * (1 - PrevH))/n) * DE_H)/PrevH,
digits = 5) #RSE.PrevH
out6 <- RSE_PrevR <- round(sqrt(((PrevR * (1 - PrevR))/n * CR * PrevH) *
DE_R)/PrevR, digits = 5) #RSE.PrevR
out_names <- c("sample.size", "Prev.HIV.and.recent", "Prev.HIV.and.nonrecent",
"RSE.I.inf.sample", "RSE.PrevH", "RSE.PrevR")
}
# IF precision IS THE OUPUT VARIABLE (SO sample size n IS FIXED)
if (sum(RSE_I == "out") > 0) {
out4 <- RSE_I_inf_ss <- round(sqrt((fot_MDRI * RSE_MDRI * MDRI)^2 + (fot_FRR *
RSE_FRR * FRR)^2)/I, digits = 5)
out5 <- RSE_PrevH <- round(sqrt(((PrevH * (1 - PrevH))/n) * DE_H)/PrevH,
digits = 5)
out6 <- RSE_PrevR <- round(sqrt(((PrevR * (1 - PrevR))/n * CR * PrevH) *
DE_R)/PrevR, digits = 5)
out1 <- RSE_I <- round(sqrt(((fot_PrevH^2) * PrevH * (1 - PrevH) * DE_H +
(fot_PrevR^2) * (PrevR * (1 - PrevR) * DE_R/(CR * PrevH)))/(n * I^2) +
RSE_I_inf_ss^2), digits = 5)
out_names <- c("RSE_I", "Prev.HIV.and.recent", "Prev.HIV.and.nonrecent",
"RSE.I.inf.sample", "RSE.PrevH", "RSE.PrevR")
}
# if (sum(RSE_I > 0.50) >0) { warning('Implied RSE of incidence is greater than
# 50%') }
if (sum(lengths(var_list)) == 15) {
# if two variables are to vary
variable.1 <- vector(length = step)
variable.2 <- vector(length = step)
for (i in 1:step) {
variable.1[i] <- paste(vary_name1, "=", vary1[i, 1])
variable.2[i] <- paste(vary_name2, "=", vary2[1, i])
}
if (length(out1) > 1) {
out1 <- data.frame(out1)
row.names(out1) <- variable.1
colnames(out1) <- variable.2
}
if (length(out2) > 1) {
out2 <- data.frame(out2)
row.names(out2) <- variable.1
colnames(out2) <- variable.2
}
if (length(out3) > 1) {
out3 <- data.frame(out3)
row.names(out3) <- variable.1
colnames(out3) <- variable.2
}
if (length(out4) > 1) {
out4 <- data.frame(out4)
row.names(out4) <- variable.1
colnames(out4) <- variable.2
}
if (length(out5) > 1) {
out5 <- data.frame(out5)
row.names(out5) <- variable.1
colnames(out5) <- variable.2
}
if (length(out6) > 1) {
out6 <- data.frame(out6)
row.names(out6) <- variable.1
colnames(out6) <- variable.2
}
# out1<-data.frame(out1) out2<-data.frame(out2) out3<-data.frame(out3)
# out4<-data.frame(out4) out5<-data.frame(out5) out6<-data.frame(out6)
# row.names(out1) <-variable.2 colnames(out1) <-variable.1 row.names(out2)
# <-variable.2 colnames(out2) <-variable.1 row.names(out3) <-variable.2
# colnames(out3) <-variable.1 row.names(out4) <-variable.2 colnames(out4)
# <-variable.1 row.names(out5) <-variable.2 colnames(out5) <-variable.1
# row.names(out6) <-variable.2 colnames(out6) <-variable.1
}
if (sum(lengths(var_list)) == 14) {
# if just one variable is allowed to vary
variable.1 <- vector(length = step)
for (i in c(1:step)) {
variable.1[i] <- paste(vary_name1, "=", vary1[i, 1])
}
if (length(out1) > 1) {
out1 <- data.frame(variable.1, out1[, 1])
names(out1) <- c(vary_name1, out_names[1])
}
if (length(out2) > 1) {
out2 <- data.frame(variable.1, out2[, 1])
names(out2) <- c(vary_name1, out_names[2])
}
if (length(out3) > 1) {
out3 <- data.frame(variable.1, out3[, 1])
names(out3) <- c(vary_name1, out_names[3])
}
if (length(out4) > 1) {
out4 <- data.frame(variable.1, out4[, 1])
names(out4) <- c(vary_name1, out_names[4])
}
if (length(out5) > 1) {
out5 <- data.frame(variable.1, out5[, 1])
names(out5) <- c(vary_name1, out_names[5])
}
if (length(out6) > 1) {
out6 <- data.frame(variable.1, out6[, 1])
names(out6) <- c(vary_name1, out_names[6])
}
}
output <- list(out1, out2, out3, out4, out5, out6)
names(output) <- out_names
return(output)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.