如何从 nls SSlogis 逻辑增长模型计算倍增时间

Question

pop.ss <- nls(MRDRSLT ~ SSlogis(TIME, phi1, phi2, phi3), data = testing, na.action = na.omit)


theta <- coef(pop.ss)  #extracting coefficients
plot(MRDRSLT ~ TIME, data = testing, main = "Logistic Growth Model", 
     xlab = "Time", ylab = "MRD")  # Census data
curve(theta[1]/(1 + exp(-(x - theta[2])/theta[3])), add = T, col = "green")  # Fitted model


summary(pop.ss)

Formula: MRDRSLT ~ SSlogis(TIME, phi1, phi2, phi3)

Parameters:
      Estimate Std. Error t value Pr(>|t|)   
phi1 9.618e-05  6.935e-06  13.869  0.00516 **
phi2 2.480e+02  1.512e+01  16.403  0.00370 **
phi3 2.896e+01  2.392e+01   1.211  0.34960   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.197e-05 on 2 degrees of freedom

Number of iterations to convergence: 1 
Achieved convergence tolerance: 1.687e-06

我希望能够计算我感兴趣的变量的倍增时间。我如何从系数中提取这个。 https://rdrr.io/r/stats/SSlogis.html

Asym/(1+exp((xmid-input)/scal))

在指数增长模型中很容易计算，因为它是 1/rate。

t_doubling <- (1 / mu) * log10(2)

Answer 1

这对CrossValidated 可能更好。加倍时间不是逻辑曲线的内在特征，因为它取决于起始值（例如，任何从 K/2 以上开始的逻辑曲线，其中 K 是渐近值，都永远不会加倍，因为它永远不能超过 K ...)。如果起始值非常小，则倍增时间的公式与指数曲线相同（您的公式是错误的：倍增时间是 2/rate 的 自然 对数（log(2)/mu ，而不是log10(2)/mu)：

x0*exp(mu*T_dbl) = 2*x0
exp(mu*T_dbl) = 2
mu*T_dbl = log(2)
T_dbl = log(2)/mu

在上面给出的参数化中（使用比例参数而不是速率参数），这将是log(2)*scal（但再次注意它仅适用于（大约）从一个低密度并且没有超过大约Asym*0.3的水平[在此之后逻辑曲线看起来越来越线性而不是指数]）

Answer 2

一般来说，倍增时间仅对于指数曲线是恒定的，而对于逻辑曲线和其他曲线，我们只能计算瞬时倍增时间，该瞬时倍增时间随点而变，而不是单个值。

情况类似于斜率，仅对一条线是恒定的。对于其他（可微分）曲线，我们只能计算每个点的切线斜率，并且该斜率因点而异，而不是单个值。

一般情况下t点的瞬时倍增时间等于

1/(log₂ f(t))'

即f(t) 的以 2 为底的对数导数相对于在 t 处求值的 t 的倒数。这里 f(t) 是逻辑曲线或其他曲线。请注意，因为 log₂(x) = log(x) / log(2) 并且由于计算函数对数导数的规则，瞬时倍增时间等于以下，其中 log 是log base e 和 f 是 f(t) 的缩写，f' 是 f(t) 对 t 在 t 处求值的导数的缩写。

log(2) * f / f'

逻辑曲线

在逻辑的情况下，它满足一个众所周知的微分方程，可以帮助我们计算倍增时间。

f = SSlogis(Time, Asym, xmid, scal) # Asym / (1 + exp((xmid - x)/scal))
f' = (1/scal) * f * (1 - f / Asym)
doubling time = log(2) * f / f' = log(2) * scal / (1 - f/Asym)

在 R 中考虑 example(SSlogis) 中使用的 Chickwt.1 数据示例。这样我们就可以进行 R 计算了：

# example data
Chick.1 <- ChickWeight[ChickWeight$Chick == 1, ] # ChickWeight comes with R
Asym <- 368; xmid <- 14; scal <- 6
Time <- Chick.1$Time

f <- SSlogis(Time, Asym, xmid, scal)

fderiv <- (1/scal) * f * (1-f/Asym)
log(2) * f / fderiv # doubling time

# or equivalently
log(2) * scal / (1 - f/Asym) # doubling time

plot(log(2) * scal / (1 - f/Asym) ~ Time, type = "o", pch = 20, 
  ylab = "Instantaneous Doubling Time")

这是向量Time中每个时间点f的瞬时倍增时间。

特别考虑三种情况：

• 当 f 相对于 Asym 较小时，log(2) * scal / (1 - f/Asym) 的分母约为 1，因此倍增时间约为 log(2) * scal。此表达式与指数表达式相同（请参阅下一节）。
• 在逻辑 f = Asym/2 的拐点处，因此倍增时间为 2 * log(2) * scal。
• 当 f 接近 Asym 时，log(2) * scal / (1 - f/Asym) 的分母接近 0，因此倍增时间接近无穷大。

（剧情后续）

指数曲线

再举个例子，指数曲线的定义和微分方程为：

f = Asym * exp((xmid - Time)/scal)
f' = (1/scal) * f

所以即使不使用 R，考虑到分子和分母中的抵消，我们可以计算倍增时间：

log(2) * f / f'
= log(2) * scal

拟合

我们可以使用 nls 和 SSlogis 来估计 Chick.1 数据的逻辑参数：

fm.nls <- nls(weight ~ SSlogis(Time, Asym, xmid, scal), Chick.1)
fm.nls  # displays coefficients
log(2) * coef(fm.nls)[["scal"]]  # doubling time in initial portion of curve

在这个数据的情况下，它实际上是指数级的。如果它是指数的，那么绘制 log2（权重）与时间的关系图会给出一条近似直线，这就是这里的情况。在这种情况下，倍增时间是最佳拟合线斜率的倒数。

fm.lm <- lm(log2(weight) ~ Time, Chick.1)
1/coef(fm.lm)[2]  # doubling time

# seems like a straight line
plot(log2(weight) ~ Time, Chick.1)   
abline(fm.lm)