R help - Apr 2010 - bootstrap confidence intervals, non iid

hello,

i need to calculate ci's for each of 4 groups within a dataset, to be able
to infere about differences in the variable "similarity". the problem
is
that data within groups is dependent, as assigned by the blocking-factor
"site". my guess was to use a block bootstrap but samples within in
these
blocks / sites are not of same length. i was not able to find a method to
deal with this -
 
i'd really appreciate any advise with this!

thanks,
kay

my data:
##########################################################
similarity<-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))))
##########################################################
-- 
View this message in context:
http://n4.nabble.com/bootstrap-confidence-intervals-non-iid-tp1751619p1751619.html
Sent from the R help mailing list archive at Nabble.com.

Kay Cichini wrote:
> 
> i need to calculate ci's for each of 4 groups within a dataset, to be
able
> to infere about differences in the variable "similarity". the
problem is
> that data within groups is dependent, as assigned by the blocking-factor
> "site". my guess was to use a block bootstrap but samples within
in these
> blocks / sites are not of same length. i was not able to find a method to
> deal with this -
> 
Have you considered bootstrapping the residuals?

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

not yet,
but it would either not solve the non iid problem, or would it?

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

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.