Как я могу оценить переменную x из переменной y в R?

Question

Вот мои данные:

# A tibble: 8 x 3
    CFU strain diltn
  <dbl> <chr>  <dbl>
1 159   aM12    8748
2 124.  aM12    2916
3  76.5 aM12     972
4  22   aM12     324
5  16.5 aM12     108
6  17   aM12      36
7  22.5 aM12      12
8  17.5 aM12       4

Это может показаться простым вопросом, но я в основном использовал R для получения основных сводок данных и их построения графиков (используя dplyr и ggplot).

Я могу построить график:

ggplot(data=data, aes(x=diltn, y=CFU))+
  geom_point()+
  geom_line()+
  scale_x_log10()

Я хотел бы оценить, что такое "diltn" (переменная x), яполучит 77 "КОЕ" (переменная у).

Я справился с этим в excel и проиллюстрировал это следующим образом, чтобы проиллюстрировать, чего я хотел бы достичь:

Answer 1

Это на самом деле гораздо более сложный вопрос (в общем), чем кажется. Дело не в том, что это невозможно сделать (есть много вариантов), но это сильно зависит от того, как ведут себя ваши данные. Например, предположим, что значение y интереса равно 20 вместо 77. Любое значение diltn между 4 и 324 теперь является «разумным» ответом.

Чтобы обойти эту проблему, мы используем статистические модели. Если я правильно угадываю, и вы работаете с моделью доза-ответ (или чем-то похожим - например, я использовал их со стандартными кривыми в анализах), вы можете проверить drm() в пакете drc, который может соответствовать этим кривым соответственно.

Что-то вроде:

mod <- drm(CFU ~ diltn, data = data, fct = LL.4())
plot(mod)

Функция ED затем используется для извлечения соответствующих данных. Я работаю со стандартными кривыми и считаю следующие настройки полезными, но вам могут понадобиться разные в зависимости от того, как работают ваши данные.

ED(mod, 77, bound = FALSE, type = 'absolute')
# Estimated effective doses
# 
#        Estimate Std. Error
# e:1:77  1103.69     176.31

Прошло некоторое время с тех пор, как я прочитал на нем виньетки, поэтому вам, вероятно, нужно будет немного почитать, чтобы убедиться, что вы получите правильный результат.

Answer 2

Основываясь на графике рассеяния, мы, вероятно, можем подогнать линию нелинейной регрессии к набору данных.Предполагая, что ваш набор данных называется dat.Мы можем использовать функцию nls, чтобы соответствовать регрессионной модели.Обратите внимание, что требуются некоторые усилия и размышления, чтобы найти возможное уравнение и начальные значения.В этом случае уравнение равно CFU ~ a * diltn/(b + diltn), а начальные значения для a и b равны 100 и 1000 соответственно.

library(tidyverse)

fit <- nls(formula = CFU ~ a * diltn/(b + diltn), 
           start = list(a = 100, b = 1000), data = dat)

summary(fit)

# Formula: CFU ~ a * diltn/(b + diltn)
# 
# Parameters:
#   Estimate Std. Error t value Pr(>|t|)    
# a   187.32      21.25   8.814 0.000118 ***
# b  1514.27     517.50   2.926 0.026420 *  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 13.17 on 6 degrees of freedom
# 
# Number of iterations to convergence: 4 
# Achieved convergence tolerance: 3.555e-06

Для визуальной проверки соответствия модели мыможно сначала создать второй кадр данных с diltn от 1 до 9000. Затем мы можем использовать функцию predict для прогнозирования КОЕ на основе diltn и модели fit.

dat2 <- data_frame(diltn = 1:9000) %>% 
  mutate(Pred = predict(fit, .))

ggplot(data = dat, aes(x = diltn, y = CFU))+
  geom_point() +
  geom_line(data = dat2, aes(x = diltn, y = Pred), color = "red")

Модель выглядит хорошо для меня.

Наконец, мы можем отфильтровать значения Pred, чтобы найти возможные значения для diltn.В этом случае, я думаю, что 1057 может быть возможным ответом.

dat2 %>% filter(Pred > 76.9, Pred < 77.1)

# # A tibble: 5 x 2
#   diltn  Pred
#   <int> <dbl>
# 1  1055  76.9
# 2  1056  77.0
# 3  1057  77.0
# 4  1058  77.0
# 5  1059  77.1

Или, поскольку мы установили модель нелинейной регрессии, и мы знаем, что подобранные параметры a и b, мыможно установить CFU = 77 и вычислить diltn.Мой расчет показывает diltn это 1056.914.