Eduard Szöcs

Data in Environmental Science and Eco(toxico-)logy

Quantitative Ecotoxicology, Page 109, Example 3.8, Bioaccumulation

This is example 3.8 on page 109 of Quantitative Ecotoxicology - reproduced with R. This example is about accumulation and elimination of bromophos from water in a guppy (Poecilia reticulata).

There are two data files for this example - one for the accumulation and on for the elimination.

Accumulation

First we will look at the accumulation phase:

require(RCurl)
# Accumulation
url_accum <- getURL("https://raw.github.com/EDiLD/r-ed/master/quantitative_ecotoxicology/data/p109_accum.csv",
ssl.verifypeer = FALSE, .opts=curlOptions(followlocation=TRUE))
ACCUM <- read.table(text = url_accum, header = TRUE, sep = ";")
head(ACCUM)
##   HOUR BRPHOS
## 1  0.5   1900
## 2  1.0   3000
## 3  2.0   5200
## 4  4.0   6900
## 5  8.0  24000
## 6 24.0  50000

Again we have two columns: One for the time and one for the concentration.

We fit can same model as in example 3.7 to this data. The uptake $(k_u)$ and elimination $(k_e)$ constants are estimated simultaneously (at the same time):

mod_accum <- nls(BRPHOS ~ KU / KE * 10.5 * (1 - exp(-KE * HOUR)),
           data = ACCUM, 
           start = list(KU = 100, KE = 0.01))

Note that I used different starting values than in the SAS-Code (must be a typo in the book). Also I didn’t specify any bounds.

summary(mod_accum)
## 
## Formula: BRPHOS ~ KU/KE * 10.5 * (1 - exp(-KE * HOUR))
## 
## Parameters:
##     Estimate Std. Error t value Pr(>|t|)    
## KU 344.79785   31.85529   10.82  4.7e-06 ***
## KE   0.00525    0.00103    5.09  0.00094 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18900 on 8 degrees of freedom
## 
## Number of iterations to convergence: 7 
## Achieved convergence tolerance: 3.72e-06
HOUR_pred <- seq(min(ACCUM$HOUR), max(ACCUM$HOUR), by = 0.1) 
# Raw data
plot(ACCUM, main = 'Accumulation')
# add model
lines(HOUR_pred, predict(mod_accum, newdata = data.frame(HOUR = HOUR_pred)))

plot of chunk plot_accum_model

So from the accumulation data we estimated the uptake and elimination constants as:

  • $k_e = 0.0053 \pm 0.0010$
  • $k_u = 344.798 \pm 31.855$

Sequential estimation

However we could also estimate the elimination constant $(k_e)$ from the elimination phase and then use this estimate for our accumulation data.

  • First estimate $k_e$ from a linear model (linear transformation)
  • Plug this estimated $k_e$ into a nonlinear model to estimate $k_u$
# Elimination data
url_elimin <- getURL("https://raw.github.com/EDiLD/r-ed/master/quantitative_ecotoxicology/data/p109_elimin.csv",
ssl.verifypeer = FALSE)
ELIMIN <- read.table(text = url_elimin, header = TRUE, sep = ";")
## Error in read.table(text = url_elimin, header = TRUE, sep = ";"): no lines available in input
head(ELIMIN)
## Error in head(ELIMIN): error in evaluating the argument 'x' in selecting a method for function 'head': Error: object 'ELIMIN' not found
plot(ELIMIN)
## Error in plot(ELIMIN): error in evaluating the argument 'x' in selecting a method for function 'plot': Error: object 'ELIMIN' not found

We will estimate $k_e$ from a linear model like in previous examples. We could also use nls for this.

First we need to transform the bromophos-concentration to linearize the relationship.

ELIMIN$LBROMO <- log(ELIMIN$BRPHOS)
## Error in eval(expr, envir, enclos): object 'ELIMIN' not found

The we can use lm() to fit the linear model:

mod_elimin_lm <- lm(LBROMO ~ HOUR, data = ELIMIN)
## Error in is.data.frame(data): object 'ELIMIN' not found
summary(mod_elimin_lm)
## Error in summary(mod_elimin_lm): error in evaluating the argument 'object' in selecting a method for function 'summary': Error: object 'mod_elimin_lm' not found

So we get an estimate of $k_e$ as $0.0147 \pm 0.0003$.

This is quite different to the $k_e$ estimated simultaneous from the accumulation data! Our linear model fits very good (R^2 = 0.998, no pattern in the residuals), so something is strange here…

par(mfrow = c(1, 2))
# plot linearized model
plot(LBROMO ~ HOUR, data = ELIMIN, main = 'Data + Model')
## Error in eval(expr, envir, enclos): object 'ELIMIN' not found
# add regression line
abline(mod_elimin_lm)
## Error in abline(mod_elimin_lm): object 'mod_elimin_lm' not found
# plot residuals
plot(fitted(mod_elimin_lm), residuals(mod_elimin_lm), main = 'Residuals')
## Error in plot(fitted(mod_elimin_lm), residuals(mod_elimin_lm), main = "Residuals"): error in evaluating the argument 'x' in selecting a method for function 'plot': Error in fitted(mod_elimin_lm) : object 'mod_elimin_lm' not found
abline(h = 0, lty = 'dotted')
## Error in int_abline(a = a, b = b, h = h, v = v, untf = untf, ...): plot.new has not been called yet

Plug $k_e$ from the elimination phase into the accumulation model

Lets take $k_e$ from the elimination phase and plug it into our accumulation model and investigate the differences:

mod_accum2 <- nls(BRPHOS ~ KU / -coef(mod_elimin_lm)[2] * 10.5 * (1 - exp(coef(mod_elimin_lm)[2] * HOUR)),
           data = ACCUM, 
           start = list(KU = 100))
## Error in nls(BRPHOS ~ KU/-coef(mod_elimin_lm)[2] * 10.5 * (1 - exp(coef(mod_elimin_lm)[2] * : parameters without starting value in 'data': mod_elimin_lm
summary(mod_accum2)
## Error in summary(mod_accum2): error in evaluating the argument 'object' in selecting a method for function 'summary': Error: object 'mod_accum2' not found

This estimates $k_u = 643.9 \pm 40.4$ which differs greatly from our initial results! Lets plot this model and the residuals:

par(mfrow=c(1,2))
HOUR_pred <- seq(min(ACCUM$HOUR), max(ACCUM$HOUR), by = 0.1) 
# Raw data
plot(ACCUM, main = 'Accumulation')
# add model
lines(HOUR_pred, predict(mod_accum2, newdata = data.frame(HOUR = HOUR_pred)))
## Error in predict(mod_accum2, newdata = data.frame(HOUR = HOUR_pred)): error in evaluating the argument 'object' in selecting a method for function 'predict': Error: object 'mod_accum2' not found
plot(fitted(mod_accum2), residuals(mod_accum2))
## Error in plot(fitted(mod_accum2), residuals(mod_accum2)): error in evaluating the argument 'x' in selecting a method for function 'plot': Error in fitted(mod_accum2) : object 'mod_accum2' not found

plot of chunk plot_accum_model2

The residuals show a clear curve pattern. But we could also look at the residual sum of squares and the AIC to see which model fit better to the accumulation data:

# Residual sum of squares
mod_accum$m$deviance()
## [1] 2870355583
mod_accum2$m$deviance()
## Error in eval(expr, envir, enclos): object 'mod_accum2' not found
# AIC
AIC(mod_accum)
## [1] 229.13
AIC(mod_accum2)
## Error in AIC(mod_accum2): object 'mod_accum2' not found

So the first model seem to better fit to the data. However see the discussion in the book for this example!

Once again we reproduced the results as in the book using R :)

Code and data are available on my github-repo under file name ‘p109’.

Written on February 24, 2013