Calculation of log-logit or log-probit models according to Finney using R

Question

I'm trying to implement the logit/probit model derivation as introduced by Finney using sample data from http://dge.stanford.edu/SCOPE/SCOPE_12/SCOPE_12.html, chapter 6 (this links to a pdf), page 130, table 6.3. Results are similar but not identical. Finney calculated LD50 = 4.85.

Here is the source code and output of my R program:

# Finney, 1952
dosis <- c(2.6,3.8,5.1,7.7,10.2)
nges <- c(50,48,46,49,50)
nok <- c(6,16,24,42,44)

edx <- function(dosis=NA, nges=NA, nok=NA, lk='logit') {
  require(MASS)
  # weights w according to Finney 1952:
  p <- nok/nges
  q <- 1-p
  w <- nges*p*q

  # logit/probit model
  edx.data <- data.frame(dosis, nges, nok)
  glm.logit <- glm(cbind(nok,nges-nok) ~ dosis, family=binomial(lk), data=edx.data, weights=w)

  # Calculation of EDx data
  r <- dose.p(glm.logit,p=seq(0.1,0.9,0.2))

  # Statistical summary
  d <- data.frame(x=c(NA,NA,NA,NA))
  rownames(d) <- c('Deg. of freedom','Deviance','1-chi.square','Significant difference between fits and observations')
  d$x[1] <- df.residual(glm.logit)
      d$x[2] <- deviance(glm.logit)
  d$x[3] <- 1-pchisq(d$x[2],d$x[1])
      d$x[4] <- ifelse(d$x[3]>0.05,'No','Yes')

  # Printing of data tables
  print(edx.data)
  writeLines('')
  print(r)
  writeLines('')
  print(d)

  return(glm.logit)
}

g <- edx(dosis=dosis, nges=nges, nok=nok, lk='logit')

The numerical output is as follows:

  dosis nges nok
1   2.6   50   6
2   3.8   48  16
3   5.1   46  24
4   7.7   49  42
5  10.2   50  44

             Dose         SE
p = 0.1: 1.171955 0.22512512
p = 0.3: 3.615213 0.12064635
p = 0.5: 5.148754 0.09792762
p = 0.7: 6.682294 0.13521132
p = 0.9: 9.125553 0.24565236

                                                                        x
Deg. of freedom                                                         3
Deviance                                                 33.3237083095658
1-chi.square                                         2.75205023214653e-07
Significant difference between fits and observations                  Yes

Thus, LD50 would be 5.149 according to this R program (vs. 4.85 as calculated by Finney). What is wrong here?

Thank you very much for any help!

Answer 1

You have made two (three?) mistakes in trying to reproduce those results.

First, as @whuber pointed out, you are using logistic regression directly on the binomial outcome instead of linear regression for the observed log-odds. While the former is the appropriate approach by today's standards, that's not what Finney's method, which was proposed before the advent of logistic regression, is doing.

Second, you are fitting a logistic model with $dosis$ as the predictor instead of $\log(dosis)$ as Finney does.

Third, you are using weights for the logistic model which are not needed, because the whole point of the logistic model is that it incorporates the variance structure appropriately.

Doing the analysis with the correct logistic regression gives almost the same results as Finney's method

library(MASS)
dosis <- c(2.6,3.8,5.1,7.7,10.2)
nges <- c(50,48,46,49,50)
nok <- c(6,16,24,42,44)
mm <- glm(cbind(nok, nges - nok) ~ log(dosis), family=binomial)
logED50 <- dose.p(mm, p=0.5)
exp(c(logED50))

Output:

p = 0.5: 
4.828918 

And here is Finney's method:

w <- nges*(nok/nges)*(1-nok/nges)
m <- lm(log(nok/(nges-nok)) ~ log(dosis), weights=w)

logED50 <- -coef(m)[1] / coef(m)[2]
exp(c(logED50))

Output:

(Intercept) 
   4.845005 

Answer 2

The formulas on that page are clearly for a least-squares estimate, not a maximum likelihood estimate (as explained three pages earlier at (6.13)). Because glm uses ML, it appears you are not actually reproducing that model. If you really want to reproduce it, you either need to understand how the formulas connect with R's lm functionality and apply lm or else you need to implement the formulas as they are given. Until you do one of those things, it's not valid to compare your output (digit by digit) with the results in that chapter.

---

It is perhaps of interest that glm's $LD_{50}$ of $5.15$ is well within the $95\%$ confidence interval of $(4.37, 5.37)$ for the Finney estimate, indicating the two methods do not significantly disagree. It is more striking that glm reports such tiny SEs (relatively) for its estimates. That high precision does not seem justified by the dataset. If your concern is with the data and the estimates, rather than the methodology, you would do well to pay attention to the deviance test and work to understand the nature of the lack of fit it indicates.