2016-12-07 50 views
1

我正在運行一些自舉置信區間,我想用均值繪製置信區間。類似這樣的:enter image description here如何繪製自舉置信區間

這是我的模型。正如你可以看到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! 

這是將數據幀

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 
+0

可否請你分享的'dput(TURN)的結果',使示例可重現? –

+0

當然,我的不好。我使用可以使用的數據框編輯了問題。乾杯 – kumbu

回答

3

首先,打開mod_lmm1_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)) 

enter image description here