两曲线间区域面积估算技术求助

要计算蓝色（t<0的代际间隔PDF区域）和绿色（两条曲线间的全部区域）的面积，需分三步实现：定位曲线交点、分区间计算积分、更新可视化代码。

1. 定位两条曲线的交点（x>0区域） #

潜伏期的对数正态分布在x≤0时概率密度为0，仅需在x>0区间寻找两条PDF的交点。用uniroot函数求解serial_pdf(x) - incubation_pdf(x) = 0的根：

# 定义两条PDF的差值函数
pdf_diff <- function(x) serial_pdf(x) - incubation_pdf(x)

# 寻找x>0区间的交点（根据分布参数，交点大致在1-15区间内）
intersection <- uniroot(pdf_diff, interval = c(1, 15))$root
cat("两条曲线的交点在x =", round(intersection, 2), "天\n")

2. 计算各区域面积 #

蓝色区域（t<0的代际间隔PDF面积） #

你已完成这部分计算，结果即变量pre_symptomatic_auc。

绿色区域（两条曲线间的全部面积） #

绿色区域需分两段计算：

• 0到交点：代际间隔PDF高于潜伏期PDF的区域，积分serial_pdf(x) - incubation_pdf(x)
• 交点到+∞：潜伏期PDF高于代际间隔PDF的区域，积分incubation_pdf(x) - serial_pdf(x)

代码实现：

# 计算0到交点的绿色区域面积
green_area1 <- integrate(function(x) serial_pdf(x) - incubation_pdf(x), lower = 0, upper = intersection)$value

# 计算交点到无穷大的绿色区域面积
green_area2 <- integrate(function(x) incubation_pdf(x) - serial_pdf(x), lower = intersection, upper = Inf)$value

# 绿色区域总面积
total_green_area <- green_area1 + green_area2

# 输出结果
cat("蓝色区域面积（t<0代际间隔）：", round(pre_symptomatic_auc, 4), "\n")
cat("绿色区域总面积（两条曲线间）：", round(total_green_area, 4), "\n")
cat("蓝色+绿色总面积：", round(pre_symptomatic_auc + total_green_area, 4), "\n")

3. 更新可视化代码，高亮绿色区域 #

用geom_ribbon填充两条曲线之间的区域，替代原有仅绘制蓝色区域的代码：

# 更新绘图逻辑
ggplot(plot_data, aes(x = Time, y = PDF, color = Type)) +
  geom_line(size = 1) +
  # 蓝色区域：t<0的代际间隔PDF下方
  geom_area(data = subset(plot_data, Type == "Serial Interval (Normal)" & Time < 0),
            aes(y = PDF), fill = "blue", alpha = 0.3) +
  # 绿色区域：两条曲线之间的全部区域
  geom_ribbon(data = data.frame(Time = time_range, 
                                min_pdf = pmin(incubation_values, serial_values),
                                max_pdf = pmax(incubation_values, serial_values)),
              aes(x = Time, ymin = min_pdf, ymax = max_pdf),
              fill = "green", alpha = 0.3) +
  labs(
    title = "Comparison of Serial Interval and Incubation Period",
    x = "Time (days)",
    y = "Probability Density",
    color = "Distribution"
  ) +
  theme_minimal()

完整修改后代码 #

library(ggplot2)
library(dplyr)
library(stats)   

# 定义分布参数
# 潜伏期（对数正态）：均值4.4，中位数4.1，25分位3.2，75分位5.3
incubation_meanlog <- log(4.4)
incubation_sdlog <- (log(5.3) - log(3.2)) / (2 * qnorm(0.75))

# 代际间隔（正态）：均值4.6，标准差4.4
serial_mean <- 4.6
serial_sd <- 4.4

# 定义概率密度函数
incubation_pdf <- function(x) dlnorm(x, meanlog = incubation_meanlog, sdlog = incubation_sdlog)
serial_pdf <- function(x) dnorm(x, mean = serial_mean, sd = serial_sd)

# 计算t<0的代际间隔面积（蓝色区域）
pre_symptomatic_auc <- integrate(serial_pdf, lower = -Inf, upper = 0)$value
total_serial_auc <- integrate(serial_pdf, lower = -Inf, upper = Inf)$value
proportion_pre_symptomatic <- pre_symptomatic_auc / total_serial_auc
cat("Proportion of pre-symptomatic transmission:", proportion_pre_symptomatic, "\n")

# 新增：计算曲线交点与绿色区域面积
pdf_diff <- function(x) serial_pdf(x) - incubation_pdf(x)
intersection <- uniroot(pdf_diff, interval = c(1, 15))$root
cat("两条曲线的交点在x =", round(intersection, 2), "天\n")

green_area1 <- integrate(function(x) serial_pdf(x) - incubation_pdf(x), lower = 0, upper = intersection)$value
green_area2 <- integrate(function(x) incubation_pdf(x) - serial_pdf(x), lower = intersection, upper = Inf)$value
total_green_area <- green_area1 + green_area2

cat("蓝色区域面积（t<0代际间隔）：", round(pre_symptomatic_auc, 4), "\n")
cat("绿色区域总面积（两条曲线间）：", round(total_green_area, 4), "\n")
cat("蓝色+绿色总面积：", round(pre_symptomatic_auc + total_green_area, 4), "\n")

# 可视化准备
time_range <- seq(-10, 20, by = 0.1)
incubation_values <- incubation_pdf(time_range)
serial_values <- serial_pdf(time_range)

plot_data <- data.frame(
  Time = rep(time_range, 2),
  PDF = c(incubation_values, serial_values),
  Type = rep(c("Incubation Period (Lognormal)", "Serial Interval (Normal)"), each = length(time_range))
)

# 绘图
ggplot(plot_data, aes(x = Time, y = PDF, color = Type)) +
  geom_line(size = 1) +
  geom_area(data = subset(plot_data, Type == "Serial Interval (Normal)" & Time < 0),
            aes(y = PDF), fill = "blue", alpha = 0.3) +
  geom_ribbon(data = data.frame(Time = time_range, 
                                min_pdf = pmin(incubation_values, serial_values),
                                max_pdf = pmax(incubation_values, serial_values)),
              aes(x = Time, ymin = min_pdf, ymax = max_pdf),
              fill = "green", alpha = 0.3) +
  labs(
    title = "Comparison of Serial Interval and Incubation Period",
    x = "Time (days)",
    y = "Probability Density",
    color = "Distribution"
  ) +
  theme_minimal()

关键说明 #

• uniroot的interval参数可根据分布参数调整，确保覆盖实际交点范围
• geom_ribbon通过pmin和pmax自动获取每个x点的最小/最大PDF值，无需手动拆分区间
• 积分时分区间计算差值，避免出现负面积，保证结果为实际区域的绝对值

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