An official website of the United States government.

This is not the current EPA website. To navigate to the current EPA website, please go to www.epa.gov. This website is historical material reflecting the EPA website as it existed on January 19, 2021. This website is no longer updated and links to external websites and some internal pages may not work. More information »

# Using R to Classify Curve Shapes

## How to Classify Curve Shapes

Curve shapes can be classified as increasing, decreasing, or unimodal, by comparing the mean responses to the confidence intervals (see the page on Placing Taxa into Tolerance Categories in the Helpful Links box).

The R script provided below classifies taxon-environment relationships by curve shape. Note that unimod.test is an R function, a series of commands that is executed when unimod.test is called.

unimod.test <- function(mnr, ubnd, lbnd) {

# Find the maximum and minimum predicted mean probabilities
lmax <- max(mnr)
lmin <- min(mnr)

# Find index locations for these probabilities
imax <- match(lmax, mnr)
imin <- match(lmin, mnr)

x.out <- F
y.out <- F

# Compare mean predicted probability to the left of maximum point
# with upper confidence bound.  Store a T in x.out if
# any point in the mean response deviates from the
# upper confidence limit
if (imax > 1) {
x.out <- sum(lmax == pmax(lmax, ubnd[1:(imax-1)])) > 0
}

# Store a T in y.out if any point in the mean probability
# to the right of the maximum point deviates from the upper
# confidence limit
if (imax < length(ubnd)) {
y.out <- sum(lmax == pmax(lmax,
ubnd[(imax+1):length(ubnd)])) > 0
}

# Perform same set of tests for lower confidence limit
a.out <- F
b.out <- F

if (imin > 1) {
a.out <- sum(lmin == pmin(lmin, lbnd[1:(imin-1)])) > 0
}
if (imin > length(lbnd)) {
b.out <- sum(lmin == pmin(lmin,
lbnd[(imin+1):length(lbnd)])) > 0
}

# The information on where the mean curve deviates from the
# confidence limits tells us its curve shape...
if (x.out & y.out) {
return("Unimodal")
}
if (a.out & b.out) {
return("Concave up")
}
if (x.out | b.out) {
return("Increasing")
}
if (y.out | a.out) {
return("Decreasing")
}
if (! (x.out | y.out | a.out | b.out)) return(NA)
}

tolcl <- rep("", times = length(taxa.names))
for (i in 1:length(taxa.names)) {
predres <- predict(modlist.gam, type= "link", se.fit = T)

# Compute upper and lower 90% confidence limits

# Convert from logit transformed values to probability.