Почему есть разница между do (lm) и geom_smooth (method = lm)?

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

df_calibration=structure(list(dilution = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 
0.8, 0.9, 1, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1), 
    area = c(1000, 2000, 3000, 4000, 5000, 6000, 7000, 7800, 
    8200, 8500, 1200, 2200, 3200, 4200, 5200, 6200, 7200, 8000, 
    8400, 8700), substance = c("A", "A", "A", "A", "A", "A", 
    "A", "A", "A", "A", "b", "b", "b", "b", "b", "b", "b", "b", 
    "b", "b")), row.names = c(NA, 20L), class = "data.frame")

df_samples=structure(list(area = c(1100, 1800, 2500, 3200, 3900, 1300, 2000, 
2700, 3400, 4100), substance = c("A", "A", "A", "A", "A", "b", 
"b", "b", "b", "b")), row.names = c(NA, 10L), class = "data.frame")

Чтобы вычислить фактические разведения по измеренным образцам, я беру параметры, полученные в результате этой подгонки:

df_fits=df_calibration %>% group_by(substance) %>% 
  do(fit = lm(area ~ poly(dilution,2), data = .))%>%
  tidy(fit) %>% 
  select(substance, term, estimate) %>% 
  spread(term, estimate)

df_fits=df_fits %>% rename(a=`poly(dilution, 2)2`,b=`poly(dilution, 2)1`,c=`(Intercept)`)

#join parameters with sample data
df_samples=left_join(df_samples,df_fits)

и эта формула формула для вычисления 

#calculate with general solution for polynomial 2nd order
df_samples$dilution_calc=
  (df_samples$b*(-1)+sqrt(df_samples$b^2-(4*df_samples$a*(df_samples$c-df_samples$area))))/(2*df_samples$a)

Однако, когда я рисую это сейчас, я замечаю кое-что очень странное. Рассчитанные значения x (разведения) не попадают на кривую из stat_smooth(). Дополнительной пунктирной линией нанесены параметры из уравнения на графике (которые соответствуют числам во фрейме данных) для вещества «А». Значит, мои расчеты должны быть правильными (или нет?) Почему есть разница? Что я делаю неправильно? Как я могу получить параметры из подгонки, выполненной stat_smooth()?

my.formula=y ~ poly(x,2)
ggplot(df_calibration, aes(x = dilution, y = area)) +
  stat_smooth(method = "lm", se=FALSE, formula = my.formula) +

  stat_function(fun=function(x){5250+(7980*x)+(-905*x^2)},      
              inherit.aes = F,linetype="dotted")+

  stat_poly_eq(formula = my.formula, 
               aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               parse = TRUE) +         
  geom_point(shape=17)+
  geom_point(data=df_samples,
           aes(x=dilution_calc,y=area),
           shape=1,color="red")+
  facet_wrap(~substance,scales = "free")

сюжет со странным поведением

Любое предложение будет высоко оценено :-)


r ggplot2 model-fitting curve-fitting data-fitting
person TobiO    schedule 31.03.2019    source источник
comment
Я не могу понять, какое отношение имеет четвертичная формула к предсказаниям на основе полинома второго порядка. Поэтому вам, вероятно, следует более подробно описать, какую логику вы используете. Если вы используете коэффициенты из poly(), тогда вам а) необходимо включить точку пересечения и б) нужно помнить, что член второго порядка был сделан ортогональным члену первого порядка. И, несмотря на выполнение library(tidyverse), я все еще получаю ошибку в tidy (., Fit): не удалось найти функцию tidy   -  person IRTFM    schedule 31.03.2019

Ответы (1)


arrow_upward
2
arrow_downward

По умолчанию poly вычисляет ортогональные многочлены. Вы можете отключить ортогонализацию с помощью аргумента raw=TRUE.

Обратите внимание, что формула имеет два вида: один раз с исходными именами переменных в соответствии с регрессиями, а затем в stat_smooth с использованием общих имен переменных x и y. Но в остальном это должна быть та же формула с raw=TRUE.

library("tidyverse")

# Define/import your data here....

df_fits <- df_calibration %>%
  group_by(substance) %>%
  do(fit = lm(area ~ poly(dilution, 2, raw = TRUE), data = .)) %>%
  broom::tidy(fit) %>%
  select(substance, term, estimate) %>%
  spread(term, estimate) %>%
  # It is simpler to rename the coefficients here
  setNames(c("substance", "c", "b", "a"))

# join parameters with sample data
df_samples <- left_join(df_samples, df_fits)

# calculate with general solution for polynomial 2nd order
df_samples <- df_samples %>%
  mutate(dilution_calc = (b * (-1) + sqrt(b^2 - (4 * a * (c - area)))) / (2 * a))

my.formula <- y ~ poly(x, 2, raw = TRUE)

df_calibration %>%
  ggplot(aes(x = dilution, y = area)) +
  stat_smooth(method = "lm", se = FALSE, formula = my.formula) +
  geom_point(shape = 17) +
  geom_point(
    data = df_samples,
    aes(x = dilution_calc, y = area),
    shape = 1, color = "red"
  ) +
  facet_wrap(~substance, scales = "free")

Создано 31 марта 2019 г. пакетом REPEX (v0.2.1)

person dipetkov    schedule 31.03.2019