Eduard Szöcs

Data in Environmental Science and Eco(toxico-)logy

Quantitative Ecotoxicology, Page 85, Example 3.2, Nonlinear Regression

Get the data from here and read it into R:

OYSTERZN <- read.table("p85.csv", header = TRUE, sep = ";")
head(OYSTERZN)
##   DAY ZINC
## 1   1  700
## 2   1  695
## 3   1  675
## 4   1  630
## 5   1  606
## 6   1  540
plot(ZINC ~ DAY, data = OYSTERZN)

plot of chunk p85_raw

First we fit a nonlinear Regression without weighting.

mod_noweight <- nls(ZINC ~ INITACT * exp((-(KE + 0.00283)) * DAY), data = OYSTERZN, 
    start = list(KE = 0.004, INITACT = 500))

So we fit the model

to our data.

In the R formula ZINC corresponds to $C_t$,
INITACT to ,
KE to ,
0.00283 is the decay rate constant for 65-Zn $k_{e2}$,
and DAY to $t$.

Were are going to estimate KE and INITACT and also supplied some start-values for the algorithm.

We can look a the summary to get the estimates and standard error:

summary(mod_noweight)
## 
## Formula: ZINC ~ INITACT * exp((-(KE + 0.00283)) * DAY)
## 
## Parameters:
##         Estimate Std. Error t value Pr(>|t|)    
## KE      2.68e-03   6.68e-04    4.01  0.00015 ***
## INITACT 4.65e+02   2.04e+01   22.86  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Residual standard error: 98.8 on 71 degrees of freedom
## 
## Number of iterations to convergence: 3 
## Achieved convergence tolerance: 5.57e-06

The resulting estimates of and are and .

We can investigate the residuals, which show a clear pattern:

res_mod_noweight <- resid(mod_noweight)
plot(OYSTERZN$DAY, res_mod_noweight)
abline(h = 0)

plot of chunk p85_residuals_nls

Secondly, we run a nonlinear regression with day-squared weighting:

We use day^2 as weights and add there a column to our data:

OYSTERZN$WNLIN <- OYSTERZN$DAY^2

We run again nls, but now we supply this new column as weights:

mod_weight <- nls(ZINC ~ INITACT * exp((-(KE + 0.00283)) * DAY), data = OYSTERZN, 
    weights = OYSTERZN$WNLIN, start = list(KE = 0.004, INITACT = 500))
summary(mod_weight)
## 
## Formula: ZINC ~ INITACT * exp((-(KE + 0.00283)) * DAY)
## 
## Parameters:
##         Estimate Std. Error t value Pr(>|t|)    
## KE      2.44e-03   3.47e-04    7.03  1.1e-09 ***
## INITACT 4.55e+02   3.82e+01   11.89  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Residual standard error: 7550 on 71 degrees of freedom
## 
## Number of iterations to convergence: 4 
## Achieved convergence tolerance: 5.42e-07

The estimates and are quite similar to the non-weighted regression.

We could plot the two models and the data:

# extract the fitted values
fit_mod_noweight <- fitted(mod_noweight)
fit_mod_weight <- fitted(mod_weight)
# plot data
plot(ZINC ~ DAY, data = OYSTERZN)
# add fitted values
lines(OYSTERZN$DAY, fit_mod_noweight)
lines(OYSTERZN$DAY, fit_mod_weight, lty = "dashed")
# add legend
legend("topright", legend = c("nonweighted", "weighted"), lty = c("solid", "dashed"))

plot of chunk p85_nls_fitted

Finally we can also fit a linear model to the transformed Zinc-Concentrations:

First we ln-transform the concentrations:

OYSTERZN$LZINC <- log(OYSTERZN$ZINC)

We see that the data has now linear trend:

plot(LZINC ~ DAY, OYSTERZN)

plot of chunk p85_raw2

And fit a linear regression:

mod_lm <- lm(LZINC ~ DAY, data = OYSTERZN)
sum_lm <- summary(mod_lm)
sum_lm
## 
## Call:
## lm(formula = LZINC ~ DAY, data = OYSTERZN)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8974 -0.2448 -0.0709  0.2958  0.5715 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.073833   0.056834   106.9   <2e-16 ***
## DAY         -0.005314   0.000243   -21.9   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Residual standard error: 0.354 on 71 degrees of freedom
## Multiple R-squared: 0.871,  Adjusted R-squared: 0.869 
## F-statistic:  478 on 1 and 71 DF,  p-value: <2e-16

which is fitting the model , with a = -0.0053 and intercept = 6.07

Now plot data and model, as well as the residuals:

# fitted values
fit_mod_lm <- fitted(mod_lm)

# data + fitted
plot(LZINC ~ DAY, OYSTERZN)
lines(OYSTERZN$DAY, fit_mod_lm)

plot of chunk p85_lm_fitted

# residuals
plot(mod_lm, which = 1)

plot of chunk p85_lm_fitted

The mean square error can be calculated from the summary:

# MSE
sum_lm$sigma^2
## [1] 0.12516

From which we can get an unbiased estimate of :

# unbiased estimate for C_0: exp(MSE/2) * exp(Intercept)
exp(sum_lm$sigma^2/2) * exp(sum_lm$coefficients[1, 1])
## [1] 462.39

where

sum_lm$coefficients[1, 1]

extracts the intercept from the summary.

The estimated in the summary output is , and

This result is similar to the weighted and non weighted nonlinear regression. Again we have the same results as with SAS :) [Small deviations may be due to rounding error]

Code and data are available at my github-repo under file name ‘p85’.

Written on December 1, 2012