---
title: "Error rates for antimicrobial performance standards"
output: word_document
author: "COST: Biofilm Regulatory Toolbox"
date: "8 Sept 2026"
editor_options: 
  chunk_output_type: inline
---


# Summary

This R Markdown file calculates pass-errors and fail-errors for the test method applied to 1 microbe given stakeholder specifications for: 

(1) what mean LR  (<=5 for STM), repeatability SD and reproducibility SD are for the the ineffective product; 

(2) what mean LR (>=8.5 for STM), repeatability SD and reproducibility SD are for highly effective products 

(3) a performance standard (PS) that an antimicrobial product must pass as imposed by a regulatory authority, which requires

* How many labs (1 for the STM), how many tests at each lab (3 for the STM, with each test applied to a different batch of product) 

* the LR that is  required per test (>=6 for the STM) 

(4) The repeatability SD and reproducibility SD for the product being tested.

# How to cite
Parker AE “KSA-SM-18—DRAFT R code to assess error rates of an antimicrobial test method,” Knowledge Sharing Articles, Center for Biofilm Engineering at Montana State University, Bozeman, MT. October 23, 2023 [E-publication] 

KSA-SM-18 is applied to the Single Tube Method 


# Background

Fitting linear mixed effects models (LMM) to control and LR data is described in:

* M. Hamilton, G. Hamilton, D. Goeres, and A. Parker. Guidelines for the Statistical Analysis of a Collaborative Study of a Laboratory Disinfectant Product Performance Test Method. JAOAC International 96(5):1138-1151, September, 2013. Appendices with R code.


Pass and Fail error calcs are described in:

* A. Parker, M. Hamilton, and S. Tomasino. A Statistical Model for Assessing Performance Standards for Quantitative and Semi-quantitative Disinfectant Test Methods. <i>JAOAC International</i>, 97(1):58-67, 2014. Appendices with R code.


Frown of variances and How to pick which SD to put into pass and fail error calculations is described by:

* Parker, Hamilton, Goeres.  Reproducibility of antimicrobial test methods.  Scientific Reports  8:12531, 2018.



# Housekeeping

## R packages needed
```{r}
library(lme4)      # lmm fits
library(mvtnorm)   # PS calcs
```
## R functions needed

GenCorrMatrix() and GenCorrMatrix_Microbe() are in Appendix to: Parker et al, 2014 (see Background section) available at 
https://math.montana.edu/parker/projects/documents/Appendices%20for%20Parker2014_JAOAC.pdf 


### GenCorrMatrix()
```{r}
# An R function to generate the correlation matrix for LRs from mulitple labs and tests  for a single microbe
GenCorrMatrix<-function(NumTests,NumLabs=1,VarLab,VarExp)
{Z = matrix(0,NumTests,(NumTests+1))
Z[,1]=1
for (i in 1:NumTests)
{Z[i,i+1]=1}
Psi=diag(c(rep(VarLab,1),rep(VarExp,NumTests)))
Vblk=Z%*%Psi%*%t(Z)
N = NumTests*NumLabs # total number of data points across all labs and exps
V=matrix(0,N,N)
for (i in 0:(NumLabs-1))
{index=(i*NumTests+1):((i+1)*NumTests)
V[index,index]=Vblk # construct the covariance matrix first
}
V = V/V[1,1]
return(V)
} ##### End of function
```

### GenCorrMatrix_Microbe()

```{r}
# An R function to generate the correlation matrix for tests of two microbe
# NumTests is a 2x1 vector, NumTests[1] specifies number of tests for Pa in each lab
# NumTests[2] sepcifies number of tests for Sa in each lab
# VarLab is a 2x1 vector, VarLab[1] is for Pa, VarLab[2] is for Sa
# VarExp is a 2x1 vector, VarExp[1] is for Pa, VarExp[2] is for Sa
GenCorrMatrix_Microbe<-function(NumTests,NumLabs,CorrMicrobe,VarLab,VarExp)
{PaBlock=GenCorrMatrix(NumTests[1],1,VarLab[1],VarExp[1])
SaBlock=GenCorrMatrix(NumTests[2],1,VarLab[2],VarExp[2])
Vblk = matrix(CorrMicrobe,sum(NumTests),sum(NumTests))
index1=1:NumTests[1]
index2=(NumTests[1]+1):sum(NumTests)
Vblk[index1,index1]=PaBlock
Vblk[index2,index2]=SaBlock
N = sum(NumTests)*NumLabs # total number of data points across all labs and exps
V=matrix(0,N,N)
for (i in 0:(NumLabs-1))
{index=(i*sum(NumTests)+1):((i+1)*sum(NumTests))
V[index,index]=Vblk}
return(V)
} ##### End of function
```


### diagRES()
```{r}
# An R function to assess standardized residuals from a mixed effects model 
diagRES <- function(m)
{
r = resid(m,type="pearson")   # standardized residuals
par(mfrow = c(2,2))

hist(r,freq=FALSE,main="Density Plot",xlab="St Residuals")
lines(density(r))
boxplot(r,main="Boxplot",ylab="St Residuals")
xy=qqnorm(r,main="Normal Plot",ylab="St Residuals")
qqline(r)
# Check: homogeneity of variance of st res
plot(fitted(m),r,main="Res Vs. Fits",ylab="St Residuals")   # raw residuals
abline(0,0)
}
```

### GetSattDF()
```{r}
# Calculate the Satterthwaite DF for reproducibility SD calculated from multi-lab study with 
# I labs 
# J tests within each lab
# Var.Lab = variance between labs (calculated by LMM, usually = (MSlab - MSE)/J when applying MEthod of Moments when data are balanced)
# Var.Exp = variance within labs (calcualted by LMM, = MSE when applying Method of Moments)

##############
GetSattDF = function(I,J,Var.Lab,Var.Exp) {
  N = I*J
  SR = sqrt(Var.Lab + Var.Exp)
  R = max(0,Var.Lab/Var.Exp) 
  gtop = (R + 1)^2
  gbot = (((R + (1/J))^2)/(I - 1)) + ((1 -(1/J))/N)
  df = gtop/gbot
  return(df)
}
```


# USER Inputs

User specifies: what data are there? LRs only?   Or LRs and controls?   What's future required PS?

```{r}
# User choose LR values for PS assessment
LR.Req = 6

LR.Ineff = 5            # Mean LR for ineffective products
Sr.Ineff = 1.1889         # Repeatability SD at ineffective level
SR.Ineff = 2.07    # Reproducibility SD at ineffective level

LR.High = 8.5       # Mean LR for highly effective products
Sr.High = 0.524         # Repeatability SD at highly effective level
SR.High = 0.784     # Reproducibility at High efficacy level

# In multi-lab study used to assess method
# User reports the number of labs and number of tests at each lab 
Num.Labs = 7    # Number of Labs in the collaborative Study of the method.  This is I in GetSattDF( )
Num.Exp = 3   # Number of tests at each lab in the collaborative study of the method.  This J in GetSattDF( ), So there's a total of Num.Labs*Num.Exp tests in study for each treatment

# In required PS, how many Labs and how many tests must be performed
L = 1 # Number of labs
K = 3 # Number of tests in each of the labs
```
 
# PS calcs

The following code is taken from Appendix to Parker et al. 2014 (see Background) available at 
https://math.montana.edu/parker/projects/documents/Appendices%20for%20Parker2014_JAOAC.pdf

$\alpha$ output below represents Type I pass error rate

$\beta$ output below represents Type II fail error rate



## 1 species, 1 test

```{r}
Vlab.Ineff = SR.Ineff^2 - Sr.Ineff^2
Vexp.Ineff = Sr.Ineff^2
df.Ineff = GetSattDF(I=Num.Labs,J=Num.Exp,Vlab.Ineff,Vexp.Ineff)
t.Ineff = (LR.Req - LR.Ineff)/SR.Ineff
alpha1 = 1-pt(t.Ineff,floor(df.Ineff))

Vlab.High = SR.High^2 - Sr.High^2
Vexp.High = Sr.High^2
df.High = GetSattDF(I=Num.Labs,J=Num.Exp,Vlab.High,Vexp.High)
t.High = (LR.Req - LR.Ineff)/SR.High
lambda1 = (LR.High-LR.Ineff)/SR.High  # non-centrality param
beta1 = pt(t.High,floor(df.High),ncp=lambda1) # for P. aeruginosa

# Summary of error rates based on 1-test PS
data.frame(ErrorType=c("Pass","Fail"),ErrorRate=c(alpha1,beta1),PS.Tests=1,PS.LR.req=LR.Req,Product=c("Ineffective","High"),mean.LR=c(LR.Ineff,LR.High),SR=c(SR.Ineff,SR.High),Sr=c(Sr.Ineff,Sr.High),t=c(t.Ineff,t.High),df=c(df.Ineff,df.High),lambda=c(0,lambda1))
```

## 1 species, mulitple tests

```{r}
ICC = 1 - c(Sr.Ineff,Sr.High)^2/c(SR.Ineff,SR.High)^2

TotTests = K*L
(R.Ineff = GenCorrMatrix(K,L,VarLab=Vlab.Ineff,VarExp=Vexp.Ineff))   # Under Ho, the correlation matrix for the tests to be conducted to satisfy a PS

alphaKL = pmvt(lower=rep(t.Ineff, TotTests),upper=rep(Inf, TotTests ),df=floor(df.Ineff),corr=R.Ineff)

(R.High = GenCorrMatrix(K,L,VarLab=Vlab.High,VarExp=Vexp.High))   # Under Ha, the correlation matrix for the tests to be conducted to satisfy a PS

betaKL = 1-pmvt(lower=rep(t.High, TotTests),upper=rep(Inf,TotTests),df=floor(df.High),corr=R.High,delta=lambda1)

# Summary of error rates based on mulitple test PS with one bug
data.frame(ErrorType=c("Pass","Fail"),ErrorRate=c(alphaKL,betaKL),PS.Bugs=1,PS.Labs=L,PS.Tests=K,PS.LR.req=LR.Req,Product=c("Ineffective","High"),mean.LR=c(LR.Ineff,LR.High),SR=c(SR.Ineff,SR.High),Sr=c(Sr.Ineff,Sr.High),ICC,t=c(t.Ineff,t.High),df=c(df.Ineff,df.High),lambda=c(0,lambda1))
```

## 2 species, mulitple tests

This example code assumes the same number of tests and variance results for species 1 and species 2, but this is not necessary.  If there are different multi-lab data for each species, then in the code below simply set:

* K_1 to the number of tests in each lab for species 1, and K_2 to the number of tests in each lab for species 2.

* Vlab.Ineff.Species1, Vlab.Ineff.Species2, Vlab.High.Species1, Vlab.High.Species2 to the estimates from the multi-lab data.

* Vexp.Ineff.Species1, Vexp.Ineff.Species2, Vexp.High.Species1, Vexp.High.Species2 to the estimates from the multi-lab data.

```{r}
K_1 = K # Number of species 1 tests in each of the labs
K_2 = K # Number of species 2 tests in each of the labs

Vlab.Ineff.Species1 = Vlab.Ineff  
Vlab.Ineff.Species2 = Vlab.Ineff
Vexp.Ineff.Species1 = Vexp.Ineff
Vexp.Ineff.Species2 = Vexp.Ineff

t.Ineff_2 = t.Ineff
lambda1_2 = lambda1
# Under Ho, the correlation matrix for the tests to be conducted to satisfy a PS
(R.Ineff = GenCorrMatrix_Microbe(c(K_1,K_2),L,CorrMicrobe=0.25,VarLab=c(Vlab.Ineff.Species1, Vlab.Ineff.Species2),VarExp=c(Vexp.Ineff.Species1, Vexp.Ineff.Species2)))

vecblk = c(rep(t.Ineff,K_1),rep(t.Ineff_2,K_2))
lower.vec = rep(vecblk,L)
deltablk = c(rep(lambda1,K_1),rep(lambda1_2,K_2))
delta.vec = rep(deltablk,L)
TotTests = (K_1 + K_2)*L
alphaKL = pmvt(lower=lower.vec,upper=rep(Inf,TotTests),df=floor(df.Ineff),corr=R.Ineff)

Vlab.High.Species1 = Vlab.High  
Vlab.High.Species2 = Vlab.High
Vexp.High.Species1 = Vexp.High
Vexp.High.Species2 = Vexp.High

t.High_2 = t.High
lambda1_2 = lambda1
(R.High = GenCorrMatrix_Microbe(c(K_1,K_2),L,CorrMicrobe=0.25,VarLab=c(Vlab.High.Species1, Vlab.High.Species2),VarExp=c(Vexp.High.Species1, Vexp.High.Species2)))

vecblk = c(rep(t.High,K_1),rep(t.High_2,K_2))
lower.vec = rep(vecblk,L)
deltablk = c(rep(lambda1,K_1),rep(lambda1_2,K_2))
delta.vec = rep(deltablk,L)

betaKL = 1-pmvt(lower=lower.vec,upper=rep(Inf,TotTests),df=floor(df.High),corr=R.High,delta=delta.vec)

# Type I pass error
alphaKL

# Type II fail error
betaKL
```