Як побудувати довірчі інтервали для основної осі, придатні за допомогою упаковки smatr в R? - r, інтервал довіри

Я підходять SMA до аллометричних даних за допомогою smatr упаковка в R, і я маю труднощі побудови 95% довірчих інтервалів, обчислених за допомогою sma() команда

Звернувши приклад даних до документації пакета, як би додати верхню і нижню 95% довірчих ліній до графіка даних xy і SMA?

# Load leaf lifetime dataset:
data(leaflife)

# Fit SMA
ft <- sma(longev~lma, data=leaflife, log="xy", method="SMA")

#plot data and fit
plot(ft, log="xy")

Як мені тепер додавати рядки для 95% довірчих інтервалів до мого сюжету?

Дякую!

Відповіді:

Я мав цю проблему нещодавно, і smatr пакет не повертає довірчі інтервали навколо прогнозів.

Проте їх можна досягти шляхомbootstrapping, де випадкові набори даних створюються знову і знову за тією ж моделлю, що підходить для них. Видобування всіх перехоплень і схилів дозволяє отримати 95% довірчі інтервали навколо ваших прогнозів.

Є, мабуть, багато способів зробити це, але я зробив це за допомогою tidyverse. При цьому використовуються пакети model, purrr, dplyr Багато інших рядків коду, ніж я б віддав перевагу тому, якщо хтось має більш елегантні рішення, дайте мені знати.

Я log10 перетворити до запуску моделі, оскільки це полегшує побудувати точки і прогнози в ggplot2.

# smatr bootstrap example

# load packages
library(smatr)
library(dplyr)
library(purrr)
library(tidyr)
library(ggplot2)

# load data
data(leaflife)

# make columns log scale
leaflife <- mutate(leaflife, log_longev = log10(longev),
log_lma = log10(lma))

# fit sma
mod <- sma(log_longev ~ log_lma, data=leaflife, method="SMA")

# plot model
plot(mod)

# create new data set of log_lma at a high resolution (200 points from min to max)
preds <- data.frame(expand.grid(log_lma = seq(min(leaflife$log_lma, na.rm = T), max(leaflife$log_lma, na.rm = T), length.out = 200), stringsAsFactors = FALSE))

# bootstrap data and get predictions
preds <- leaflife %>%
# create new bootstrapped data sets
modelr::bootstrap(n = 1000, id = "boot_num") %>%
# fit sma to every bootstrap
group_by(boot_num) %>%
mutate(., fit = map(strap, ~ sma(log_longev ~ log_lma, data=data.frame(.), method="SMA"))) %>%
ungroup() %>%
# extract intercept and slope from each fit
mutate(., intercept = map_dbl(fit, ~coef(.x)[1]),
slope = map_dbl(fit, ~coef(.x)[2])) %>%
select(., -fit) %>%
# get fitted values for each bootstrapped model
# uses the preds dataframe we made earlier
group_by(boot_num) %>%
do(data.frame(fitted = .$intercept + .$slope*preds$log_lma,
log_lma = preds$log_lma)) %>%
ungroup() %>%
# calculate the 2.5% and 97.5% quantiles at each log_lma value
group_by(., log_lma) %>%
dplyr::summarise(., conf_low = quantile(fitted, 0.025),
conf_high = quantile(fitted, 0.975)) %>%
ungroup() %>%
# add fitted value of actual unbootstrapped model
mutate(., log_longev = coef(mod)[1] + coef(mod)[2]*log_lma)

# plot with ggplot
ggplot(leaflife, aes(log_lma, log_longev)) +
geom_point() +
geom_line(data = preds) +
geom_ribbon(aes(ymin = conf_low, ymax = conf_high), alpha = 0.1, preds) +
theme_bw()

Це дає цей сюжет