Q & K Box 13.3

Calculations in ecotoxicology. The data used here are summary vaules taken from the Bliss data used by Finney in his description of probit analysis. Data are in 2 forms - individual live/dead values, and summary proportions for each dose.

The data files are here and here

Preliminaries

First, load the required packages (vcd)

source("../R/libraries.R")   #This is the common library
library(performance)
library(vcd)

Disaggregated data file

finndisagg <- read_csv("../data/finneydisagg.csv")

Rows: 292 Columns: 4── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
dbl (4): dose, dead, indiv, ldose
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(finndisagg)

plot(dead~dose, data=finndisagg)

Fit GLM with binomial response variable

fdisagg.lm<-glm(dead~dose,family=binomial, data=finndisagg)
tidy(fdisagg.lm, conf.int=TRUE)

plot(fdisagg.lm$fitted.values~finndisagg$dose)

coef(fdisagg.lm)

(Intercept)        dose 
 -3.2798094   0.6117899 

summary(fdisagg.lm) #get model fit

Call:
glm(formula = dead ~ dose, family = binomial, data = finndisagg)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.4539  -0.6157  -0.2718   0.6548   1.8742  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.27981    0.37646  -8.712   <2e-16 ***
dose         0.61179    0.06897   8.870   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 402.11  on 291  degrees of freedom
Residual deviance: 249.19  on 290  degrees of freedom
AIC: 253.19

Number of Fisher Scoring iterations: 5

Calc LC50 Use parameter values from the dose-response, solve for 50%

3.28/.612

[1] 5.359477

Original analysis used log 10 transform of dose

fdisagg.lm3<-glm(dead~ldose,family=binomial, data=finndisagg, na.action=na.omit)
tidy(fdisagg.lm3, conf.int=TRUE)

summary(fdisagg.lm3) #get model fit

Call:
glm(formula = dead ~ ldose, family = binomial, data = finndisagg, 
    na.action = na.omit)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.2191  -0.8573   0.4221   0.6445   1.9749  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -5.1177     0.6648  -7.698 1.38e-14 ***
ldose         7.4269     0.9248   8.031 9.70e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 335.05  on 242  degrees of freedom
Residual deviance: 240.14  on 241  degrees of freedom
  (49 observations deleted due to missingness)
AIC: 244.14

Number of Fisher Scoring iterations: 4

solve for y=0; use coefficients to calculate LC50

coef(fdisagg.lm3)

(Intercept)       ldose 
  -5.117703    7.426909 

5.12/7.43

[1] 0.6890983

backtransform to get LC50 on raw dose scale

10^0.689

[1] 4.886524

Fit same model as probit, rather than logit

fdisagg.lm4<-glm(dead~ldose,family=binomial(link="probit"), data=finndisagg, na.action=na.omit)
tidy(fdisagg.lm4, conf.int=TRUE)

summary(fdisagg.lm4) #get model fit

Call:
glm(formula = dead ~ ldose, family = binomial(link = "probit"), 
    data = finndisagg, na.action = na.omit)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.2347  -0.8651   0.4145   0.6614   1.9755  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -3.0308     0.3648  -8.308   <2e-16 ***
ldose         4.3827     0.4994   8.775   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 335.05  on 242  degrees of freedom
Residual deviance: 240.48  on 241  degrees of freedom
  (49 observations deleted due to missingness)
AIC: 244.48

Number of Fisher Scoring iterations: 4

Use coefficients to calculate LC50

coef(fdisagg.lm4)

(Intercept)       ldose 
  -3.030785    4.382655 

3.03/4.38

[1] 0.6917808

10^0.692

[1] 4.920395

For comparison to original analyses, we’ll use aggregated values, generate probit and logit values for each dose and fit linear regression.

finney <- read_csv("../data/finneyagg.csv")

Rows: 6 Columns: 8── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
dbl (8): dose, n, responding, propkill, l10dose, probit_finney, probit_raw, logit
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(finney)

Fit logit model

finn.lm<-lm(logit~l10dose, data=finney, na.action=na.omit)
tidy(finn.lm)

#Use coefficients to calculate LC50

coef(finn.lm)

(Intercept)     l10dose 
  -4.805097    7.003509 

4.805/7.004

[1] 0.6860366

# Back-transform to untransformed dose

10^.686

[1] 4.852885

Probit transformed values

finn.lm2<-lm(probit_raw~l10dose,data=finney, na.action=na.omit)
tidy(finn.lm2)

# use coefficients to calculate LC50
coef(finn.lm2)

(Intercept)     l10dose 
  -2.850557    4.150930 

2.85/4.15

[1] 0.686747

# Back-transform
10^0.687

[1] 4.864072

Note: Finney used probits with 5 added

finn.lm3<-lm(probit_finney~l10dose,data = finney, na.action=na.omit)
tidy(finn.lm3)

# solve for y=5
(5-2.14)/4.16

[1] 0.6875

10^.688

[1] 4.875285