R语言：基于昆虫解剖总数加权的三参数Logistic模型构建咨询

我是R语言新手，收集了多篇文献中的昆虫寄生虫感染数据，数据符合三参数Logistic函数（x为寄生虫浓度计数，y为昆虫感染比例，由感染数/解剖总数计算得出）。希望以昆虫解剖总数（解剖数量越多，结果越可靠）为权重拟合模型，获取三参数；尝试过nlme包，但不想按寄生虫计数分组（会丢失大量数据细节），询问是否可通过optim、lme或nlme等包实现。当前尝试的nls代码如下：

library(nlme)
# 3 params model are choosen based on visual interpretations by using SSlogis()
model_sslogis <- nls(Insect_larvae_infected_proportion_fromtotaldissected ~ SSlogis(Parasite_count_per_1uL, A, B, C),
                     data = dat,
                     algorithm = "port",
                     # weights = 1/sd_Insect_larvae_infected
)

summary(model_sslogis)

数据结构示例：

> dat <- read_excel('Data.xlsx') %>% 
+   # view() %>% 
+   glimpse()
Rows: 158
Columns: 15
$ Reference                                              <chr> "(Bryan & Southgate, 1988…
$ Parasite_count_per_1uL                                 <dbl> 0.9313223, 1.0464999, 1.1…
$ Insect_totaldissected                                  <dbl> 30, 50, 20, 36, 32, 40, 3…
$ Insect_infected_count                                  <dbl> 1, 4, 3, 7, 3, 6, 2…
$ Insect_larvae_infected_proportion_fromtotaldissected   <dbl> etc...

解决方案 #

方法1：在nls中直接加入权重（最简方案） #

感染比例服从二项分布，方差为p(1-p)/n（n为解剖总数），因此权重应与解剖总数成正比，拟合时更重视样本量大的观测。修改你的nls代码，将weights参数设为解剖总数：

library(nlme)
# 带权重的三参数Logistic模型
model_weighted <- nls(
  Insect_larvae_infected_proportion_fromtotaldissected ~ SSlogis(Parasite_count_per_1uL, A, B, C),
  data = dat,
  algorithm = "port",
  weights = Insect_totaldissected  # 用解剖总数作为权重
)

summary(model_weighted)

方法2：用optim手动拟合加权模型 #

如果nls出现收敛问题，可以手动定义目标函数（加权残差平方和最小化），用optim进行拟合：

1. 定义三参数Logistic函数：

logistic_3p <- function(x, A, B, C) {
  # A=渐近上限，B=中点x值，C=斜率参数
  A / (1 + exp(-(x - B)/C))
}

2. 定义加权残差平方和目标函数：

weighted_resid_ss <- function(params, x, y, weights) {
  A <- params[1]
  B <- params[2]
  C <- params[3]
  pred <- logistic_3p(x, A, B, C)
  sum(weights * (y - pred)^2)  # 加权残差平方和
}

3. 用SSlogis获取初始参数，再调用optim拟合：

# 先获取初始参数估计
init_params <- coef(nls(
  Insect_larvae_infected_proportion_fromtotaldissected ~ SSlogis(Parasite_count_per_1uL, A, B, C),
  data = dat,
  algorithm = "port"
))

# 带边界约束的optim拟合
optim_fit <- optim(
  par = init_params,
  fn = weighted_resid_ss,
  x = dat$Parasite_count_per_1uL,
  y = dat$Insect_larvae_infected_proportion_fromtotaldissected,
  weights = dat$Insect_totaldissected,
  method = "L-BFGS-B",
  lower = c(0, min(dat$Parasite_count_per_1uL), 0),  # 参数边界：A>0，B在x范围内，C>0
  upper = c(1, max(dat$Parasite_count_per_1uL), Inf)
)

# 查看拟合参数
optim_fit$par

方法3：用nlme包拟合无分组加权模型 #

若想用nlme但不按寄生虫计数分组，可以指定无随机效应的固定模型，同时加入权重：

library(nlme)
# nlme无分组加权拟合
nlme_fit <- nlme(
  Insect_larvae_infected_proportion_fromtotaldissected ~ SSlogis(Parasite_count_per_1uL, A, B, C),
  data = dat,
  fixed = A + B + C ~ 1,  # 仅固定效应，无分组
  random = pdDiag(A + B + C ~ 1),  # 随机效应设为对角矩阵（等价于无随机效应）
  weights = varFixed(~ 1/Insect_totaldissected)  # 权重与解剖总数成反比，效果同加权nls
)

summary(nlme_fit)

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