Я собирался написать комментарий, но я выдвигал ограничение на число символов.
Я не знаком с пакетом inflection
, поэтому я не один, чтобы судить, является ли 2020-02-08
истинным перегибом , Однако я скажу, что это трудно ответить с помощью R, потому что R не обязательно хорош в вычислении производных. Если бы у вас было оценочное уравнение линии - тогда вы могли бы потенциально использовать его для построения первой и второй производных. Вычисление грубой дельты с помощью разницы в (Y_n+1-Y_n)/(X_n+1-X_n)
никогда не является оптимальным, поскольку теоретически производный представляет собой дельту двух точек, бесконечно близко расположенных друг к другу. Вы принципиально не можете получить отличную оценку производной. Вы даже можете увидеть это, потому что вы вынуждены изменить эту оценку на n
или n+1
. Кроме того, можно ожидать, что точка перегиба x_0
будет локальным минимумом / максимумом в первой производной и равна нулю во второй производной. Так что я не думаю, что ваш второй заговор помогает. Но это может быть связано с подсчетом дельты.
Что бы я сделал, это сначала подгонил ваши данные к какой-то модели. В этом примере я собираюсь использовать пакет dr4pl
, чтобы смоделировать ваши данные для модели logisti c с 4 параметрами. Поскольку функция четырехпараметрической модели хорошо известна, я могу написать, какой должна быть первая и вторая производные функции, затем построить эти значения, используя stat_function
в пакете ggplot2
.
library(ggplot2)
library(dr4pl)
df<-structure(list(date = structure(c(18277, 18278, 18279, 18280,
18281, 18282, 18283, 18284, 18285, 18286, 18287, 18288, 18289,
18290, 18291, 18292, 18293, 18294, 18295, 18296, 18297, 18298,
18299, 18300, 18301, 18302, 18303, 18304, 18305, 18306, 18307),
class = "Date"),
cases = c(45, 62, 121, 198, 258, 363, 425,
495, 572, 618, 698, 1590, 1905, 2261, 2639, 3125, 4109, 5142,
6384, 8351, 10117, 11618, 13603, 14982, 16903, 18454, 19558,
20630, 21960, 22961, 23621)),
class = "data.frame", row.names = c(NA, -31L))
xlb_0<- structure(c(18281, 18285, 18289, 18293,
18297, 18301, 18305,
18309), class = "Date")
df$dat_as_num <- as.numeric(df$date)
dr4pl_obj <- dr4pl(cases~dat_as_num, data = df, init.parm = c(30000, 18300, 2, 0))
#first derivative derivation
d1_dr4pl <- function(x, theta, scale = F){
if (any(is.na(theta))) {
stop("One of the parameter values is NA.")
}
if (theta[2] <= 0) {
stop("An IC50 estimate should always be positive.")
}
f <- -theta[3]*((theta[4]-theta[1])/((1+(x/theta[2])^theta[3])^2))*((x/theta[2])^(theta[3]-1))
if(scale) {
f <- scales::rescale(x = f, to = c(theta[4],theta[1]))
}
return(f)
}
#Second derivative derivation
d2_dr4pl <- function(x, theta, scale = F){
if (any(is.na(theta))) {
stop("One of the parameter values is NA.")
}
if (theta[2] <= 0) {
stop("An IC50 estimate should always be positive.")
}
f <- 2*((theta[3]*(x/theta[2])^(theta[3]-1))^2)*((theta[4]-theta[1])/((1+(x/theta[2])^(theta[3]))^3))-theta[3]*(theta[3]-1)*((x/theta[2])^(theta[3]-2))*((theta[4]-theta[1])/((1+(x/theta[2])^theta[3])^2))
if(scale) {
f <- scales::rescale(x = f, to = c(theta[4],theta[1]))
f <- f - f[1]
}
return(f)
}
ggplot(df, aes(x = dat_as_num)) +
geom_hline(yintercept = 0) +
[![enter image description here][1]][1]geom_point(aes(y = cases), color = "grey", alpha = .6, size = 5) +
stat_function(fun = d1_dr4pl, args = list(theta = dr4pl_obj$parameters, scale = T), color = "red") +
stat_function(fun = d2_dr4pl, args = list(theta = dr4pl_obj$parameters, scale = T), color = "blue") +
stat_function(fun = dr4pl::MeanResponse, args = list(theta = dr4pl_obj$parameters), color = "gold") +
geom_vline(xintercept = dr4pl_obj$parameters[2], linetype = "dotted") +
theme_classic()
Как вы можете видеть, точка перегиба, которая является значением IC50 (тета 2) модели logisti c с 4 параметрами, хорошо выстраивается, когда мы приближаемся к этому way.
summary(dr4pl_obj)
#$call
#dr4pl.formula(formula = cases ~ dat_as_num, data = df, init.parm = c(30000, 18300, 2, 0))
#
#$coefficients
# Estimate StdErr 2.5 % 97.5 %
#Upper limit 25750.61451 4.301008e-05 25750.59681 25750.63221
#Log10(IC50) 18298.75347 4.301008e-09 18298.67889 18298.82806
#Slope 5154.35449 4.301008e-05 5154.33678 5154.37219
#Lower limit 58.48732 4.301008e-05 58.46962 58.50503
#
#attr(,"class")
#[1] "summary.dr4pl"
Кроме того, используя dr4pl
, он говорит, что значение IC50 примерно 18298.8
, что позднее 2020-02-06
. Не слишком далеко от значения inflection
. Я уверен, что может быть лучше использовать модель, чем 4pl, но я знал, что это была единственная модель, которую я мог бы написать первое и второе производные для ответа на этот вопрос.
Я уверен другие языки кодирования более специализированы, когда дело доходит до производных, и могут даже рассчитать их для вас, если вы начнете с начальной функции. Я думаю, что один из этих языков - mathematica.
A disclaimer , я закончил масштабирование первого и второго производных так, чтобы они могли быть построены вместе. Их фактические значения намного больше, чем показано здесь.