Estimating the ratio of two parameters and confidence bounds

Part 1 of this post is available here.

The dose at which some proportion of mortality occurs (e.g. 50%) has classically been a common summary statistic as it can be convenient to talk about mortality in terms of doses rather than model coefficients. This can be computed algebraically from the regression equation for each population as , which can be solved as . The lethal dose required for 50% mortality occurs when or .

-coef(m1)[1]/coef(m1)[2]

This is a compound parameter of two model parameters that each have their own variance and co-variance and it possible to estimate the distribution of this quantity analytically, but not straightforward (at least not for me).

The delta method

The delta method can be used to estimate variance of by assuming the quantity is multivariate normally distributed. The dose.p function in the MASS package can do this for us using the following code.

library(MASS)
# we can simply run dose.p(m1, cf = c(1,3) or more verbosely
cf = c(1, 3) # The terms in the coefficient vector giving the intercept and coefficient of (log-)dose
p = 0.5 # Probabilities at which to predict the dose needed.
eta <- family(m1)$linkfun(p)
b <- coef(m1)[cf]
x.p <- (eta - b[1L])/b[2L]
names(x.p) <- paste("p = ", format(p), ":", sep = "")
pd <-  -cbind(1, x.p)/b[2L]
SE <- sqrt(((pd %*% vcov(m1)[cf, cf]) * pd) %*% c(1, 1))

> # LD50
> exp(c(x.p - 1.96 * SE, x.p + 1.96 * SE))
[1]  7.294287 13.992988

Using the posterior distribution

I find the solution to this problem using the Bayesian framework more intuitive. The parameter covariation structure is captured in the posterior distribution, so we can sample parameter values from the posterior distribution and then do some algebra (solving for ) as before and estimate the LC50 for the control population

exp(c(x.p - 1.96 * SE, x.p + 1.96 * SE))
post <- extract(fit_mat)
logx50 <- -post$beta[, 1] / post$beta[, 3]
exp(mean(logx50))

> exp(quantile(logx50, c(0.025, 0.975)))
     2.5%     97.5% 
 7.330893 13.955303 

In general, solving for the distribution of arbitrary compound parameters is much more straightforward in a Bayesian context.

The ratio of two LC50 estimates

But perhaps the most compelling case for a Bayesian approach with dose responses is when it comes to estimating whether or not the LC50s of both populations are different (and then how different the responses are expressed as a resistance factor). You see many papers getting around this issue by somewhat informally inspecting the overlap of confidence intervals, which is likely to be an overly conservative test. The resistance factor is usually provided as the ratio of the LC50 values, which usually lacks a confidence interval.

# LC50 using maximum liklihood
# control
con = dose.p(m1, c(1, 3), p = 0.5)
exp(con)
t95 <- qnorm(0.025, lower.tail = FALSE)
exp(c(con - t95 * attr(con, "SE"), con + t95 * attr(con, "SE")))

# resistant
res = dose.p(m1, c(2, 4), p = 0.5)
res
exp(c(res - t95 * attr(con, 'SE'), res + t95 * attr(con, 'SE')))

# resistance factor
> exp(res[1])/exp(con[1])
4.369501 

When using this test statistic, I circumvent the issue of estimating a resistance ratio by performing a log-likelihood ratio test on the full model with a simplified models that assume universal slope, or intercept parameters, but this approach is a little less direct than ideal.

m1_full = glm(cbind(dead, alive) ~ pop*log(dose), family = binomial('logit'), 
         data = d )
m1_no_pop = glm(cbind(dead, alive) ~ log(dose), family = binomial('logit'), 
         data = d )
anova(m1_full, m1_no_pop, test = 'LRT')

> anova(m1_full, m1_no_pop, test = 'LRT')
Analysis of Deviance Table

Model 1: cbind(dead, alive) ~ pop * log(dose)
Model 2: cbind(dead, alive) ~ log(dose)
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1        10     28.145                          
2        12     59.612 -2  -31.468 1.469e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

But the issue of confidence bounds around the resistance factor remains.

Sampling from our posterior distribution makes this all look like a walk in the park.

# control LC50
exp(mean(-post$beta[, 1] / post$beta[, 3]))

# resistant LC50
exp(mean(-post$beta[, 2] / post$beta[, 4]))

# resistance factor with 95% credible intervals
x50_diff <- exp(-post$beta[, 2] / post$beta[, 4] - 
    -post$beta[, 1] / post$beta[, 3])

> mean(x50_diff)
[1] 4.551865

> quantile(x50_diff, c(0.025, 0.975))
    2.5%    97.5% 
2.737778 7.041024