# R/bgtvs.R In binGroup: Evaluation and Experimental Design for Binomial Group Testing

```"bgtvs" <- function
(n, s, y, conf.level = 0.95, alternative="two.sided", maxiter=100)
{

alternative<-match.arg(alternative, choices=c("two.sided", "less", "greater"))

# error checking:

if (length(y)!=length(n)||length(n)!=length(s))
{stop("vectors n, s, y must have exactly the same length")}

if ( any(y>n) )
{stop("values of y must be smaller than or equal to the corresponding values n") }

if(any(abs(round(n)-n) > 1e-07) | any(n<1))
{stop("number of groups n must be a vector of integer values > 0")}

if(any(abs(round(s)-s) > 1e-07) | any(s<1))
{stop("group sizes s must be a vector of integer values > 0")}

if(any(abs(round(y)-y) > 1e-07) | any(y<0))
{stop("number of positive groups y must be a vector of non-negative integer values > 0")}

if(conf.level<=0 | conf.level>=1)
{stop("conf.level must be a numeric value between 0, and 1, usually e.g. 0.95")}

m <- length(y)
alpha <- 1-conf.level

# # # Likelihood according to Hepworth 1996 (1):

likHepI <- function(n, s, y, p)
{
prod( choose(n,y) * (( 1-(1-p)^s )^y) * ( (1-p)^(s*(n-y)) )  )
}

# function to estimate the MLE:

estBGTvgs2<-function(n, s, y)
{
if (length(y)!=length(n)||length(n)!=length(s)) {stop("vectors n, s, y must have exactly the same length")}

# total number of individuals in the trial

# the current iteration method is sensitive
# to give wrong results for too small tolerances

tol=1e-10
total <- sum(n*s)

# Hepworth(1996), equation 2

mindiff <- function(n,s,y,p, total)
{total - sum( (s*y) / (1-(1-p)^(s)) )}

# check the case all y=0 (equa 2 not defined),
# else iterate p

if(all(y==0))
{esti<-0}

else
{
maxiter <- 100
dir<-1
crit<-numeric(length=maxiter)

if(any(y==0))
{esti<-0+tol}
else
{esti<-0}

crit[1] <- mindiff(n=n, s=s, y=y, p=esti, total=total)
step=0.5

for(i in 2:maxiter)
{
crit[i] <- mindiff(n=n, s=s, y=y, p=esti, total=total)
if(sign(crit[i-1])!=sign(crit[i]))
{dir <- dir*(-1); step <- step/2}
if(esti>1 | esti<0)
{dir<-dir*(-1)}
esti<-esti+dir*step
}
}
return(esti)
}

# # # Point estimate for the observed outcome:
# Pobs (Hepworth)

point.esti <- estBGTvgs2(n=n, s=s, y=y)

# # # calculate the MLE for all possible combinations of y1, ..., ym

possibley<-list()

for(i in 1:length(n))
{
possibley[[i]]<-0:n[i]
}
allComb=expand.grid(possibley)

MLEComb<-numeric(length=nrow(allComb) )

for (comb in 1:nrow(allComb))
{ytemp<-as.numeric(allComb[comb,])
MLEComb[comb] <- estBGTvgs2(n=n, s=s, y=ytemp)
}

out=cbind(allComb,MLEComb)

# # # order the event space by their associated MLE

outlower<-as.matrix( out[which(out\$MLEComb >= point.esti), 1:m] )

outupper <- as.matrix( out[which(out\$MLEComb <= point.esti), 1:m])

# Iteration functions for the upper and lower bound:

itupperlimit <- function(pstart, outupper, alpha, n, s, maxiter=maxiter)
{
dir<-1
step<-0.1

crit<-numeric(length=maxiter)

crit[1] <- alpha - sum(apply(X=outupper, MARGIN=1, FUN=likHepI, n=n, s=s, p=pstart))

p <- pstart+dir*step

for(i in 2:maxiter)
{
crit[i] <- alpha - sum(apply(X=outupper, MARGIN=1, FUN=likHepI, n=n, s=s, p=p))
if(sign(crit[i-1])!=sign(crit[i]))
{dir <- dir*(-1); step <- step/2}

p <- min(p+dir*step,1)
}

return(pu=p)
}

itlowerlimit <- function(pstart, outlower, alpha, n, s, maxiter=maxiter)
{
dir<-(-1)
step<-0.1

crit<-numeric(length=maxiter)

crit[1] <- alpha - sum(apply(X=outlower, MARGIN=1, FUN=likHepI, n=n, s=s, p=pstart))

p<-pstart+dir*step

for(i in 2:maxiter)
{
crit[i] <- alpha - sum(apply(X=outlower, MARGIN=1, FUN=likHepI, n=n, s=s, p=p))
if(sign(crit[i-1])!=sign(crit[i]))
{dir <- dir*(-1); step <- step/2}

p <- max(p+dir*step,0)
}

return(pl=p)
}

# # # Calculation of Bounds

if(alternative=="two.sided")
{
upper <- itupperlimit(pstart=point.esti, outupper=outupper, alpha=alpha/2, n=n, s=s, maxiter=100)
lower <- itlowerlimit(pstart=point.esti, outlower=outlower, alpha=alpha/2, n=n, s=s, maxiter=100)
}

if(alternative=="less")
{
upper <- itupperlimit(pstart=point.esti, outupper=outupper, alpha=alpha, n=n, s=s, maxiter=100)
lower <- 0
}

if(alternative=="greater")
{
upper=1
lower <- itlowerlimit(pstart=point.esti, outlower=outlower, alpha=alpha/2, n=n, s=s, maxiter=100)
}

out <- list(conf.int=c(lower, upper),
estimate=point.esti,
conf.level=conf.level,
alternative=alternative,
input=rbind("number of groups"=n,
"group size"=s,
"number of positive groups"=y)
)

class(out)<-"bgtvs"

return(out)

}

"print.bgtvs"<-function
(x,...)
{
args<-list(...)

conf<-round(100*x\$conf.level,2)

cat("Input","\n")
print(x\$input)
cat("\n")
cat("Exact",conf,"percent confidence interval","\n")
cat("\n")

if(is.null(args\$digits))
args\$digits <- 4
args\$x<-x[["conf.int"]]
do.call("print", args)
cat("\n")
cat("Point estimate:","\n")
cat("\n")
args\$x<-x[["estimate"]]
do.call("print", args)
}
```

## Try the binGroup package in your browser

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

binGroup documentation built on March 18, 2018, 2:04 p.m.