Search code examples
rlogistic-regressionglm

Lethal Dose (LD50) for the different factors in a logistic regression

I have a logistic regression model that contains a continuous variable and a qualitative variable with various levels/factors. Is there a way to calculate LD50 for each of the factors?

library(dplyr)
library(ggeffects)
library(MASS)

# Produce a table with data
set.seed(321)

x1 <- rnorm(400)
x2 <- 3*x1^3 + 2*x1^2 + x1 + 1
x3 <- 1/(1 + exp(-x2))
x4 = rbinom(400, 1, x3) 

dt <- tibble(binary = x4,
             continuous = x1,
             qualitative = sample(c("white", "pink", "lime"), 400, TRUE))

# Logistic model
mdl <- glm(binary ~ continuous + qualitative, family = "binomial", data = dt)

# Calculate LD50, which gives the lethal dose for the whole model
dose.p(mdl) # -0.58 ± 0.11 se

# Plot regression model by factor of the qualitative variable
newdat <- ggpredict(mdl, terms = c("continuous[all]", "qualitative"))
plot(newdat)

How can I calculate LD50 ± SE for each of the factors (LD50 for lime, LD50 for pink, etc.)?

enter image description here

Solution

  • This below fits a model for binary vs continuous for every factor and pulls out the LD50 and SE:

    get_LD50 = function(fit){
   data.frame(
   LD50 = dose.p(fit)[1],
   SE = attributes(dose.p(fit))$SE[,1]
   )
}

dt %>% group_by(qualitative) %>% 
do(get_LD50(glm(binary ~ continuous, family = "binomial", data = .)))

# A tibble: 3 x 3
# Groups:   qualitative [3]
  qualitative   LD50    SE
  <chr>        <dbl> <dbl>
1 lime        -0.578 0.115
2 pink        -0.479 0.104
3 white       -0.392 0.116

    If you would need something actually more complicated on the statistical side, you have to clarify