2015-09-16 5 views
5

में एक बूटस्ट्रैप उत्पादन का विश्वास अंतराल मैं एक dataframe df है (नीचे देखें)प्लॉट मंझला, ggplot2

dput(df) 
    structure(list(x = c(49, 50, 51, 52, 53, 54, 55, 56, 1, 2, 3, 
    4, 5, 14, 15, 16, 17, 2, 3, 4, 5, 6, 10, 11, 3, 30, 64, 66, 67, 
    68, 69, 34, 35, 37, 39, 2, 17, 18, 99, 100, 102, 103, 67, 70, 
    72), y = c(2268.14043972082, 2147.62290922552, 2269.1387550775, 
    2247.31983098201, 1903.39138268307, 2174.78291538358, 2359.51909126411, 
    2488.39004804939, 212.851575751527, 461.398994384333, 567.150629704352, 
    781.775113821961, 918.303706148872, 1107.37695799186, 1160.80594193377, 
    1412.61328924168, 1689.48879626486, 260.737164468854, 306.72700499362, 
    283.410379620422, 366.813913489692, 387.570173754128, 388.602676983443, 
    477.858510450125, 128.198042456082, 535.519377609133, 1028.8780498564, 
    1098.54431357711, 1265.26965941035, 1129.58344809909, 820.922447928053, 
    749.343583476846, 779.678206156474, 646.575242339517, 733.953282899613, 
    461.156280127354, 906.813018662913, 798.186995701282, 831.365377249207, 
    764.519073183124, 672.076289062505, 669.879217186302, 1341.47673353751, 
    1401.44881976186, 1640.27575962036)), .Names = c("x", "y"), row.names = c(NA, 
    -45L), class = "data.frame") 

मैं एक गैर रेखीय प्रतिगमन (NLS) मेरी डाटासेट के आधार पर बनाया है।

nls1 <- nls(y~A*(x^B)*(exp(k*x)), 
      data = df, 
      start = list(A = 1000, B = 0.170, k = -0.00295), algorithm = "port") 

मैंने पैरामीटर (ए, बी और के) के कई सेट प्राप्त करने के लिए इस फ़ंक्शन के लिए बूटस्ट्रैप की गणना की।

library(nlstools) 
Boo <- nlsBoot(nls1, niter = 200) 

मैं अब मंझला वक्र के साथ ही ऊपरी और निचले विश्वास अंतराल घटता एक ggplot2 में एक साथ बूटस्ट्रैप वस्तु से गणना की साजिश करना चाहते हैं। प्रत्येक वक्र के पैरामीटर (ए, बी और के) Boo_Gamma$bootCI में निहित है। क्या कोई मेरी मदद कर सकता है? अग्रिम में धन्यवाद।

+2

में आपका स्वागत है अतः करने के लिए और एक _perfectly formatted_ पहला सवाल के लिए +1000! क्या 'nlstools' में एक समकक्ष साजिश का निर्माण है जिसे आप ggplot2 पर "पोर्ट" करने का प्रयास कर रहे हैं? – hrbrmstr

+0

वास्तव में दुर्भाग्य से वास्तव में नहीं। मुझे मध्यवर्ती वक्र के चारों ओर आत्मविश्वास अंतराल बैंड की साजिश के लिए कोई फंक्शन नहीं मिला। – SimonB

उत्तर

2

AFAIK, पैकेज nlstools केवल बूटस्ट्रैप पैरामीटर अनुमान नहीं भविष्यवाणी मूल्यों रिटर्न ...

इसलिए, यहाँ है एक त्वरित समाधान, मैन्युअल रूप से चलाया पैरामीटर अनुमान का उपयोग कर भविष्यवाणियों गणना करने के लिए और उसके बाद आकड़े recomputing भविष्यवाणियों का, क्योंकि यहां मॉडल गैर-रैखिक है। यह सबसे सुंदर नहीं है, लेकिन यह क्या करना चाहिए यह :)

# Matrix with the bootstrapped parameter estimates 
Theta_mat <- Boo$coefboot 

# Model 
fun <- function(x, theta) theta["A"] * (x^theta["B"]) * (exp(theta["k"] * x)) 

# Points where to evaluate the model 
x_eval <- seq(min(df$x), max(df$x), length.out = 100) 

# Matrix with the predictions 
Pred_mat <- apply(Theta_mat, 1, function(theta) fun(x_eval, theta)) 

# Pack the estimates for plotting 
Estims_plot <- cbind(
    x = x_eval, 
    as.data.frame(t(apply(Pred_mat, 1, function(y_est) c(
     median_est = median(y_est), 
     ci_lower_est = quantile(y_est, probs = 0.025, names = FALSE), 
     ci_upper_est = quantile(y_est, probs = 0.975, names = FALSE) 
    )))) 
) 

library(ggplot2) 
ggplot(data = Estims_plot, aes(x = x, y = median_est, ymin = ci_lower_est, ymax = ci_upper_est)) + 
    geom_ribbon(alpha = 0.7, fill = "grey") + 
    geom_line(size = rel(1.5), colour = "black") + 
    geom_point(data = df, aes(x = x, y = y), size = rel(4), colour = "red", inherit.aes = FALSE) + 
    theme_bw() + labs(title = "Bootstrap results\n", x = "x", y = "y") 
ggsave("bootpstrap_results.pdf", height = 5, width = 9) 

Bootstrap results

+0

कूल। बहुत बहुत धन्यवाद। – SimonB

+0

इसलिए, आप बस चाल करने के लिए 'Boo_Gamma $ bootCI' ऑब्जेक्ट को दिए गए मध्य और सीआई के पैरामीटर का उपयोग नहीं कर सकते हैं? – SimonB

संबंधित मुद्दे