ggplot, um Konfidenzintervalle von Bootstrapping-Kurvenanpassung zu zeigen

Question

Dank der Hilfe, die ich bekam here, konnte ich ein Spaghetti-Plot Kurvenanpassungen mit Bootstrapping bekommen. Ich versuche, Vertrauensbänder aus diesen angepassten Modellen abzuleiten. Ich habe wie

quants <- apply(fitted_boot, 1, quantile, c(0.025, 0.5, 0.975)) 

arbeiten mit folgendem kein Glück immer etwas hatte:

library(dplyr) 
library(broom) 
library(ggplot2) 

xdata <- c(-35.98, -34.74, -33.46, -32.04, -30.86, -29.64, -28.50, -27.29, -26.00, 
      -24.77, -23.57, -22.21, -21.19, -20.16, -18.77, -17.57, -16.47, -15.35, 
      -14.40, -13.09, -11.90, -10.47, -9.95,-8.90,-7.77,-6.80, -5.99, 
      -5.17, -4.21, -3.06, -2.29, -1.04) 
ydata <- c(-4.425, -4.134, -5.145, -5.411, -6.711, -7.725, -8.087, -9.059, -10.657, 
      -11.734, NA, -12.803, -12.906, -12.460, -12.128, -11.667, -10.947, -10.294, 
      -9.185, -8.620, -8.025, -7.493, -6.713, -6.503, -6.316, -5.662, -5.734, -4.984, 
      -4.723, -4.753, -4.503, -4.200) 

data <- data.frame(xdata,ydata) 
x_range <- seq(min(xdata), max(xdata), length.out = 1000) 

fitted_boot <- data %>% 
    bootstrap(100) %>% 
    do({ 
    m <- nls(ydata ~ A*cos(2*pi*((xdata-x_0)/z))+M, ., start=list(A=4,M=-7,x_0=-10,z=30)) 
    f <- predict(m, newdata = list(xdata = x_range)) 
    data.frame(xdata = x_range, .fitted = f) 
    }) 

ggplot(data, aes(xdata, ydata)) + 
    geom_line(aes(y=.fitted, group=replicate), fitted_boot, alpha=.1, color="blue") + 
    geom_point(size=3) + 
    theme_bw() 

Ich dachte, vielleicht geom_ribbon() wäre ein schöner Weg zu gehen, aber ich Ich weiß einfach nicht, wo ich von hier aus hingehen soll.

Danke an Axeman für die Hilfe auf dem anderen Beitrag!

Answer 1

Ein Ansatz wäre es, das Konfidenzintervall bei jedem x-Wert zu berechnen und dann einfach zu plotten. Hier verwende ich den ersten Wert außerhalb des 2.5. Perzentils und des 97.5. Perzentils, obwohl Sie den Code nach Bedarf anpassen könnten.

Zuerst wechsle ich zu group_by die xdata Standorte (anstelle von Replikaten). Dann, ich arrange durch die .fitted Werte, so dass ich slice die Werte, die ich will (die erste außerhalb der Perzentil-Cutoffs). Schließlich markiere ich sie, mit welchem ​​Band ich komme (sie gehen immer tiefer als oben, weil wir sortiert haben).

forConfInt <- 
    fitted_boot %>% 
    ungroup() %>% 
    group_by(xdata) %>% 
    arrange(.fitted) %>% 
    slice(c(floor(0.025 * n()) 
      , ceiling(0.975 * n()))) %>% 
    mutate(range = c("lower", "upper")) 

Dies gibt:

replicate  xdata .fitted range 
     <int>  <dbl>  <dbl> <chr> 
1   9 -35.98000 -4.927462 lower 
2   94 -35.98000 -4.249348 upper 
3   9 -35.94503 -4.927248 lower 
4   94 -35.94503 -4.257776 upper 
5   9 -35.91005 -4.927228 lower 
6   94 -35.91005 -4.266334 upper 
7   9 -35.87508 -4.927401 lower 
8   94 -35.87508 -4.275020 upper 
9   9 -35.84010 -4.927766 lower 
10  94 -35.84010 -4.283836 upper 
# ... with 1,990 more rows 

Und wir können dann nur eine zusätzliche Leitung zum ggplot Anruf hinzufügen:

ggplot(data, aes(xdata, ydata)) + 
    geom_line(aes(y=.fitted, group=replicate), fitted_boot, alpha=.1, color="blue") + 
    # Added confidence interval: 
    geom_line(aes(y=.fitted, group=range), forConfInt, color="red") + 
    geom_point(size=3) + 
    theme_bw() 

Gibt dieses Grundstück: