An R program for constructing pointwise and simultaneous functional tolerance bands using the methodology proposed in a Biometrics article titled "Tolerance bands for functional data" by L.N Rathnayake and P. K. Choudhary.

# Download the program:

Right click on the link TB_functions.R and save the associated file as TB_functions.R (an R file) in the R working directory. The file contains all the functions that are needed to compute tolerance bands. The computation is actually done by the TB.calc function. See the example below for how to run the program.

# Requires:

The program requires "mgcv" and "refund" packages. They can be obtained from the R website.

# Usage:

TB.calc(y, alpha = 0.05, prob = 0.90, simul = TRUE, n.boot)

# Arguments:

y = functional data that will be used to construct tolerance bands. It must have the same format as "Y" in the ccb.fpc function in the refund package.

alpha = one minus confidence level so that 1-alpha is the  confidence level (default = 0.05)

prob = desired probability content for the tolerance band (default = 0.90)

simul = logical argument to indicate whether a simultaneous band is desired in addition to the pointwise band(default = TRUE)

n.boot = number of bootstrap repetitions to compute the tolerance factor (should be at least 500 to get a reasonable estimate)

# Values:

A list containing values of the desired tolerance band(s) and the estimated model parameters.

# Author(s):

L.N Rathnayake and P.K. Choudhary,[email protected]), University of Texas at Dallas

# Reference(s):

Rathnayake, L.N, Choudhary, P.K.,(2015). Tolerance bands for functional data. Submitted.

# Illustration using CD4 data from refund package:

Step 1. Launch R.

Step 2. Change the working directory to the directory where TB_functions.R file is stored.

Step 3. Run source("TB_functions.R") to read the functions in the file. This will also load the mgcv and refund packages.

Step 4. Run the following commands in R (takes about 30 minutes on a windows computer with 2.4 GHz processor and 2 GB RAM)

data(cd4)

y <-  cd4

result <- TB.calc(y, alpha = .05, prob = 0.90, n.boot = 500, simul = T)

 

Step 5. Run the following commands in R to plot the tolerance bands.

 

# two-sided band

 

plot(y[1,], ylim = c(-114, 3500), xaxt = "n", xlab = "months since seroconversion", ylab = "cd4 count")

axis(1, at = seq(0, 65, 10), labels = seq(-20, 45, 10))

 

for(j in 1:dim(y)[1]){

    lines(x = which(y[j,] != 'NA'), y = y[j,][which(y[j,] != 'NA')], type = 'o', col = "gray80")

  }

 

  lines(result[[3]][,1], type = 'l', lty = 2, col = 'red', lwd = 2)

  lines(result[[3]][,2], type = 'l', lty = 2, col = 'red', lwd = 2)

  lines(result[[4]][,1], type = 'l', lty = 1, col = 'blue',lwd = 2)

  lines(result[[4]][,2], type = 'l', lty = 1, col = 'blue',lwd = 2)

  lines(result[[5]][,1], type = 'l', lty = 3, col = 'black',lwd = 2)

 

  legend("topright", legend = c("estimated mean", "pointwise band", "simultaneous band"), lty = c(3, 2, 1), lwd = c(2, 2, 2), col = c("black", "red", "blue"), cex = .75)

 

# output

# one-sided upper band

 

plot(y[1,], ylim = c(0,3500), xaxt = "n", xlab = "months since seroconversion", ylab = "cd4 count")

axis(1, at = seq(0,65,10), labels = seq(-20,45,10))

 

for(j in 1:dim(y)[1]){

  lines(x = which(y[j,] != 'NA'), y = y[j,][which(y[j,] != 'NA')], type = 'o', col = "gray80")

}

 

lines(result[[1]][,1], type = 'l', lty = 2, col = 'red', lwd = 2)

lines(result[[2]][,1], type = 'l', lty = 1, col = 'blue',lwd = 2)

lines(result[[5]][,1], type = 'l', lty = 3, col = 'black',lwd = 2)

legend("topright", legend = c("estimated mean", "pointwise band", "simultaneous band"), lty = c(3, 2, 1), lwd = c(2, 2, 2), col = c("black", "red", "blue"), cex = .75)

 

# output