如何在R中用Trimmed Spearman-Karber法计算蜜蜂24/48小时LD50及置信区间

实验背景与数据 #

针对不同农药暴露浓度（0、3、30、300ng）的蜜蜂，每4小时测定一次存活率，用于计算LD50（致死浓度，即杀死50%受试生物的剂量）。实验数据如下：

Control <- c(100, 100, 100, 96, 96, 96, 96, 72, 60, 60, 60, 60, 60, 52, 48, 48, 40, 40)
`300ng` <- c(100.00, 100.00, 100.00, 96.30, 96.30, 92.59, 92.59,70.37, 62.96, 44.44, 40.74, 37.04, 29.63, 25.93, 25.93,22.22, 11.11, 11.11)
`30ng` <- c(100.00, 96.30, 96.30, 96.30, 96.30, 96.30, 96.30, 85.19, 81.48, 77.78, 74.07, 74.07, 74.07, 70.37, 70.37, 70.37, 70.37, 62.96)
`3ng` <- c(100.00, 100.00, 100.00, 100.00, 100.00, 96.30, 85.19, 74.07, 70.37, 66.67, 66.67, 66.67, 66.67, 59.26, 59.26, 59.26, 59.26, 55.56)
HoursExp <- c(0, 4, 8, 12, 16, 20, 24, 28, 32, 36, 40, 44, 48, 52, 56, 60, 64, 68)

已绘制存活率随暴露时间变化曲线：

plot(`300ng`~HoursExp, type="l", col=1)
points(`30ng`~HoursExp, type="l", col=2)
points(`3ng`~HoursExp, type="l", col=3)
points(Control~HoursExp, type="l", col=4)

需求 #

计算蜜蜂暴露24小时和48小时的LD50及置信区间，优先使用Trimmed Spearman-Karber法（tsk函数），也可采用其他适用方法，已安装并加载tsk和drc包。

解决方案 #

步骤1：提取目标时间点数据 #

提取24小时（HoursExp=24，对应第7个元素）和48小时（HoursExp=48，对应第13个元素）的存活率数据：

# 提取24小时存活率
surv_24 <- data.frame(
  dose = c(0, 3, 30, 300),
  surv = c(Control[7], `3ng`[7], `30ng`[7], `300ng`[7])
)

# 提取48小时存活率
surv_48 <- data.frame(
  dose = c(0, 3, 30, 300),
  surv = c(Control[13], `3ng`[13], `30ng`[13], `300ng`[13])
)

步骤2：校正死亡率（Abbott公式） #

用Abbott公式校正对照组自然死亡的影响，校正死亡率 = max(0, (对照组存活率 - 处理组存活率)/对照组存活率)：

# 24小时校正死亡率
control_surv_24 <- surv_24$surv[surv_24$dose == 0]
surv_24$corr_mort <- pmax(0, (control_surv_24 - surv_24$surv)/control_surv_24)

# 48小时校正死亡率
control_surv_48 <- surv_48$surv[surv_48$dose == 0]
surv_48$corr_mort <- pmax(0, (control_surv_48 - surv_48$surv)/control_surv_48)

方法1：Trimmed Spearman-Karber法（tsk包） #

tsk函数需要剂量对数、校正死亡率（比例）和每组样本量。假设每组样本量为27（从存活率百分比反推），代码如下：

library(tsk)

# 计算24小时LD50及置信区间
tsk_result_24 <- tsk(
  x = log10(surv_24$dose[-1]),  # 去掉对照组，取剂量对数
  p = surv_24$corr_mort[-1],    # 去掉对照组的校正死亡率
  n = rep(27, 3)                # 每组样本量
)

# 计算48小时LD50及置信区间
tsk_result_48 <- tsk(
  x = log10(surv_48$dose[-1]),
  p = surv_48$corr_mort[-1],
  n = rep(27, 3)
)

# 输出结果
print("24小时LD50结果：")
print(tsk_result_24)
print("48小时LD50结果：")
print(tsk_result_48)

方法2：剂量-反应模型拟合（drc包） #

使用drc包的drm函数拟合log-logistic模型，再计算LD50：

library(drc)

# 整理24小时死亡数据（死亡数=总样本量*(1-存活率/100)）
mort_24 <- data.frame(
  dose = surv_24$dose[-1],
  dead = round((1 - surv_24$surv[-1]/100)*27),
  total = rep(27, 3)
)

# 拟合模型并计算24小时LD50
model_24 <- drm(dead/total ~ dose, weights = total, data = mort_24, fct = LL.2(), type = "binomial")
ld50_24 <- ED(model_24, 50, interval = "confidence")

# 整理48小时死亡数据
mort_48 <- data.frame(
  dose = surv_48$dose[-1],
  dead = round((1 - surv_48$surv[-1]/100)*27),
  total = rep(27, 3)
)

# 拟合模型并计算48小时LD50
model_48 <- drm(dead/total ~ dose, weights = total, data = mort_48, fct = LL.2(), type = "binomial")
ld50_48 <- ED(model_48, 50, interval = "confidence")

# 输出结果
print("24小时LD50结果：")
print(ld50_24)
print("48小时LD50结果：")
print(ld50_48)

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