如何在ggplot的logit模型曲线中高亮50%概率点及置信区间？

Let's adjust your existing code to include the 50% probability marker, vertical reference line, and confidence interval for the corresponding temperature. Here's the modified code with clear explanations:

Step 1: Calculate the 50% Probability Temperature & Its Confidence Interval #

For a logit model, a 50% success probability corresponds to a linear predictor value of 0 (since logit(p) = ln(p/(1-p)) = 0 when p=0.5). We solve for the temperature that hits this value, then compute its confidence interval using the delta method:

data <- data.frame(MaxTemp = c(53.2402, 59.01004,51.42602,41.53883,44.70763,53.90285,51.130318,54.5929,43.697559,49.772446,54.902222,52.720528,58.782608,47.680374,48.30313,56.10921,57.660324,46.387924,60.503147,53.803177,52.27771,58.58555,55.74136,49.04505,46.816269,52.58295,52.751373,56.209747,51.733894,51.424305,50.74564,47.046513,53.030407,56.68752,56.639351,53.526585,51.562313), 
                   Success=c(1,1,1,0,0,1,1,1,0,0,1,1,1,0,0,1,1,0,1,1,1,1,1,1,0,1,1,1,1,1,1,0,1,1,1,1,1))

TempProbitModel <- glm(Success ~ MaxTemp, data=data, family=binomial(link="logit"))
temp.data <- data.frame(MaxTemp = seq(40, 62, 0.5))
predicted.data <- as.data.frame(predict(TempProbitModel, newdata = temp.data, type="link", se=TRUE))
new.data <- cbind(temp.data, predicted.data)
std <- qnorm(0.95 / 2 + 0.5)
new.data$ymin <- TempProbitModel$family$linkinv(new.data$fit - std * new.data$se)
new.data$ymax <- TempProbitModel$family$linkinv(new.data$fit + std * new.data$se)
new.data$fit <- TempProbitModel$family$linkinv(new.data$fit)

# Calculate temperature where success probability = 50%
temp_50 <- -coef(TempProbitModel)[1]/coef(TempProbitModel)[2]

# Compute confidence interval for temp_50 using delta method
vc <- vcov(TempProbitModel)
se_temp_50 <- sqrt( vc[1,1]/(coef(TempProbitModel)[2]^2) + 
                    (2*coef(TempProbitModel)[1]*vc[1,2])/(coef(TempProbitModel)[2]^3) + 
                    (coef(TempProbitModel)[1]^2*vc[2,2])/(coef(TempProbitModel)[2]^4) )
ci_temp_50 <- temp_50 + c(-1,1)*qnorm(0.975)*se_temp_50

Step 2: Update the Plot with Highlighted Elements #

Now we'll add the 50% point, vertical line, and confidence interval annotations to your ggplot:

library(ggplot2)

TempProb <- ggplot(data, aes(x=MaxTemp, y=Success)) +
  geom_point() +
  geom_ribbon(data=new.data, aes(y=fit, ymin=ymin, ymax=ymax), alpha=0.5) +
  geom_line(data=new.data, aes(y=fit)) +
  # Highlight the 50% probability point
  geom_point(aes(x=temp_50, y=0.5), color="red", size=4, shape=19) +
  # Vertical line to x-axis for the 50% temperature
  geom_vline(xintercept=temp_50, color="red", linetype="dashed", linewidth=1) +
  # Add confidence interval for the temperature (dotted lines + annotation)
  geom_vline(xintercept=ci_temp_50, color="darkred", linetype="dotted") +
  annotate("text", x=mean(ci_temp_50), y=0.05, 
           label=paste0("95% CI: ", round(ci_temp_50[1],2), " to ", round(ci_temp_50[2],2)),
           color="darkred", size=3.5) +
  # Label the 50% point
  annotate("text", x=temp_50 + 1, y=0.5, 
           label=paste0("50% Probability\nTemp: ", round(temp_50,2)),
           color="red", hjust=0) +
  labs(x="Peak Temperature", y="Probability of Success") +
  theme_minimal()

TempProb

Key Notes: #

• 50% Temperature Calculation: We solve intercept + slope*Temp = 0 (since logit(0.5)=0) to get temp_50 = -intercept/slope.
• Confidence Interval: The delta method approximates the standard error of our derived temperature value using the model's variance-covariance matrix, then computes the 95% CI with qnorm(0.975).
• Plot Annotations: Red elements make the 50% point stand out, dashed/dotted lines clarify the temperature and its CI, and text labels ensure viewers understand each component.

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