[R] bootstrap confidence intervals, non iid

[Kay Cichini]

hi glen,

i need conf.intervals for blocked data, as described in the first place.  
i've learned in the meantime, that the boot() function can handle this.

i had to formulate the function for the boot command, 
put "sites" to the strata argument and resample from each subsetted level of
the factor "stage".
with boot.ci() i then yielded the boot ci's for each stage.

hi glen,

i need conf.intervals for blocked data, as described in the first place.  
i've learned in the meantime, that the boot() function can handle this.

i needed to formulate the function for the boot command, 
put "sites" to the strata argument and to resample from each subsetted level
of the factor "stage".
with boot.ci() i then yielded the boot ci's for each stage.


here's the worked example, for anyone who is interested:

#######################################################################
#my data:
#######################################################################
sim<-data.frame(list(structure(list(stage = structure(c(1L, 1L, 1L, 1L, 1L,
1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L), .Label = c("A", "B", "C", "D"), class = "factor"), site =
structure(c(1L,
1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 5L, 6L, 6L,
6L, 6L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 10L, 10L,
10L, 10L, 11L, 11L, 12L, 12L, 12L, 13L, 13L, 13L, 14L, 14L, 14L,
14L, 15L, 15L, 15L, 15L, 16L, 16L, 16L, 17L, 17L, 17L, 17L, 18L,
18L, 19L, 19L, 19L, 19L, 20L, 20L, 20L, 20L, 21L, 21L, 21L, 21L,
22L, 22L, 22L, 22L, 23L, 23L, 23L, 24L, 24L, 24L, 24L, 25L, 25L,
25L, 25L, 26L, 26L, 26L, 26L, 27L, 27L, 27L, 27L, 28L, 28L, 28L,
28L, 29L, 29L, 29L, 30L, 30L, 30L, 30L, 31L, 31L, 32L, 32L, 32L,
32L, 33L, 33L, 33L, 33L, 34L, 34L, 34L, 34L, 35L, 35L, 35L, 35L,
36L, 36L, 36L, 36L, 37L, 37L, 38L, 38L, 38L, 38L, 39L, 39L, 39L
), .Label = c("A11", "A12", "A14", "A15", "A16", "A17", "A18",
"A19", "A20", "A5", "A7", "A8", "B1", "B12", "B13", "B14", "B15",
"B17", "B18", "B2", "B4", "B7", "B8", "B9", "C1", "C10", "C11",
"C15", "C17", "C18", "C19", "C2", "C20", "C3", "C4", "C6", "D1",
"D4", "D7"), class = "factor"), MH.Index = c(0.392156863, 0.602434077,
0.576923077, 0.647482014, 0.989010989, 0.857142857, 1, 1, 1,
0, 1, 0.378378378, 0.839087948, 0.252915554, 1, 0.22556391, 0.510366826,
0.476190476, 0.555819477, 0.961538462, 0.666666667, 0.089285714,
0.923076923, 0.571428571, 0, 0.923076923, 0.617647059, 0.599423631,
0, 0.727272727, 0.998112812, 0, 0, 0, 1, 0.565656566, 0.75, 0.923076923,
0.654545455, 0.14084507, 0.617647059, 0.315789474, 0.179347826,
0.583468021, 0.165525114, 0.817438692, 0.455551457, 0.49548886,
0.556127703, 0.707431246, 0.506757551, 0.689655172, 0.241433511,
0.379232506, 0.241935484, 0, 0.30848329, 0.530973451, 0.148148148,
0, 0.976744186, 0.550218341, 0.542168675, 0.769230769, 0.153310105,
0, 0, 0.380569406, 0.742174733, 0.222222222, 0.046925432, 0,
0.068076328, 0.772727273, 0.830039526, 0.503458415, 0.863910822,
0.39401263, 0.081818182, 0.368421053, 0.088607595, 0, 0.575499851,
0.605657238, 0.714854232, 0.855881172, 0.815689401, 0.552207228,
0.81708081, 0.583228133, 0.334466349, 0.259477365, 0.194711538,
0.278916707, 0.636304805, 0.593715432, 0.661016949, 0.626865672,
0.420219245, 0.453535143, 0.471243706, 0.462427746, 0.56980057,
0.453821155, 0.052828527, 0.926829268, 0.51988266, 0.472200264,
0.351219512, 0.290030211, 0.765258974, 0.564894108, 0.789699571,
0.863378215, 0.525181559, 0.803061458, 0.260164645, 0.477265792,
0.265889379, 0.317791411, 0.107623318, 0.279181709, 0.471953363,
0.463724265, 0.241966696, 0.403647213, 0.693087992, 0.494259925,
0.68904453, 0.39329147, 0.498161213, 0.376225983, 0.407001046,
0.825016633, 0.718991658, 0.662995912)), .Names = c("stage",
"site", "MH.Index"), class = "data.frame", row.names = c(NA,
-136L))))
#######################################################################
#my code:
#######################################################################
library(boot)
library(Hmisc)

####function:
mean.fun <- function(x, index){mean(x[index])}

X<-"A"
###the mean:
mean.fun(sim$MH.Index[sim$stage==X])

###the boot sample and the ci's: 
boot<-boot(sim$MH.Index[sim$stage==X], mean.fun, R=1000,
strata=sim$site[sim$stage==X])
ci<-boot.ci(boot,type = c("norm", "basic", "perc"))

###get ci's (method: normal)
ci[2]->meanA
data.frame(ci[4])[1,2]->lowA
data.frame(ci[4])[1,3]->uppA

X<-"B"
###the mean:
mean.fun(sim$MH.Index[sim$stage==X])

###the boot sample and the ci's:
boot<-boot(sim$MH.Index[sim$stage==X], mean.fun, R=1000,
strata=sim$site[sim$stage==X])
ci<-boot.ci(boot,type = c("norm", "basic", "perc"))

###get ci's (method: normal)
ci[2]->meanB
data.frame(ci[4])[1,2]->lowB
data.frame(ci[4])[1,3]->uppB

X<-"C"
###the mean:
mean.fun(sim$MH.Index[sim$stage==X])

###the boot sample and the ci's:
boot<-boot(sim$MH.Index[sim$stage==X], mean.fun, R=1000,
strata=sim$site[sim$stage==X])
ci<-boot.ci(boot,type = c("norm", "basic", "perc"))

###get ci's (method: normal)
data.frame(ci[2])[1]->meanC
data.frame(ci[4])[1,2]->lowC
data.frame(ci[4])[1,3]->uppC

X<-"D"
###the mean:
mean.fun(sim$MH.Index[sim$stage==X])

###the boot sample and the ci's:
boot<-boot(sim$MH.Index[sim$stage==X], mean.fun, R=1000,
strata=sim$site[sim$stage==X])
ci<-boot.ci(boot,type = c("norm", "basic", "perc"))

###get ci's (method: normal)
ci[2]->meanD
data.frame(ci[4])[1,2]->lowD
data.frame(ci[4])[1,3]->uppD


pldat<-data.frame(data.frame(mean=as.numeric(c(meanA,meanB,meanC,meanD))),
data.frame(low=as.numeric(c(lowA,lowB,lowC,lowD))),
data.frame(upp=as.numeric(c(uppA,uppB,uppC,uppD))))
rownames(pldat)<-c("A","B","C","D")
pldat

###plotting:
windows(4,6)
errbar(c(1:4),pldat$mean,pldat$upp,pldat$low,ylab="MH-Similarity",pch=15,xaxt="n",xlab="Stage")
axis(1,c(1:4),labels=row.names(pldat))
legend("top","errorbars =\n95% normal-\nbootstrap\nconf. int",bty="n")
#######################################################################


-- 
View this message in context: http://n4.nabble.com/bootstrap-confidence-intervals-non-iid-tp1751619p1819182.html
Sent from the R help mailing list archive at Nabble.com.