मैं एसएमए को एलामेट्रिक डेटा का उपयोग करके फिटिंग कर रहा हूं smatr आर में पैकिंग, और मुझे 95% विश्वास अंतराल की गणना करने में कठिनाई हो रही है sma() आदेश।
पैकेज दस्तावेज़ीकरण में उदाहरण डेटा लेते हुए, मैं xy डेटा और एसएमए फिट के प्लॉट में ऊपरी और निचले 95% विश्वास लाइनों को कैसे जोड़ूंगा?
# 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% विश्वास अंतराल के लिए लाइनें कैसे जोड़ूं
धन्यवाद!
उत्तर:
उत्तर № 1 के लिए 1मुझे हाल ही में यह समस्या हुई थी और smatr पैकेज भविष्यवाणियों के आसपास आत्मविश्वास अंतराल को वापस नहीं करता है।
हालांकि, इन के माध्यम से हासिल किया जा सकता हैबूटस्ट्रैपिंग, जहाँ एक ही मॉडल पर यादृच्छिक डेटासेट बनाये जाते हैं और उन्हें फिट किया जाता है। सभी अनुमानों और ढलानों को निकालने से आप अपने पूर्वानुमानों के आसपास 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()
इसके बाद यह भूखंड देता है
