如何绘制自举置信区间
how to plot bootstraping confidence intervals
我 运行 一些 bootstrap 置信区间，我想用均值绘制置信区间。像这样：
这是我的模型。如您所见，DISTANCE 和 POS 是因素。
lmm1 <- lmer((total) ~ DISTANCE+POS + (1|NO_UNIT),data=TURN)
TURN$POS<-as.factor(TURN$POS)#Change position and distance to factors
TURN$DISTANCE<-as.factor(TURN$DISTANCE)
TURN$NO_UNIT <- as.factor(TURN$NO_UNIT)
这是我正在使用的功能：
mySumm <- function(.) { s <- sigma(.)
c(beta =getME(., "beta"), sigma = s, sig01 = unname(s * getME(., "theta"))) }
# run bootstrap analysis for calculation of confidence intervals of parameter estimates
mod_lmm1_boot <- bootMer(lmm1,mySumm, nsim=300)
# confidence interval forcoefficient DISTANCE1
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=2)
# confidence interval forcoefficient DISTANCE2
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=3)
# confidence interval forcoefficient DISTANCE3
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=4)
# confidence interval forcoefficient DISTANCE4
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=5)
# confidence interval forcoefficient DISTANCECONTROL
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=6)
# confidence interval forcoefficient POS2
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=7)
# confidence interval forcoefficient POS3
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=8)
# confidence interval forcoefficient POS4
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=9)
# confidence interval forcoefficient POS5
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=10)#I want to plot the confidence intervals and means corresponding to these indexes!
这是TURN数据框
TREAT NO_UNIT POS DISTANCE total
NBNA 1 1 Control 2525837593
NBNA 1 2 Control 2127040915
NBNA 1 3 Control 1387070744
NBNA 1 4 Control 1458541776
NBNA 1 5 Control 1858881414
NBNA 2 1 Control NA
NBNA 2 2 Control 1481932857
NBNA 2 3 Control 2037767853
NBNA 2 4 Control 2056201434
NBNA 2 5 Control 2265049369
NBNA 3 1 Control 1474510932
NBNA 3 2 Control 1801028575
NBNA 3 3 Control 1542852992
NBNA 3 4 Control 2210304532
NBNA 3 5 Control 2557068715
NBA1 4 1 0 640224197.7
NBA1 4 2 1 704589654.3
NBA1 4 3 2 558153711.5
NBA1 4 4 3 2263088969
NBA1 4 5 4 1779212490
NBA1 5 1 0 1483773056
NBA1 5 2 1 1539240152
NBA1 5 3 2 1987871365
NBA1 5 4 3 2555767964
NBA1 5 5 4 2040386118
NBA1 6 1 0 1855829526
NBA1 6 2 1 1868973099
NBA1 6 3 2 1345878086
NBA1 6 4 3 1651096675
NBA1 6 5 4 1513168820
NBA5 7 1 4 1832017482
NBA5 7 2 3 1858590718
NBA5 7 3 2 1445652450
NBA5 7 4 1 1544741442
NBA5 7 5 0 1330161829
NBA5 8 1 4 1896550338
NBA5 8 2 3 2016638692
NBA5 8 3 2 1333723238
NBA5 8 4 1 1570307025
NBA5 8 5 0 1666068148
NBA5 9 1 4 NA
NBA5 9 2 3 1990898325
NBA5 9 3 2 1675553680
NBA5 9 4 1 1556562879
NBA5 9 5 0 1910142673
BNA 10 1 Control 370990618.1
BNA 10 2 Control 484424075.2
BNA 10 3 Control 659926517.8
BNA 10 4 Control NA
BNA 10 5 Control 1532572993
BNA 11 1 Control 590917346
BNA 11 2 Control 1826109624
BNA 11 3 Control 318371884.5
BNA 11 4 Control 3046708581
BNA 11 5 Control NA
BNA 12 1 Control 755992256.9
BNA 12 2 Control 457416788.1
BNA 12 3 Control 874742750.4
BNA 12 4 Control 67841374.52
BNA 12 5 Control 2933480357
BA1 13 1 0 12067695.33
BA1 13 2 1 10083668.91
BA1 13 3 2 6416950.819
BA1 13 4 3 7398691.286
BA1 13 5 4 10860980.63
BA1 14 1 0 11892423.38
BA1 14 2 1 27004799.05
BA1 14 3 2 597357273.2
BA1 14 4 3 1572190656
BA1 14 5 4 1249666803
BA1 15 1 0 38998930.08
BA1 15 2 1 279361097.3
BA1 15 3 2 350236872
BA1 15 4 3 931806452.5
BA1 15 5 4 NA
BA5 16 1 4 623714889
BA5 16 2 3 481547462
BA5 16 3 2 992879231.3
BA5 16 4 1 2287090423
BA5 16 5 0 1742484997
BA5 17 1 4 786425214.1
BA5 17 2 3 899998604.5
BA5 17 3 2 1244515257
BA5 17 4 1 2432196221
BA5 17 5 0 386404680.3
BA5 18 1 4 597970711
BA5 18 2 3 781618489.7
BA5 18 3 2 2024931390
BA5 18 4 1 1663706249
BA5 18 5 0 1231669010
首先，要将 mod_lmm1_boot 对象转换为数据框，您可以对 "boot" 对象使用 tidy 方法（来自我的 broom 包）：
library(broom)
tidy(mod_lmm1_boot)
给出输出：
term statistic bias std.error
1 beta1 511476574 -2340765.59 250202968
2 beta2 264511629 8063563.41 239403518
3 beta3 104454362 5408454.64 237085206
4 beta4 391743711 12945670.41 231843390
5 beta5 188803173 -6839936.62 243065434
6 beta6 479700475 6934270.75 308630904
7 beta7 171444397 -87209.96 44159725
8 sigma 566047522 -10557609.04 50260359
9 sig01 476939097 10975306.03 121196856
您还可以计算您列出的置信区间，例如：
ci <- do.call(rbind, lapply(1:9, function(i) {
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=i)$percent
}))
它给出了一个包含 5 列的矩阵，其中第 4 和第 5 是您的百分位数置信区间。
你用 ggplot2 标记了你的问题，所以这里是在 ggplot2 中绘制结果间隔的代码：
library(broom)
library(ggplot2)
td <- tidy(mod_lmm1_boot)
td$conf.low <- ci[, 4]
td$conf.high <- ci[, 5]
ggplot(td, aes(term, statistic)) +
geom_point() +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high))
我 运行 一些 bootstrap 置信区间，我想用均值绘制置信区间。像这样：
这是我的模型。如您所见，DISTANCE 和 POS 是因素。
lmm1 <- lmer((total) ~ DISTANCE+POS + (1|NO_UNIT),data=TURN)
TURN$POS<-as.factor(TURN$POS)#Change position and distance to factors
TURN$DISTANCE<-as.factor(TURN$DISTANCE)
TURN$NO_UNIT <- as.factor(TURN$NO_UNIT)
这是我正在使用的功能：
mySumm <- function(.) { s <- sigma(.)
c(beta =getME(., "beta"), sigma = s, sig01 = unname(s * getME(., "theta"))) }
# run bootstrap analysis for calculation of confidence intervals of parameter estimates
mod_lmm1_boot <- bootMer(lmm1,mySumm, nsim=300)
# confidence interval forcoefficient DISTANCE1
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=2)
# confidence interval forcoefficient DISTANCE2
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=3)
# confidence interval forcoefficient DISTANCE3
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=4)
# confidence interval forcoefficient DISTANCE4
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=5)
# confidence interval forcoefficient DISTANCECONTROL
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=6)
# confidence interval forcoefficient POS2
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=7)
# confidence interval forcoefficient POS3
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=8)
# confidence interval forcoefficient POS4
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=9)
# confidence interval forcoefficient POS5
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=10)#I want to plot the confidence intervals and means corresponding to these indexes!
这是TURN数据框
TREAT NO_UNIT POS DISTANCE total
NBNA 1 1 Control 2525837593
NBNA 1 2 Control 2127040915
NBNA 1 3 Control 1387070744
NBNA 1 4 Control 1458541776
NBNA 1 5 Control 1858881414
NBNA 2 1 Control NA
NBNA 2 2 Control 1481932857
NBNA 2 3 Control 2037767853
NBNA 2 4 Control 2056201434
NBNA 2 5 Control 2265049369
NBNA 3 1 Control 1474510932
NBNA 3 2 Control 1801028575
NBNA 3 3 Control 1542852992
NBNA 3 4 Control 2210304532
NBNA 3 5 Control 2557068715
NBA1 4 1 0 640224197.7
NBA1 4 2 1 704589654.3
NBA1 4 3 2 558153711.5
NBA1 4 4 3 2263088969
NBA1 4 5 4 1779212490
NBA1 5 1 0 1483773056
NBA1 5 2 1 1539240152
NBA1 5 3 2 1987871365
NBA1 5 4 3 2555767964
NBA1 5 5 4 2040386118
NBA1 6 1 0 1855829526
NBA1 6 2 1 1868973099
NBA1 6 3 2 1345878086
NBA1 6 4 3 1651096675
NBA1 6 5 4 1513168820
NBA5 7 1 4 1832017482
NBA5 7 2 3 1858590718
NBA5 7 3 2 1445652450
NBA5 7 4 1 1544741442
NBA5 7 5 0 1330161829
NBA5 8 1 4 1896550338
NBA5 8 2 3 2016638692
NBA5 8 3 2 1333723238
NBA5 8 4 1 1570307025
NBA5 8 5 0 1666068148
NBA5 9 1 4 NA
NBA5 9 2 3 1990898325
NBA5 9 3 2 1675553680
NBA5 9 4 1 1556562879
NBA5 9 5 0 1910142673
BNA 10 1 Control 370990618.1
BNA 10 2 Control 484424075.2
BNA 10 3 Control 659926517.8
BNA 10 4 Control NA
BNA 10 5 Control 1532572993
BNA 11 1 Control 590917346
BNA 11 2 Control 1826109624
BNA 11 3 Control 318371884.5
BNA 11 4 Control 3046708581
BNA 11 5 Control NA
BNA 12 1 Control 755992256.9
BNA 12 2 Control 457416788.1
BNA 12 3 Control 874742750.4
BNA 12 4 Control 67841374.52
BNA 12 5 Control 2933480357
BA1 13 1 0 12067695.33
BA1 13 2 1 10083668.91
BA1 13 3 2 6416950.819
BA1 13 4 3 7398691.286
BA1 13 5 4 10860980.63
BA1 14 1 0 11892423.38
BA1 14 2 1 27004799.05
BA1 14 3 2 597357273.2
BA1 14 4 3 1572190656
BA1 14 5 4 1249666803
BA1 15 1 0 38998930.08
BA1 15 2 1 279361097.3
BA1 15 3 2 350236872
BA1 15 4 3 931806452.5
BA1 15 5 4 NA
BA5 16 1 4 623714889
BA5 16 2 3 481547462
BA5 16 3 2 992879231.3
BA5 16 4 1 2287090423
BA5 16 5 0 1742484997
BA5 17 1 4 786425214.1
BA5 17 2 3 899998604.5
BA5 17 3 2 1244515257
BA5 17 4 1 2432196221
BA5 17 5 0 386404680.3
BA5 18 1 4 597970711
BA5 18 2 3 781618489.7
BA5 18 3 2 2024931390
BA5 18 4 1 1663706249
BA5 18 5 0 1231669010
首先，要将 mod_lmm1_boot 对象转换为数据框，您可以对 "boot" 对象使用 tidy 方法（来自我的 broom 包）：
library(broom)
tidy(mod_lmm1_boot)
给出输出：
term statistic bias std.error
1 beta1 511476574 -2340765.59 250202968
2 beta2 264511629 8063563.41 239403518
3 beta3 104454362 5408454.64 237085206
4 beta4 391743711 12945670.41 231843390
5 beta5 188803173 -6839936.62 243065434
6 beta6 479700475 6934270.75 308630904
7 beta7 171444397 -87209.96 44159725
8 sigma 566047522 -10557609.04 50260359
9 sig01 476939097 10975306.03 121196856
您还可以计算您列出的置信区间，例如：
ci <- do.call(rbind, lapply(1:9, function(i) {
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=i)$percent
}))
它给出了一个包含 5 列的矩阵，其中第 4 和第 5 是您的百分位数置信区间。
你用 ggplot2 标记了你的问题，所以这里是在 ggplot2 中绘制结果间隔的代码：
library(broom)
library(ggplot2)
td <- tidy(mod_lmm1_boot)
td$conf.low <- ci[, 4]
td$conf.high <- ci[, 5]
ggplot(td, aes(term, statistic)) +
geom_point() +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high))