R help - Oct 2006 - Variance Component/ICC Confidence Intervals via Bootstrap or Jackknife

I'm using the lme function in nmle to estimate the variance components
of a fully nested two-level model:

Y_ijk = mu + a_i + b_j(i) + e_k(j(i))

lme computes estimates of the variances for a, b, and e, call them v_a,
v_b, and v_e, and I can use the intervals function to get confidence
intervals. My understanding is that these intervals are probably not
that robust plus I need intervals on the intraclass correlation
coefficients:

v_a/(v_a + v_b + v_e)

and

(v_a + v_b)/(v_a + v_b + v_e).

I would like to use a bootstrap or a jackknife estimate of these
confidence intervals. Is there any package available?

I've tried to write the R code myself but I'm not sure if I understand
exactly how to do it. The way I've tried to do a bootstrap estimate is
to sample with replacement for i units, then sample with replacement the
j(i) units, and then finally sample with replacement the k(j(i)) units.

But the few references I've been able to track down (Arvesen, Biometrcs,
1970 is one), seem to say that I should just sample with replacement the
i units. Plus they seem to indicate that a log transform is needed. The
Arvesen reference used something like using log(v_a/v_e) as an estimator
for sigma_a^2/sigma_e^2 and getting an interval and then transforming to
get to an interval for the ICC (although it's not clear to me how to get
the other ICC in a two-level nested design).

Any insights would be appreciated.

Rick B.

I can think of two ways to get confidence intervals on intraclass 
correlations (ICCs) and more accurate intervals for variance 
components:  (1) modifying 'simulate.lme' to store the estimated 
variance components as well as "logLik" and (2) using 'lmer'
and
'mcmcsamp' in library(lme4). 

      The difficulty with (1) is that you have to make a local copy of 
'simulate.lme', then figure out where and how to modify it.  I've
just
looked at the code, and it looks like the required modifications should 
be fairly straightforward.  The problem with (2) is that you would have 
to learn the syntax for a nested model for 'lmer'.  It's different
from
that for 'lme' but not difficult.  The 'lmer' function is newer
and
better in many ways but is not as well documented and does not have as 
many helper functions yet.  The best documentation available for it may 
be the 'MlmSoftRev' vignette in the 'mlmRev' package plus the
"R News
5/1" article from May 2005.  If you are not familiar with vignettes, 
RSiteSearch("graves vignette") produced 74 hits for me just now.  Find
'vignette' in the first hit led me to an earlier description of
vignettes.

      If it were my problem, I'd probably try the second, though I might 
try both and compare.  If you try them both, I'd be interested in the 
comparison.  If the answers were substantially different, I'd worry. 

      Hope this helps. 
      Spencer Graves

Rick Bilonick wrote:
> I'm using the lme function in nmle to estimate the variance components
> of a fully nested two-level model:
>
> Y_ijk = mu + a_i + b_j(i) + e_k(j(i))
>
> lme computes estimates of the variances for a, b, and e, call them v_a,
> v_b, and v_e, and I can use the intervals function to get confidence
> intervals. My understanding is that these intervals are probably not
> that robust plus I need intervals on the intraclass correlation
> coefficients:
>
> v_a/(v_a + v_b + v_e)
>
> and
>
> (v_a + v_b)/(v_a + v_b + v_e).
>
> I would like to use a bootstrap or a jackknife estimate of these
> confidence intervals. Is there any package available?
>
> I've tried to write the R code myself but I'm not sure if I
understand
> exactly how to do it. The way I've tried to do a bootstrap estimate is
> to sample with replacement for i units, then sample with replacement the
> j(i) units, and then finally sample with replacement the k(j(i)) units.
>
> But the few references I've been able to track down (Arvesen,
Biometrcs,
> 1970 is one), seem to say that I should just sample with replacement the
> i units. Plus they seem to indicate that a log transform is needed. The
> Arvesen reference used something like using log(v_a/v_e) as an estimator
> for sigma_a^2/sigma_e^2 and getting an interval and then transforming to
> get to an interval for the ICC (although it's not clear to me how to
get
> the other ICC in a two-level nested design).
>
> Any insights would be appreciated.
>
> Rick B.
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>