# Quantitative Ecotoxicology, Page 178, Example 4.9, Time-to-death

This is example 4.9 on page 178 of Quantitative Ecotoxicology - time-to-death data.

Thankfully, Prof. Newman provided me the data for this example. You can get it from the github-repo (TOXICTY.csv).

The data consists of 5 columns:

• TTD : Time to death
• TANK : Tank
• PPT : NaCl Concentration
• WETWT : wet weight
• STDLGTH : Standard length

Columns 4 and 5 have 70 NA’s (no data available due to measurement error), but we won’t use these in this example. The observations with TTD = 97 are ‘survivors’, since the experiment run only 96 hours.

First we need to create a column FLAG for the status of the animal (dead/alive):

So 1 denotes alive and 2 dead.

Then we can plot the data. Each line is a tank and colors denote the NaCl concentrations.

We see a clear relationship between concentration and the survival curves. In this example we are interested in differences between the duplicates. We see that the two curves for the 11.6 g/L concentration are quite similar, while there is more divergence between tanks in the 13.2 g/L treatment.

We can test for differences using the survdiff function. With the rho argument we can specify the type of test: rho = 0 is a log-rank test and rho = 1 is equivalent to the Peto & Peto modification of the Gehan-Wilcoxon test.

First the log-rank test for each concentration:

We could also run this in a for loop (here the Wilcoxon test):

Basically we get the same results as in the book:

None of the log-rank tests is statistically significant (at the 0.05 level). The wilcoxon test for the 13.2 g/L treatment shows a p < 0.05. This is also in agreement with the plot.

The $\chi^2$ values differ slightly but share the same trend - I suspect this is due to different data used.

With this dataset we can do much more. We already saw that there might be a relationship between survival time and concentration, but more on this later (example 4.10).

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

# Quantitative Ecotoxicology, Page 162, Example 4.7, Duplicate Treatments

This is a short one (example 4.7 on page 1621 of Quantitative Ecotoxicology).

First we create the data as matrix:

The we can easily run fisher’s Exact test for this 2x2 table using the function fisher.test():

The results are identical to the one in the book.

# Quantitative Ecotoxicology, Page 159, Example 4.6, Control Mortality

This is example 4.6 on page 159 of Quantitative Ecotoxicology. It is about how to deal with control mortalities.

First we need the data:

The data consists of number of dead animals (DEAD) from all animals (TOTAL) exposed to different concentrations (CONC). First we create a new column with the proportion of dead animals:

Here is a plot of the data. Note the use of expression() (greek letters in the axis labels).

### Control mortality

We can estimate the mean control mortality and the confidence interval for the mean using the t.test function:

These can be also easily extracted from the t.test object:

This gives nearly the same values as in the book. I don’t know what the SAS option OPTC is doing or computing, however it seems like it is the mean +- CI for the control group.

### Abbott’s formula

We could adjust for control mortality using Abbott’s formula:

with $p_c$ = the corrected, p = original and $p_0$ = control mortality.

The mean control mortality can be calculated as:

And the corrected mortalities using Abbotts formula as:

### Dose-Response-Models

#### Ignoring control mortality

As in the previous example we can fit a dose-response-model to this data using the drc package:

Comparing with other model this models performs quite good.

#### Using the corrected mortalities

We can also fit a model to the corrected mortalities PROP_c.

Abbotts correction resulted to some negative mortalities, therefore I set the control and all negative mortalities to zero:

Then we fit a dose-response model:

However a Weibull model fits slightly better the data, so I change to a two-parameter Weibull model (using the update function).

#### A model without fixed lower limit

The two-parameter log-logistic model from above (mod1) performs quite good. However its lower limit is fixed to 0 and the upper limit to 1. Since we have a small amount of control mortality we could check if a model with varying lower limit (will be estimated) makes sense.

Let’s fit a three parameter log-logistic function, where the lower limit is an additional parameter:

However looking at the summary we see that the lower limit is estimated as $0.007 \pm 0.02$ and is statistically not significant.

Since the lower limit (=control mortality) is so low we could also stick with mod1 since it brings no improvement.

All three considered models give nearly the same $LC_{50}$ around 2100:

Code and data is also available on my github-repo under file name ‘p159’.

# R, Knitr and Octopress

Recently I moved this blog from tumblr to github since I feeled I need more liberty. I have no experience with web programming, so I chose Octopress as blogging framework since it nicely integrates with github.

However after relocation my workflow for generating posts was quite laborious - it was similar to the tumblr-workflow:

• Create an .Rmd file
• knit the file
• move the file to octopress/source/_posts/
• move the figures to octopress/source/figure/
• and so on…

This was non-satisfying and I thought that this must be automated into a smooth workflow!

Last week I took a look at the new Rcpp Gallery and noticed that this is also hosted on github. I skimmed through their code and found a lot of useful stuff. Shamelessly I copied some files from the and modified them for my octopress needs.

Basically I have now a folder for my .Rmd files (octopress/source/src). Invoking make in this folder creates the markdown-files in _posts and sets up the right figure paths (in _posts/figure). My workflow now:

• rake new_rmd[test], this sets up an empty .Rmd file in /src, according to the octopress naming-scheme and header.
• add my post to this .Rmd file
• invoke make in /src, this sets up the .markdown file in _posts and figures
• the usual rake gen_deploy

I am quite happy at the moment with this workflow. I can also use R-Studio.

Here are the changes I made to my octopress:

1) I added a new function rake new_rmd to my Rakefile. This mimics the behavior of rake new_post and creates a .Rmd-file in source/src. Simply add these lines to your Rakefile:

2) I modified and cleaned knit.sh from the Rcpp Gallery for my needs. This runs knitr on the .Rmd files and saves the output to _posts. Put this file into source/_scripts!

3) Modified and cleaned the Makefile from the Rcpp Gallery for my needs. This runs knit.sh on the /src folder, returns the .markdown files to _posts/ and cleans up the .md and .hmtl files created by R-Studio. Put this file into source/src!

You can skim through my github-repo for this site and copy the files from there. I have not tested this very extensively - in case you find any bugs or have improvements, please let me know.

Edit 01.03.2013

 to source/_includes/custom/head.html so that the figures are displayed properly, see here. 
 
 Quantitative Ecotoxicology, Page 147, Example 4.3, LC50 Feb 25th, 2013 This is about example 4.3 on page 147 of Quantitative Ecotoxicology. This example is about Calculations of $LC_{50}$ values. In this post I won’t reproduce the SAS-Code since I do not have any experience with SAS PROC PROBIT and I do not fully understand whats happening there. Instead I will fit Dose-Response-Models using the drc-package to this data. Get the data from here and read it into R: 1 2 3 4 require(RCurl) url <- getURL("https://raw.github.com/EDiLD/r-ed/master/quantitative_ecotoxicology/data/p147.csv", ssl.verifypeer = FALSE) salt <- read.table(text = url, header = TRUE, sep = ";") 1 head(salt) 1 2 3 4 5 6 7 ## DEAD TOTAL CONC ## 1 16 76 10.3 ## 2 22 79 10.8 ## 3 40 77 11.6 ## 4 69 76 13.2 ## 5 78 78 15.8 ## 6 77 77 20.1 So the data consists of number of dead animals (DEAD) from all animals (TOTAL) exposed to a concentration (CONC). First we create a new column with the proportion of dead animals: 1 salt$prop <- salt$DEAD/salt$TOTAL Lets have a look at the raw data (note that I use a logarithmic scale for the x-axis): 1 2 plot(salt$CONC, salt$prop, xlab = "Concentration", ylab = "Proportion dead", log = "x") I will use the drc-package of Christian Ritz and Jens Strebig to fit dose-response-curves to this data. The main function of this package is drm: Here I fit a two-parameter log-logistic model to the data (see Ritz (2010) for a review of dose-response-models used in ecotoxicology): 1 2 require(drc) mod <- drm(prop ~ CONC, data = salt, fct = LL.2()) So the usage is similar to lm() or nls(), except the fct argument. This argument defines the model that is fitted to the data. We can compare this model with other models using the AIC (the smaller the better). Here I compare the 2-parameter log-logistic model with a two-parameter Weibull and a 2-parameter Gompertz model. drc has for this purpose the mselect() function that shows some diagnostics for every model: likelihood (the higher the better; however use this only when your models are nested) AIC (the lower the better) residual variance (the lower the better) 1 mselect(mod, fctList = list(W1.2(), G.2())) 1 2 3 4 ## logLik IC Lack of fit Res var ## LL.2 14.613 -23.226 NA 0.00067322 ## W1.2 NA NA NA NA ## G.2 NA NA NA NA The LL.2-model has the lowest AIC so I will keep this. Lets see how the model looks like: 1 2 3 4 5 6 # raw data plot(prop ~ CONC, data = salt, xlim = c(9, 21), ylim = c(0, 1), ylab = "Proportion dead", xlab = "Concentration") conc_pred <- seq(9, 21, 0.1) lines(conc_pred, predict(mod, newdata = data.frame(CONC = conc_pred))) We can get the$LC_{50}$with confidence interval from the model using the ED() function: 1 ED(mod, 50, interval = "delta") 1 2 3 4 5 6 ## ## Estimated effective doses ## (Delta method-based confidence interval(s)) ## ## Estimate Std. Error Lower Upper ## 1:50 11.4768 0.0575 11.3171 11.6 Code and data is available on my github-repo under file name ‘p147’. References Ritz C (2010). “Toward A Unified Approach to Dose-Response Modeling in Ecotoxicology.” Environmental Toxicology And Chemistry, 29. ISSN 07307268, http://dx.doi.org/10.1002/etc.7. ← Older Blog Archives   Tag Cloud Recent Posts Quantitative Ecotoxicology, Page 178, Example 4.9, Time-to-death Quantitative Ecotoxicology, Page 162, Example 4.7, Duplicate Treatments Quantitative Ecotoxicology, Page 159, Example 4.6, Control Mortality R, knitr and Octopress Quantitative Ecotoxicology, Page 147, Example 4.3, LC50 GitHub Repos Status updating… @EDiLD on GitHub$.domReady(function(){ if (!window.jXHR){ var jxhr = document.createElement('script'); jxhr.type = 'text/javascript'; jxhr.src = '/javascripts/libs/jXHR.js'; var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(jxhr, s); } github.showRepos({ user: 'EDiLD', count: 0, skip_forks: false, target: '#gh_repos' }); }); 
 
 Copyright © 2013 - Eduard Szöcs - Powered by Octopress var disqus_shortname = 'redgh'; var disqus_script = 'count.js'; (function () { var dsq = document.createElement('script'); dsq.type = 'text/javascript'; dsq.async = true; dsq.src = 'http://' + disqus_shortname + '.disqus.com/' + disqus_script; (document.getElementsByTagName('head')[0] || document.getElementsByTagName('body')[0]).appendChild(dsq); }()); MathJax.Hub.Config({ jax: ["input/TeX", "output/HTML-CSS"], tex2jax: { inlineMath: [ ['$', '$'] ], displayMath: [ ['$$', '$$']], processEscapes: true, skipTags: ['script', 'noscript', 'style', 'textarea', 'pre', 'code'] }, messageStyle: "none", "HTML-CSS": { preferredFont: "TeX", availableFonts: ["STIX","TeX"] } });