0

我正在尝试找到适合我的数据的方法。但到目前为止还没有运气。尝试了对数的,与 drc 包不同的.. 但我确信一定有更好的我只是不知道类型。另一方面,我将不胜感激有关如何进行一般曲线狩猎的建议。

library(drc)    
df<-structure(list(x = c(10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 
        20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 
        36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 
        52, 53, 54, 55), y = c(0.1066, -0.6204, -0.2028, 0.2621, 0.4083, 
        0.4497, 0.6343, 0.7762, 0.8809, 1.0029, 0.8089, 0.7845, 0.8009, 
        0.9319, 0.9414, 0.9505, 0.9323, 1.0321, 0.9381, 0.8975, 1.0929, 
        1.0236, 0.9589, 1.0644, 1.0411, 1.0763, 0.9679, 1.003, 1.142, 
        1.1049, 1.2868, 1.1569, 1.1952, 1.0802, 1.2125, 1.3765, 1.263, 
        1.2507, 1.2125, 1.2207, 1.2836, 1.3352, 1.1311, 1.2321, 1.4277, 
        1.1645), w = c(898, 20566, 3011, 1364, 1520, 2376, 1923, 1934, 
        1366, 1010, 380, 421, 283, 262, 227, 173, 118, 113, 95, 69, 123, 
        70, 80, 82, 68, 83, 76, 94, 101, 97, 115, 79, 98, 84, 92, 121, 
        97, 102, 93, 92, 101, 74, 124, 64, 52, 63)), row.names = c(NA, 
        -46L), class = c("tbl_df", "tbl", "data.frame"), na.action = structure(c(`47` = 47L), class = "omit"))
    
    fit <- drm(data = df,y ~ x,fct=LL.4(), weights = w)
    plot(fit)

在此处输入图像描述

4

2 回答 2

3

1)如果我们忽略权重,那么 y = a + b * x + c/x^2 似乎拟合并且在系数中是线性的,因此很容易拟合。这似乎是向上倾斜的,所以我们从一条线开始,但后来我们需要抑制它,所以我们添加了一个倒数项。基于残差平方和,倒数二次函数的效果略好于普通倒数,所以我们改用它。

fm <- lm(y ~ x + I(1 / x^2), df)
coef(summary(fm))
##                  Estimate   Std. Error   t value     Pr(>|t|)
## (Intercept)  1.053856e+00  0.116960752  9.010341 1.849238e-11
## x            4.863077e-03  0.002718613  1.788808 8.069195e-02
## I(1/x^2)    -1.460443e+02 16.518887452 -8.841049 3.160306e-11

x 项的系数在 5% 的水平上并不显着——上表中的 p 值为 8%——所以我们可以删除它,它几乎同样适合给出只有两个参数的模型。在下图中,带有 3 个参数的 fm 拟合为实线,带有 2 个参数的 fm2 拟合为虚线。

fm2 <- lm(y ~ I(1 / x^2), df)

plot(y ~ x, df)
lines(fitted(fm) ~ x, df)
lines(fitted(fm2) ~ x, df, lty = 2)

截屏

2)另一种方法是使用两条直线。这仍然是连续的,但在过渡点有一个不可微分的点。该模型有 4 个参数,即每条线的截距和斜率。下面我们确实使用了权重。它具有基于数据外观的明显动机的优势。两条线相交处的断点可能具有重要意义,因为它是较高坡度的初始生长和较低坡度的后续生长之间的过渡点。

# starting values use lines fitted to 1st ten and last 10 points
fm_1 <- lm(y ~ x, df, subset = 1:10)
fm_2 <- lm(y ~ x, df, subset = seq(to = nrow(df), length = 10))
st <- list(a = coef(fm_1)[[1]], b = coef(fm_1)[[2]],
           c = coef(fm_2)[[1]], d = coef(fm_2)[[2]])

fm3 <- nls(y ~ pmin(a + b * x, c + d * x), df, start = st, weights = w)

# point of transition
X <- with(as.list(coef(fm3)), (a - c) / (d - b)); X
## [1] 16.38465
Y <- with(as.list(coef(fm3)), a + b * X); Y
## [1] 0.8262229

plot(y ~ x, df)
lines(fitted(fm3) ~ x, df)

截屏

于 2021-11-20T23:24:40.957 回答
2

基本思想是了解所选功能的执行方式。取一个你知道的函数(例如逻辑)并修改它。或者(甚至更好)查阅文献,看看人们在您的特定领域中使用了哪些功能。然后创建一个用户定义的模型,使用它来理解参数,定义好的起始值,然后拟合它。

她是一个快速而肮脏的用户定义函数示例(带有包增长率)。它肯定可以用 drc 类似地制作。

library("growthrates")

grow_userdefined <- function (time, parms) {
  with(as.list(parms), {
    y <- (K * y0)/(y0 + (K - y0) * exp(-mumax * time)) + shift 
    return(as.matrix(data.frame(time = time, y = y)))
  })
}

fit <- fit_growthmodel(FUN=grow_userdefined, 
  p = c(y0 = -1, K = 1, mumax = 0.1, shift = 1), 
  time = df$x, y = df$y)

plot(fit)
summary(fit)

当然可以做得更好。由于我们在开始时没有指数开始,例如可以从简单的饱和函数而不是逻辑开始,例如类似于 Monod 的东西。如前所述,首选方式是使用与应用程序域相关的功能。

改良后勤

于 2021-11-20T23:09:12.823 回答