Como traçar intervalos de confiança para um ajuste do eixo principal usando o pacote smatr em R? - r, intervalo de confiança

Estou ajustando o SMA aos dados alométricos usando o smatr embalagem em R, e estou tendo dificuldade em traçar os intervalos de confiança de 95% calculados pelo sma() comando.

Tomando os dados de exemplo na documentação do pacote, como eu adicionaria as linhas de confiança de 95% superior e inferior ao gráfico de dados xy e ajuste 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")

Como adiciono agora linhas para os intervalos de confiança de 95% ao meu gráfico?

Obrigado!

Respostas:

Eu tive esse problema recentemente e o pacote smatr não retorna intervalos de confiança em torno das previsões.

No entanto, estes poderiam ser alcançadosbootstrapping, onde conjuntos de dados aleatórios são feitos sobre e sobre o mesmo modelo ajustado a eles. A extração de todos os interceptos e inclinações permite obter intervalos de confiança de 95% em torno de suas previsões.

Há provavelmente muitas maneiras de fazer isso, mas eu fiz isso usando o tidyverse. Isso usa os pacotes model, purrr, dplyr etc Muito mais linhas de código do que eu preferiria, se alguém tiver mais elegantes soluções, me avise.

Eu log10 transforma antes da execução do modelo, pois facilita a plotagem dos pontos e previsões em 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()

Isso então dá este enredo