[S] Summary Bootstrap Question

[Cox.Bill@epamail.epa.gov]

     Thanks toTim Hesterberg and Angelo Canty for responding to my post
(5/11/200) asking for help in performing block bootstrapping.   Both suggested
that I need to write my own sampling function that would accomodate sampling by
blocks within groups.   The two responses and suggested functions for doing
this:

Tim Hesterberg.

You wrote:
>       I am fitting a generalized model to daily data that extends
>over a number of years.  I want to be able to bootstrap the
>regression coefficients (and perhaps other stats) within groups
>(years) but also by blocks of say consecutive three day periods to
>preserve serial correlation known to exit.  The data looks something
>like this:
>
>   Date       Yr        Ozone         Temp           Wspd      Other variables
>   June 1     94           100             72                3.4
>........
>   June 2     94             85             88                4.3
>........
>   June 3     94             .                   .                    .
>.....
>     ........
>   Oct 30     99            120            76                4.5
>.........
>
> One form of the model being evaluated is:
>
>  Ozfit <- glm(Ozone ~ ns(temp,5) + ns(wspd,5 ) + .... + ns(yr,5),
>family=quasi(link=log),  data= ..  )
>I can bootstrap the reg coefficients nicely with
>
>  Ozfit.boot <-     bootstrap(data = ..., statistic = coef(eval(Ozfit$call)), B
>= 200, group = yr)
>
>Is there some way I can use the bootstrap function (or the boot
>function from statlib), to handle the serial correlation,
>e.g. bootsrap blocks of days so that each block of days would be
>sampled rather than each day?  Any other suggestions for handling
>this issue are welcome.  I am using Splus Pro 2000 release 2.
>Thanks.

You may pass a "sampler" of your choice to bootstrap(), including
a sampler that does block bootstrapping.
For example,

samp.boot.block <-
function(n, B, blockLength = 20)
{
  # Simple block bootstrap, overlapping blocks, no wrap-around,
  # no matching of ends.
  # To change blockLength use make.samp.boot.block()
  nblocks <- ceiling(n/blockLength)
  x <- sample(1:(n-blockLength+1), B*nblocks, replace=T)
  dim(x) <- c(nblocks, B)
  apply(x, 2, function(y, L) (0:(L-1)) + rep(y, each=L), L=blockLength)[1:n,]
}

make.samp.boot.block <- function(blockLength){
  # return a copy of samp.boot.block with the desired blockLength
  f <- samp.boot.block
  f$blockLength <- blockLength
  f
}
make.samp.boot.block(5)

samp.boot.block(5, 10, 3)
samp.boot.block(6, 10, 3)
bootstrap(1:25, mean)
bootstrap(1:25, mean, sampler=make.samp.boot.block(5)) # larger standard error

Angelo Canty

You essentially want to bootstrap a time series.  The function tsboot
(in my library from statlib) will do this.  One potential difficulty
is that you cannot specify strata (years) to this function but I am
not sure that you should be doing that in your case anyway.  If you
really want to include this then I'm afraid that you need to write
your own sampling function and use sim="parametric" in function boot.
One possible sampling function could be

blockboot.gen <- function(orig, mle) {
# mle is the block length.
    strata <- unique(orig$year)
    out <- orig
    for (s in strata) {
       dat <- orig[orig$year==s,]
       ns <- nrow(dat)
       nb <- ceiling(ns/mle)
       starts <- sample(ns-mle+1,nb,replace=T)
       inds <- apply(cbind(starts,mle-1),1,function(x)
x[1]:sum(x))[1:ns]
       out[out$year==s,] <- dat[inds,]
    }
    out
}

Note that since you are fitting a generalized linear model, you may
want to consider some form of model-based bootstrap.  You can still use
a block bootstrap as above but now just do this for some form of
residuals and add them back to the original model.  Please see chapters
7 and 8 of Davison and Hinkley (1997) for more information.

Thank to each of you.  Each suggested sampling function does the job nicely.


Bill Cox


-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu.  To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message:  unsubscribe s-news