R function conf_int_Accuracy

 

R function

conf_int_Accuracy <- function(GMTobs, s, nDays = 2, nAnalysts = 2, nPlates = 2, nReplicates = 3, FoldDilution = 2, Threshold = 1, alpha = 0.05) {

  # INPUT:

  # GMTobs: vector of observed geometric means at each fold dilution

  # s: vector of the standard deviations of the log-transformed replicates

  # nDays: number of experimental days (default = 2)

  # nAnalysts: number of analysts performing the tests in validation (default = 2)

  # nPlates: number of plates used per analyst per day (default = 2)

  # nReplicates: number of measurements that each analyst performs per plate per day

  # FoldDilution: step of the serial dilution (default = 2)

  # Threshold: critical threshold of the log-difference between the observed and the true mean (default=1)

  # alpha: significance level (default=0.05)

  #

  # OUTPUT:

  # vector of Relative Accuracy calculated at each fold dilution and its confidence interval. In addition, a plot of Relative Accuracy vs. fold dilution (log2 units)

 

  # calculate the expected geometric mean vector

  GMTexp <- round(GMTobs[1] * 2^(seq(0, -(length(GMTobs) - 1), -1)), 0)

 

  # Relative Accuracy (or recovery)

  RA <- log(GMTobs / GMTexp)

 

  N <- nDays * nAnalysts * nPlates * nReplicates

  n1 = n2 = N

 

  # standard deviations

  s1 <- s  # vector of standard deviations of the log-transformed replicates

  fd <- round(FoldDilution^-seq(0, length(GMTobs) - 1), 4)

  s2 <- rep(s1[1], length(GMTobs))  # fd*s1[1]

 

  # pooled variance

  sp2 = ((n1 - 1) * s1^2 + (n2 - 1) * s2^2) / (n1 + n2 - 2)

 

  # standard error

  se <- sqrt((sp2 / n1) + (sp2 / n2))

 

  # t-statistic

  tstat <- qt(1 - alpha / 2, n1 + n2 - 2)

 

  # confidence interval

  lower <- RA - tstat * se

  upper <- RA + tstat * se

 

  # critical thresholds

  theta.low <- FoldDilution^-Threshold

  theta.up <- FoldDilution^Threshold

 

  # exponentiate to get back to original scale

  RA <- exp(RA)

  lower <- exp(lower)

  upper <- exp(upper)

 

  # plot RA versus fold dilution

  plot(log2(fd), RA, type = "o", pch = 16, lty = 2, xlab = "fold dilution (log2)", ylab = "Relative Accuracy", ylim = c(0, 3))

  lines(log2(fd), lower, col = "blue", lty = 4)

  lines(log2(fd), upper, col = "blue", lty = 4)

  grid()  

 

  # return RA and its confidence interval

  return(cbind(fd,RA, lower, upper))

}

 

This function calculates the Relative Accuracy (RA) and its confidence interval at each fold dilution and plots RA versus fold dilution (log2 units). The function takes the following inputs:

• GMTobs: vector of observed geometric means at each fold dilution
• s: vector of the standard deviations of the log-transformed replicates
• nDays: number of experimental days (default = 2)
• nAnalysts: number of analysts performing the tests in validation (default = 2)
• nPlates: number of plates used per analyst per day (default = 2)
• nReplicates: number of measurements that each analyst performs per plate per day
• FoldDilution: step of the serial dilution (default = 2)
• Threshold: critical threshold of the log-difference between the observed and the true mean (default=1)
• alpha: significance level (default=0.05)

The function returns a matrix containing the RA and its confidence interval for each fold dilution. The confidence interval is calculated using the t-distribution with n1+n2-2 degrees of freedom and a significance level of alpha. The function also plots RA versus fold dilution (log2 units) and the confidence interval.