当前位置 : 主页 > 网页制作 > Bootstarp >

在ggplot2中绘制bootstrap输出的中位数,置信区间

来源:互联网 收集:自由互联 发布时间:2021-06-12
我有一个数据帧df(见下文) 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 =
我有一个数据帧df(见下文)

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")

然后,我为此函数计算了一个引导程序,以获取多组参数(A,B和k).

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

我现在想要绘制中值曲线以及在一个ggplot2中从引导对象一起计算的上下置信区间曲线.每条曲线的参数(A,B和K)包含在Boo_Gamma $bootCI中.有人能帮帮我吗?提前致谢.

AFAIK,包nlstools只返回bootstrapped参数估计值,而不是预测值……

因此,这是一个快速的解决方案,手动使用自举参数估计来计算预测,然后重新计算预测中的统计数据,因为这里的模型是非线性的.它不是最优雅的,但应该这样做:)

# 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

网友评论