Analyse test repeatability

This utility evaluates the repeatability of a test, based on the results of repeated testing on a series of 1 or more samples for a test producing a continuous outcome.

Data required is a series of test results for a panel of one or more samples. The first row of the data is expected to contain column headers identifying each of the sample series. The first column is assumed to contain identifiers for each test occasion (for example plate numbers, dates or just a sequential number). Testing data is assumed to be temporally ordered and matched across samples. Additional columns (after the first) are assumed to contain test results for the samples to be evaluated.

Outputs include (for each sample):

Numerical and graphical summaries of testing results;
Calculation of warning and critical limits and of numbers and percentages of occasions when these limits are exceeded;
Time series analysis, including plots of raw data, moving average and lowess trend line;
Summary of regression analysis for univariate linear trend; and
Autocorrelation plot.

The width of the window for calculating moving averages can be set by the user.
Warning and critical limits can be either:

Calculated from the data (mean +/- 2 and 3 standard deviations); or
Calculated from mean and standard deviation values specified by the user (mean +/- 2 and 3 standard deviations); or
Specified explicitly by the user.

Example data available here

No results

No example available
Repeatability_data_demo.xls

No references available

				######################################
# Program to summarise continuous data from time series for test evaluation
######################################

rm(list = ls())
# cat("
Test:",length(commandArgs()))
test<- ifelse(length(commandArgs()) < 3, TRUE, FALSE)
fpath<- ifelse(test, "webRootUrl", "rtoolsPath")
# cat(test)

# load header scripts
  source(paste(fpath, "R/epi_head.R", sep = ""))
  source(paste(fpath, "R/HTMLStream.R", sep = ""))
  source(paste(fpath, "R/epitools_functions.r", sep = ""))

# extract command arguments 1 = conf, 2 = digits, 3 = ma window, 4  = limits option (1 = from data, 2 = input mean&sd, 3 = input limits
# cat("
", a0[8:13])
    a1<- type.convert(a0[9:12])

# cat("
", a1)

lim<- a1[1]
ma.window<- a1[2]
conf<- a1[3]
digits<- a1[4]
# cat("
", a0[14])
spec.lim<- 0

if (lim > 1) {
#  if (test) {
      spec.lim<- c(0.5, 1, 1.5, 0.3, 0.4, 0.5)
#  } else { 
#      spec.lim<- as.numeric(strsplit(a0[14],",")[[1]])
#  }
}             
# spec.lim<- c(1.2, 1, 0.5, 0, -0.2, 1.5, 1.3, 1, 0.8, 0.5, 2, 1.8, 1.5, 1.2, 1)
filename<- digest(Sys.time())
test.name<- ifelse(test, "Test data", a0[8])              #test name
data.file<- ifelse(test, paste(fpath, "docs/ts_data_1.txt", sep = ""), a0[13])              #test name
tmp.path<- paste(fpath, "tmp/", sep = "")
tmp.file<- paste(fpath, "tmp/", filename, sep = "")
sinkfile<- paste(fpath, "tmp/", filename, ".txt", sep="")   # fpath,
# cat("
", data.file)

# read data from file
size<- file.info(data.file)$size
# cat("
", size)
if (length(spec.lim) == 1 && lim > 1) cat("
No critical limits entered")
if (size == 0) cat("
No Data entered")
data.list<- read.csv(data.file)

#	cat("
", spec.lim, lim)
# read testing data
#if (test) { 
#  data.file<- paste(fpath, "ts_data.csv", sep = "")              
#  data.list<- read.csv(data.file)
#} else {
#  data.file<- a0[12]             
#  data.list<- read.delim(data.file)
#}
# str(data.list)
samples<- ncol(data.list)-1
tests<- nrow(data.list)
#	cat("
", tests)
if (length(spec.lim) < samples*5 && lim == 3) cat("
Critical limits incomplete")
if (length(spec.lim) < samples*2 && lim == 2) cat("
Critical limits incomplete")
if (lim > 1) spec.lim<- matrix(spec.lim, nrow = length(spec.lim)/samples, ncol = samples, byrow = T)
# cat("
", spec.lim[1,], "
", spec.lim[2,])
 sink(sinkfile)

# summarise data
# values<- na.omit(data.list[,2])
# temp<- summarise.cont(data.list[,2], pclist = c(0.05, 0.25, 0.5, 0.75, 0.95), conf = conf, digits = digits)
if (samples == 1) {
    temp<- summarise.cont(data.list[,2], pclist = c(0.05, 0.25, 0.5, 0.75, 0.95), conf = conf, digits = digits)
    summary.results<- array(temp, dim = c(length(temp), 1))
    rownames(summary.results) = names(temp)
    colnames(summary.results)<- colnames(data.list)[2]
} else {
    summary.results<- apply(data.list[,2:ncol(data.list)], FUN = summarise.cont, MARGIN = 2, pclist = c(0.05, 0.25, 0.5, 0.75, 0.95), conf = conf, digits = digits)
}
# array(temp, dim = c(length(temp), 1))

# calculate numbers exceeding limits
limits<- array(0, dim = c(5, 3, samples))
colnames(limits)<- c("Value", "Number exceeding", "Proportion exceeding")
if (lim == 3) {
    rownames(limits)<- c("Upper critical limit", "Upper warning limit", "Expected value", "Lower warning limit", "Lower critical limit")
} else {
    rownames(limits)<- c("+3 sd", "+2 sd", "Mean", "-2 sd", "-3 sd")
}
dimnames(limits)[[3]]<- colnames(summary.results)
temp<- c(3, 2, 0, -2, -3)
n<- summary.results[13,]
mean.vals<- summary.results[8,]
sd.vals<- summary.results[9,]
# cat("
", test)

for (s in 1:samples) {
    if (lim == 1) {
        limits[1:5, 1, s]<- mean.vals[s] + temp * sd.vals[s]
    } else if (lim == 2) {
        limits[1:5, 1, s]<- spec.lim[1, s] + temp * spec.lim[2, s]
    } else {
        limits[1:5, 1, s]<- spec.lim[, s]
    } # end of if else
        
    for (x in 1:5) {
        if (temp[x] > 0) {
            limits[x, 2, s]<- length(which(data.list[,s+1] > limits[x, 1, s]))
        } else if (temp[x] < 0) {
            limits[x, 2, s]<- length(which(data.list[,s+1] < limits[x, 1, s]))
        }
    }
    limits[, 3, ]<- limits[, 2, s]/n[s]
} # end of s loop

# zoo.data<- zoo(data.list[,2])

# graph and regression summary for each column
lm.summary<- vector("list", samples)
for (s in 1:samples) {
    # create time-series object and plot
    ts.data<- ts(data.list[,s+1], start = 1, frequency = 1)
    graphfile<- paste(tmp.file, "_scc_", colnames(data.list)[s+1], ".png", sep="")
    OpenGraphOutput(graphfile, pointsize = 12, ht = 6, wd = 8)
        plot(ts.data, x = 1:length(ts.data), type = "n", ylim = c(min(limits[, 1, s], ts.data, na.rm = T), 
                      max(limits[,1, s], ts.data, na.rm = T)), xlab = "Run", ylab = "Test result", 
                      main = "Statistical control chart", xaxt = "n")
        axis(side = 1, labels = data.list[,1], at = 1:length(ts.data))
        mtext(colnames(data.list)[s+1], side = 3, line = 0.5, cex = 0.8)
        rect(xleft = 1, xright = length(ts.data), ytop = max(limits[, 1, s], ts.data, na.rm = T), ybottom = min(limits[, 1, s], ts.data, na.rm = T), col = "pink", border = "red")
        rect(xleft = 1, xright = length(ts.data), ytop = limits[1, 1, s], ybottom = limits[5, 1, s], col = "yellow", border = "orange")
        rect(xleft = 1, xright = length(ts.data), ytop = limits[2, 1, s], ybottom = limits[4, 1, s], col = "lightgreen", border = "green")
        abline(a = limits[3, 1, s], b = 0, col = "blue")
        abline(a = limits[1, 1, s], b = 0, col = "red")
        abline(a = limits[5, 1, s], b = 0, col = "red")
        abline(a = limits[2, 1, s], b = 0, col = "orange")
        abline(a = limits[4, 1, s], b = 0, col = "orange")
        lines(ts.data, col = "darkgreen", ylim = c(min(limits[, 1, s], ts.data, na.rm = T), max(limits[, 1, s], ts.data, na.rm = T)), type = "o", pch = "x", cex = 0.6)
        
        # regression line
        index<- c(1:nrow(data.list))
        lm.ts<- lm(data.list[,s+1] ~ index)
        # abline(lm.ts, lty = "dashed", col = "purple", lwd = 2)
        lines(lowess(na.omit(data.list[,s+1]), f = 0.5), lty = "dashed", col = "purple", lwd = 2)
# cat("
", test)
        # moving average
        mov.avg<- rollapply(ts.data, ma.window, FUN = mean, na.pad = T)
        # ma.start<- ceiling(ma.window/2)
        lines(y = mov.avg, x = 1:nrow(data.list), col = "brown", lty = "solid")
        legend("bottomright", c("Values", paste(ma.window, "-point moving avge", sep = ""), "Trend line (lowess)"),
              col = c("darkgreen", "brown", "purple"), lty = c("solid", "solid", "dashed"), lwd = c(1, 1, 2), cex = 0.7)

    CloseGraphOutput("R")
    
    # plot ACFs 
    graphfile<- paste(tmp.file, "_acf_", colnames(data.list)[s+1], ".png", sep="")
    OpenGraphOutput(graphfile, pointsize = 12, ht = 6, wd = 8)
        ts.acf<- acf(ts.data, na.action = na.pass, main = "Autocorrelation plot")
        mtext(colnames(data.list)[s+1], side = 3, line = 0.5, cex = 0.8)
    CloseGraphOutput("R")

# summarise results
    # save regression summary
    lm.summary[[s]]<- capture.output(summary(lm.ts))
    lm.summary[[s]]<- gsub(" ", " ", lm.summary[[s]])
} # end of s loop

# boxplot of summary results
graph.file<- paste(tmp.file, "_boxplot.png", sep = "")
OpenGraphOutput(graph.file, pointsize = 12, ht = 6, wd = 4+(samples-1))
    boxplot(data.list[,2:(samples+1)], col = "green", border = "darkgreen", ylab = "Test Result",
            main = "Distribution of sample results")
  mtext(test.name, side = 3, line = 0.5, cex = 1)
CloseGraphOutput("R")

# strip plot of infected and uninfected
  Title<-c("Distribution of sample results")
  graphfile<- paste(fpath, "tmp/", filename, "_stripchart.png", sep="")
  OpenGraphOutput(graphfile, pointsize = 12, ht = 6, wd = 6)
  stripchart(data.list[, 2:(samples+1)], xlab = "", col = "blue", jitter = 0.2, 
          method = "jitter",  main = Title, vertical = T,
          pch = 1, cex = 0.8, ylab = "Test result", xlim = c(0.5, samples+0.5))
  # mtext(Title, side = 3, line = 1.5, cex = 1.2)
  mtext(test.name, side = 3, line = 0.5, cex = 1)
  CloseGraphOutput("R")

# output results
d1<- paste(substr(date(), 1, 10), substr(date(), 20,24), " @", substr(date(), 11, 16))
# reformat headers
output<- paste("
Repeatability analysis for", test.name, "
\n")
output<- paste(output, "
Results
\n")

# start building HTML output
output<- paste(output, "", "Analysed: ", d1, "\n", sep="")
output<- paste(output, "Summary of testing data\n", sep="")
# summary results
tmp1<- HTMLStream(t(summary.results), cellborder = 0, classfirstline = "mbg", classfirstcolumn = "mbg", classcellinside = "left_mar",  cellalign = "center", align = "left", digits = digits)
    graphs<- array("", dim = c(1, 2))
    colnames(graphs)<- c("Box plots", "Stripcharts")
    graph<- paste(fpath, "tmp/", filename, "_boxplot.png", sep="")    
    graphs[1,1]<- paste("\n", sep="")
    graph<- paste(fpath, "tmp/", filename, "_stripchart.png", sep="")    
    graphs[1,2]<- paste("\n", sep="")
    tmp2<- HTMLStream(graphs, cellborder = 0, classfirstline = "mbg", classfirstcolumn = "mbg", classcellinside = "left_mar",  cellalign = "center", align = "left", digits = digits)
output<- paste(output, tmp1, tmp2, sep="")
# cat("
", test)

# statistical control data and charts for each sample
output<- paste(output, "Summary of performance for each sample\n", sep="")
for (s in 1:samples) {
    graphs<- array("", dim = c(1, 2))
    colnames(graphs)<- c("Time series plot", "Autocorrelation plot")
    graph<- paste(fpath, "tmp/", filename, "_scc_", colnames(data.list)[s+1], ".png", sep="")    
    graphs[1,1]<- paste("\n", sep="")
#    graph<- paste("tmp/", filename, "_LR_summary_", colnames(data.list)[s+1], ".png", sep = "")
#    graphs[1,2]<- paste("\n", sep="")
    graph<- paste(fpath, "tmp/", filename, "_acf_", colnames(data.list)[s+1], ".png", sep="")    
    graphs[1,2]<- paste("\n", sep="")
# cat("
", test)

    tmp1<- paste("", colnames(data.list)[s+1], "
\nNumbers and proportions of testing occasions exceeding calculated limits
\n", sep="")
    tmp2<- HTMLStream(limits[,,s], cellborder = 0, classfirstline = "mbg", classfirstcolumn = "mbg", classcellinside = "left_mar",  cellalign = "center", align = "left", digits = digits)
    tmp3<- HTMLStream(graphs, cellborder = 0, classfirstline = "mbg", classfirstcolumn = "mbg", classcellinside = "left_mar",  cellalign = "center", align = "left", digits = digits)
    tmp4<- "
Summary linear regression model
" 
    for (r in 1:length(lm.summary[[s]])) {
        tmp4<- paste(tmp4, lm.summary[[s]][r], "
")
    }
        tmp4<- paste(tmp4, "")

    output<- paste(output, tmp1, tmp2, tmp3, tmp4, sep="")
} # end of s loop


sink()
# write output to file
file.name<- file(paste(tmp.file, "_result.html", sep = ""), open = "wt")
cat(output, file=file.name)  
close(file.name)           #)
cat(output)
# 

# calculate cusum
cusum<- array(0, dim = c(tests+1, samples))
target<- limits[3, 1,]
cusum[1,]<- 0
for (r in 1:tests) {
    vals<- ifelse(is.na(data.list[r, 2:(samples+1)]), target, data.list[r, 2:(samples+1)])
    cusum[r+1,]<- cusum[r,] + (vals - target)
} # end of r loop

# plot(cusum, type = "l", col = "darkblue", x = 0:tests, xaxt = "n")
# axis(side = 1, labels = c(0, data.list[,1]), at = 0:tests)
# mtext(colnames(data.list)[s+1], side = 3, line = 0.5, cex = 0.8)
# abline(a = 0, b = 0, col = "blue")