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

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



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


theta <- coef(pop.ss)  #extracting coefficients
plot(MRDRSLT ~ TIME,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,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)

解决方法

这可能对 CrossValidated 更好。倍增时间不是逻辑曲线的内在特征，因为它取决于起始值（例如，任何开始于 K/2 以上的逻辑曲线，其中 K 是渐近值，将永远翻倍，因为它永远不能超过 K ......）。如果起始值非常小，则倍增时间的公式与指数曲线的公式相同（您的公式是错误的：倍增时间是 2/rate ({{1 }}，而不是 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 的水平 [超过这个水平，逻辑曲线看起来越来越线性而不是指数])

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

这种情况类似于斜率，它只对一条直线保持不变。对于其他（可微分）曲线，我们只能计算每个点处切线的斜率，并且该斜率因点而异，而不是单个值。

一般来说，t 点的瞬时倍增时间等于

1/(log₂ f(t))'

即f(t) 的以 2 为底的对数导数相对于在 t 处求得的 t 的倒数。这里 f(t) 是逻辑或其他曲线。请注意，因为 log₂(x) = log(x) / log(2) 并且由于评估函数对数的导数的规则，瞬时倍增时间等于以下，其中 log 是对数底数 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,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,scal),Chick.1)
fm.nls  # displays coefficients
log(2) * coef(fm.nls)[["scal"]]  # doubling time in initial portion of curve

在这个数据的情况下，它实际上是指数级的。如果它是指数的，那么绘制 log2(weight) vs. Time 会给出大约一条直线，这里就是这种情况。在这种情况下，倍增时间是最佳拟合线斜率的倒数。

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)