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 »

CADDIS Volume 4

Using R to Estimate Environmental Limits

How to Estimate Environmental Limits

We can use a for loop similar to that used for weighted averages to compute environmental limits (see Environmental Limits page, Equation 2) for different taxa.

In this method, you must select a value of the cumulative percentile that you wish to use for computing the environmental limit. In the script below, the cutoff value is selected to be 0.75.

First make sure that you have loaded the sample biological and environmental data and merged them into a single data frame called dfmerge (see Download Scripts and Sample Data in the Helpful Links box).

Also make sure that you have selected the taxa for which you wish to calculate environmental limits and saved them in the vector taxa.names (as described in the R Script for Central Tendencies  in the Helpful Links box).

# Define a storage vector for the cumulative percentile

CP <- rep(NA, times = length(taxa.names))

# Sort sites by the value of the environmental variable
dftemp <- dfmerge[order(dfmerge$temp), ]

# Select a cutoff percentile
cutoff <- 0.75

# Specify three plots per page
par(mfrow = c(1,3), pty = "s")
for (i in 1:length(taxa.names)) {
  # Compute cumulative sum of abundances
  csum <- cumsum(dftemp[, taxa.names[i]])/
          sum(dftemp[,taxa.names[i]])

  # Make plots like Figure 4
  plot(dftemp$temp, csum, type = "l", xlab = "Temperature",
       ylab = "Proportion of total", main = taxa.names[i])

  # Search for point at which cumulative sum is 0.75
  ic <- 1
  while (csum[ic] < 0.75) ic <- ic + 1

  # Save the temperature that corresponds to this percentile.
  CP[i] <- dftemp$temp[ic]

}
names(CP) <- taxa.names
print(CP)

Top of Page