2016-12-07 7 views
1

Ich führe einige Bootstrap-Konfidenzintervalle und möchte die Konfidenzintervalle mit dem Mittelwert darstellen. Etwas wie dieses: enter image description hereWie zeichne ich Bootstrapping-Konfidenzintervalle auf

Das ist mein Modell. Wie Sie sehen können, sind DISTANCE und POS Faktoren.

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) 

Dies ist die Funktion ich verwende:

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! 

Dies ist der TURN-Datenrahmen

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

Könnten Sie bitte die Ergebnisse der 'dput (TURN)' teilen, um die Beispiel reproduzierbar? –

+0

Sicher, mein Schlechter. Ich habe die Frage mit einem Datenrahmen bearbeitet, den Sie verwenden können. Prost – kumbu

Antwort

3

Zuerst das mod_lmm1_boot Objekt in einen Datenrahmen drehen Sie die verwenden können, tidy Methode für "Boot" -Objekte (von meinem broom Paket):

library(broom) 
tidy(mod_lmm1_boot) 

, die den Ausgang gibt:

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 

Sie können auch berechnen, die Konfidenzintervalle Sie mit so etwas wie Liste:

ci <- do.call(rbind, lapply(1:9, function(i) { 
    boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=i)$percent 
})) 

, die eine Matrix mit fünf Spalten gibt, wo die 4. und 5. sind Ihre prozentualen Konfidenzintervalle.

Sie etikettierte Ihre Frage mit ggplot2, der Code ist so hier die resultierenden Intervalle in ggplot2 plotten:

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

Verwandte Themen