# R/reduction_power.R In mdenwood/bayescount: Statistical Analyses and Power Calculations for Count Data and Faecal Egg Count Reduction Tests (FECRT)

#### Documented in reduction_powerreduction_pval

```#' @name reduction_power
#' @aliases reduction.power fecrt.power FECRT.power
#' @title Count reduction (eg FECRT) data power calculations
#'
#' @description
#' Power calculation based on BNB code and Christian approximation.
#'

reduction_pval <- function(pre_data, post_data, H0_1=0.9, H0_2=0.95, edt_pre=1, edt_post=1, prob_priors=c(1,1), k_change=1, true_k=NA, delta_method=TRUE, beta_iters=10^6){

stopifnot(beta_iters <= 10^8)
# Max 10^8 - 10^6 is pretty instant and gives the same results

stopifnot(all(abs(round(pre_data)-pre_data) < 10^-6, na.rm=TRUE))
stopifnot(all(abs(round(post_data)-post_data) < 10^-6, na.rm=TRUE))
stopifnot(length(prob_priors)==2 && all(prob_priors > 0))

edt_change <- edt_pre / edt_post
stopifnot(length(edt_change)==1 && !is.na(edt_change) && edt_change > 0)

stopifnot(length(k_change)==1 && !is.na(k_change) && k_change > 0)

usetruek <- TRUE
if(all(is.na(true_k))){
true_k <- 1
usetruek <- FALSE
}
stopifnot(length(true_k)==1 && !is.na(true_k) && true_k > 0)

results <- .C(C_fecrt_pee_wrap, pre=as.integer(pre_data), preN=as.integer(length(pre_data)), post=as.integer(post_data), postN=as.integer(length(post_data)), H0_1=as.numeric(H0_1), H0_2=as.numeric(H0_2), edt_change=as.numeric(edt_change), prob_priors=as.numeric(prob_priors), kchange=as.numeric(k_change), truek=as.numeric(true_k), usetruek=as.integer(usetruek), delta=as.integer(delta_method), beta_iters=as.integer(beta_iters), p_1=as.numeric(0), p_2=as.numeric(0))

return(list(p_1=results\$p_1, p_2=results\$p_2))

}

reduction_pvals <- function(N, pre_mean, reduction, pre_k=1, post_k=pre_k, paired=FALSE, iterations=10000, animal_k=1, H0_1=0.9, H0_2=0.95, preN=N, postN=N, edt_pre=1, edt_post=1, prob_priors=c(1,1), k_change=1, true_k=NA, delta_method=TRUE, beta_iters=10^6){

stopifnot(beta_iters <= 10^8)
# Max 10^8 - 10^6 is pretty instant and gives the same results
stopifnot(iterations <= 10^8)

usetruek <- TRUE
if(all(is.na(true_k))){
true_k <- 1
usetruek <- FALSE
}
stopifnot(length(true_k)==1 && !is.na(true_k) && true_k > 0)

if(paired){
stopifnot(preN==postN)
}else{
animal_k <- 1000  # Can be >100 -> means unpaired
}
edt_change <- edt_pre / edt_post

stopifnot(length(H0_1)==length(H0_2))
H0_N <- length(H0_1)
maxN <- max(preN, postN)

p_1=p_2 <- rep(0, iterations*H0_N)

pre <- rep(0, maxN)
post <- rep(0, maxN)

powres <- .C(C_fecrt_pvals, iters=as.integer(iterations), preN=as.integer(preN), postN=as.integer(postN), maxN=as.integer(maxN), premean=as.numeric(pre_mean), reduction=as.numeric(reduction), edt_change=as.numeric(edt_change), animalk=as.numeric(animal_k), prek=as.numeric(pre_k), postk=as.numeric(post_k), H0_1=as.numeric(H0_1), H0_2=as.numeric(H0_2), H0_N=as.integer(H0_N), prob_priors=as.numeric(prob_priors), kchange=as.numeric(k_change), truek=as.numeric(true_k), usetruek=as.integer(usetruek), delta=as.integer(delta_method), beta_iters=as.integer(beta_iters), predata=as.integer(pre), postdata=as.integer(post), p_1=as.numeric(p_1), p_2=as.numeric(p_2))

p_1 <- powres\$p_1
p_2 <- powres\$p_2
dim(p_1) <- c(iterations, H0_N)
dim(p_2) <- c(iterations, H0_N)

return(data.frame(p1=p_1, p2=p_2))

}

reduction_power <- function(N, pre_mean, reduction, pre_k=1, post_k=pre_k, paired=FALSE, iterations=10000, animal_k=1, H0_1=0.9, H0_2=0.95, preN=N, postN=N, edt_pre=1, edt_post=1, prob_priors=c(1,1), k_change=1, true_k=NA, delta_method=TRUE, beta_iters=10, lci_c=0.005, uci_c=0.995, lci_b=0.025, uci_b=0.975, dobson_priors=c(1,1), tval=2.048, psig=0.025){

stopifnot(beta_iters <= 10^8)
# Max 10^8 - 10^6 is pretty instant and gives the same results
stopifnot(iterations <= 10^8)

edt_change <- edt_pre / edt_post
stopifnot(length(edt_change)==1 && !is.na(edt_change) && edt_change > 0)

stopifnot(length(k_change)==1 && !is.na(k_change) && k_change > 0)

usetruek <- TRUE
if(all(is.na(true_k))){
true_k <- 1
usetruek <- FALSE
}
stopifnot(length(true_k)==1 && !is.na(true_k) && true_k > 0)

if(paired){
stopifnot(preN==postN)
}else{
animal_k <- 1000  # Can be >100 -> means unpaired
}

# Needs other checks

maxN <- max(preN, postN)
predata=postdata <- rep(0, maxN)
classifications <- rep(0, 3*3*3*3*2)

results <- .C(C_fecrt_power_comparison, iters=as.integer(iterations), preN=as.integer(preN), postN=as.integer(postN), maxN=as.integer(maxN), premean=as.numeric(pre_mean), reduction=as.numeric(reduction), edt_change=as.numeric(edt_change), animalk=as.numeric(animal_k), prek=as.numeric(pre_k), postk=as.numeric(post_k), H0_1=as.numeric(H0_1), H0_2=as.numeric(H0_2), lci_c=as.numeric(lci_c), uci_c=as.numeric(uci_c), lci_b=as.numeric(lci_b), uci_b=as.numeric(uci_b), dobson_priors=as.numeric(dobson_priors), tval=as.numeric(tval), psig=as.numeric(psig), prob_priors=as.numeric(prob_priors), kchange=as.numeric(k_change), truek=as.numeric(true_k), usetruek=as.integer(usetruek), delta=as.integer(delta_method), beta_iters=as.integer(beta_iters), predata=as.integer(predata), postdata=as.integer(postdata), classifications=as.integer(classifications))

classifications <- results\$classifications
dim(classifications) <- c(3,3,3,3,2)
dimnames(classifications) <- list(paste('WAAVP', c('Susceptible','Inconclusive','Resistant'), sep='_'), paste('Dobson', c('Susceptible','Inconclusive','Resistant'), sep='_'), paste('Beta', c('Susceptible','Inconclusive','Resistant'), sep='_'), paste('Delta', c('Susceptible','Inconclusive','Resistant'), sep='_'), c('100pcObs', 'Lt100pcObs'))

return(classifications)
}

reduction.power <- reduction_power
fecrt.power <- reduction_power
FECRT.power <- reduction_power
```
mdenwood/bayescount documentation built on Sept. 2, 2017, 11:33 p.m.