如何在含两水平因子的逻辑回归中求概率0.5对应的协变量值？

Ah, I see the issue—dose.p() is only giving you a single value because it defaults to using the reference level of your PadLen factor. Since PadLen has two distinct levels, you'll need to calculate the speed value that corresponds to a 50% probability separately for each level. Here's how to do it properly:

1. Understand the Model's Logit Equation #

Your binomial GLM uses the logit link, so the equation looks like this:

log(p/(1-p)) = β₀ + β₁*speed + β₂*PadLen250

(Note: This assumes PadLen50 is the reference level; adjust if your reference is PadLen250 instead.) When the probability p=0.5, the left-hand side simplifies to log(1) = 0, so we can rearrange to solve for speed:

speed = -(β₀ + β₂*PadLen_level) / β₁

2. Calculate Exact Values Using Model Coefficients #

First, pull out your model's coefficients:

coefs <- coef(model)
coefs

Suppose your output looks like this (example values):

(Intercept)       speed    PadLen250 
   -12.3          0.34         -2.1 

Then compute the 50% probability speed for each level:

• For PadLen=50 (reference level, so PadLen250=0):

  speed_p50_pad50 <- -coefs["(Intercept)"] / coefs["speed"]
  

• For PadLen=250 (where the factor coefficient is active):

  speed_p50_pad250 <- -(coefs["(Intercept)"] + coefs["PadLen250"]) / coefs["speed"]
  

3. Use dose.p() with Factor Level Adjustments #

You can also adapt dose.p() to target each factor level by specifying the relevant coefficient combination:

# For PadLen=50 (reference level)
dose.p(model, p=0.5, cf = c(coefs["(Intercept)"], coefs["speed"], 0))

# For PadLen=250
dose.p(model, p=0.5, cf = c(coefs["(Intercept)"], coefs["speed"], coefs["PadLen250"]))

4. Inverse Prediction with predict() (For More Complex Models) #

If you want a more flexible approach (useful if your model has interactions or non-linear terms), you can generate predicted probabilities across a range of speeds and find where they cross 0.5:

library(dplyr)

# Create a dense grid of speed values
speed_grid <- seq(min(myData$speed), max(myData$speed), length.out = 1000)

# Predict probabilities for both PadLen levels
pred_df <- expand.grid(speed = speed_grid, PadLen = c("50", "250")) %>%
  mutate(prob = predict(model, newdata = ., type = "response"))

# Find the speed closest to 50% probability for each level
speed_at_p50 <- pred_df %>%
  group_by(PadLen) %>%
  slice(which.min(abs(prob - 0.5))) %>%
  select(PadLen, speed)

print(speed_at_p50)

Why Your Original dose.p() Call Only Gave One Result #

By default, dose.p() uses the model's coefficients as they are, which corresponds to the reference level of your factor (since the factor's coefficient is 0 for the reference group). To get the value for the non-reference level, you need to include that level's coefficient in the calculation, as shown in the methods above.

内容的提问来源于stack exchange，提问作者Benjamin