如何对 y = a*(1- exp(-b*t) [关闭]

Question

我有时间序列数据，其中值随着时间（t）增加然后变得稳定：我想拟合指数曲线

NPQ_Lss = a*(1- exp(-b*t)

。稍后我想比较人群（

P1

和

P2

）之间的拟合曲线，看看它们是否有显着差异。

我怎样才能在 R 中做到这一点？我了解我的数据点很少的限制。

我的数据如下所示：

t	NPQ_Lss	人口
-2	1.777	P2
0	2.224	P2
3	2.492	P2
5	2.384	P2
7	2.407	P2
-2	1.719	P1
0	2.191	P1
3	2.418	P1
5	2.303	P1
7	2.403	P1

Answer 1

使用最后注释中的可重现数据，我们发现原始模型不太适合，因此我们将使用 3 参数自启动逻辑模型

SSlogis

，如图所示。

fm1

和

fm2

是人口

P1

和

P2

的模型。

fm3

是针对每个群体单独拟合的模型，

fm4

是对两个群体都具有单一拟合的模型。

我们可以使用

anova

来测试

fm3

和

fm4

是否有显着差异。
p 值为 0.732，它们没有显着差异。

DF$Population <- factor(DF$Population)  # important

fm1 <- nls(NPQ_Lss ~ SSlogis(t, Asym, mid, scal), DF, subset = Population == "P1")
fm2 <- update(fm1, subset = Population == "P2")

fo <- NPQ_Lss ~ Asym[Population] / 
  (1 + exp((mid[Population] - t)/scal[Population]))
st <- Map(c, coef(fm1), coef(fm2))
fm3 <- nls(fo, DF, start = st) # fit each pop separately
fm4 <- update(fm1, subset = TRUE) # single fit to both

anova(fm3, fm4)
## Analysis of Variance Table
## 
## Model 1: NPQ_Lss ~ Asym[Population]/(1 + exp((mid[Population] - t)/scal[Population]))
## Model 2: NPQ_Lss ~ SSlogis(t, Asym, mid, scal)
##   Res.Df Res.Sum Sq Df     Sum Sq F value Pr(>F)
## 1      4   0.018712                             
## 2      7   0.025006 -3 -0.0062939  0.4485  0.732

plot(NPQ_Lss ~ t, DF, col = Population, pch = 20, cex = 2)
lines(fitted(fm3) ~ t, DF, subset = Population == "P1", col = 1)
lines(fitted(fm3) ~ t, DF, subset = Population == "P2", col = 2)
lines(fitted(fm4) ~ t, DF, subset = Population == "P1", col = 3)

点就是数据。线条是拟合值。黑色是P1。红色是P2。绿色对两种人群均适用 (

fm4

)。

旧

这是原来的答案，但看了之后似乎不太合适。

假设最后的注释中显示的数据可重复显示，请使用单独的参数拟合每个数据，然后使用相同的参数拟合两个数据，然后进行方差分析。它给出的 P 值为 0.998，这并不显着，因此我们可以对两者使用相同的参数。

# separate parameters
DF$Population <- factor(DF$Population)
fo1 <- NPQ_Lss ~ a[Population]*(1 - exp(-b[Population]*t))
fm1 <- nls(fo1, DF, start = list(a = c(-2, -2), b = c(-0.2, -0.2)))

# same parameters
fo2 <- NPQ_Lss ~ cbind(a = (1 - exp(-b*t)))
fm2 <- nls(fo2, DF, start = list(b = -0.2), alg = "plinear")

anova(fm2, fm1)

## Analysis of Variance Table
##
## Model 1: NPQ_Lss ~ cbind(a = (1 - exp(-b * t)))
## Model 2: NPQ_Lss ~ a[Population] * (1 - exp(-b[Population] * t))
##   Res.Df Res.Sum Sq Df   Sum Sq F value Pr(>F)
## 1      8     25.136                           
## 2      6     25.134  2 0.001791   2e-04 0.9998

注意

DF <- data.frame(
  t = rep(c(-2L, 0L, 3L, 5L, 7L), 2),
  NPQ_Lss = c(1.777, 2.224, 2.492, 2.384, 2.407, 1.719, 2.191, 2.418, 2.303, 2.403),
  Population = rep(c("P2", "P1"), each = 5L)
)

Answer 2

原则上，

library(nlme)
fit <- gnls(NPQ_Lss ~ SSasympOrig(t, a, log_b),
            data = DF,
            params = a + log_b ~ Population)

应该使用特定于人群的参数来拟合模型。 （

SS

代表“自启动”，因此如果您的数据表现相当良好，则不需要提供起始值。）

但是，实际上并非如此：这里的主要问题是，您的数据实际上与曲线

a*(1-exp(-b*t))

非常匹配，该曲线假设曲线在时间 0 处通过 y=0。（您的曲线是大部分到达时间

t=0

的渐近线。）

使用

SSasymp

将数据拟合为一个整体，不假设曲线穿过原点，工作正常：

fit <- gnls(NPQ_Lss ~ SSasymp(t, a, y0, log_b),
            data = DF)

不幸的是，这个规范，它应该工作，却不起作用（“提供了3个起始值，需要6个”）

fit2 <- gnls(NPQ_Lss ~ SSasymp(t, a, y0, log_b),
            data = DF,
            params = a + y0 + log_b ~ Population)

让我们尝试给它初始值

fit2 <- gnls(NPQ_Lss ~ SSasymp(t, a, y0, log_b),
            data = DF,
            params = a + y0 + log_b ~ Population,
            start = c(rbind(coef(fit), 0)))

start

值是一种将起始值设置为

{a, 0, y0, 0, log_b, 0}

的偷偷摸摸的方法，知道该模型将根据

{"a for population 1", "difference in a between pops", "y0 for P1", "y0 diff", "log_b for P1", "log_b diff"}

进行参数化...

型号比较：

newdf <- expand.grid(Population = c("P1", "P2"), t = seq(-3, 10, length = 51))
df_pred <- (list(orig_0=fm1, orig_pop = fm2, flex_0 = fit, flex_pop = fit2)
    |> purrr::map_dfr(~ data.frame(newdf, NPQ_Lss = predict(., newdata = newdf)), .id = "model")
)

library(ggplot2); theme_set(theme_bw())
ggplot(df_pred,
       aes(t, NPQ_Lss, colour = Population)) +
    geom_line(aes(linetype = model)) +
    geom_point(data = DF)
ggsave("SO77667833.png")

坚持通过 (0,0) 拟合曲线的模型给出了技术上正确的预测（上凹/加速曲线）（这是该模型对这些数据的最小二乘拟合），但很荒谬；不限制拟合通过原点的模型 (

SSasymp

) 会给出合理的结果。

两种群体模型 (

summary()

) 的 

fit2

将为您提供各个参数差异的 p 值；

anova(fit, fit2)

按照@G.Grothendieck的建议将联合测试参数的差异。