[R-sig-ME] calculation of confidence intervals for random slope model

[Ben Bolker]

On 15-11-28 06:18 PM, W. Duncan Wadsworth wrote:
> Dr. Bolker,
> 
> Could you possibly expand on the statement "The other alternative is
> to use bootMer
> + predict to get confidence intervals ..."? In particular, what do you mean
> by using `predict` here?

  It took me a while to remember what the issue is here -- by default
`bootMer` *resamples* the random effects rather than conditioning on
them.  You need to use "use.u" or "re.form" to specify the
conditioning, as follows ...

sleepy = lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy)
## bootstrap summary function
sumfun = function(.) {
    coef(.)$Subject[,"Days"]
}
## where the fun happens
booty = as.data.frame(bootMer(sleepy, sumfun, nsim = 99,
       re.form=~Days|Subject))
## for plotting
rr <- ranef(sleepy)$Subject
colnames(booty) = rownames(rr)
slope_bs_distros = booty %>% tidyr::gather("Subject", "Value")
## bootstrap distributions of individual slope parameters (??)
## indiv_ests <- rr %>% add_rownames("Subject") %>%
##    mutate(Days=Days+fixef(sleepy)["Days"])
indiv_ests <- add_rownames(coef(sleepy)$Subject,"Subject")
ggplot(slope_bs_distros, aes(x = Value)) + geom_histogram() +
  facet_wrap(~ Subject, ncol = 5) +
  geom_vline(xintercept = 0, color = "red", size = 1.5)+
  geom_vline(aes(xintercept=Days),colour="blue",data=indiv_ests)


> 
> My first stab at this would look like:
> 
> library(lme4)
> library(ggplot2)
> library(dplyr)
> data("sleepstudy")
> ## visualize individual slopes
> ggplot(sleepstudy, aes(x = Days, y = Reaction, group = factor(Subject))) +
>   geom_smooth(aes(color = factor(Subject)), method = "lm", se = F) +
>   theme(legend.position = "none")
> sleepy = lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy)
> ## bootstrap summary function
> sumfun = function(.){
>   overall = fixef(.)
>   individuals = ranef(.)$Subject
>   return(c(overall["Days"] + individuals[,"Days"]))
> }
> ## where the fun happens
> booty = as.data.frame(bootMer(sleepy, sumfun, nsim = 99))
> ## for plotting
> colnames(booty) = rownames(ranef(sleepy)$Subject)
> slope_bs_distros = booty %>% tidyr::gather("Subject", "Value")
> ## bootstrap distributions of individual slope parameters (??)
> ggplot(slope_bs_distros, aes(x = Value)) + geom_histogram() +
>   facet_wrap(~ Subject, ncol = 5) +
>   geom_vline(xintercept = 0, color = "red", size = 1.5)
> 
> But this doesn't seem correct since Subject 335 appears to have a negative
> slope. Are we seeing some kind of shrinkage effect here or am I just doing
> the bootstrapping wrong?
> 
> 
> 
> On Mon, Nov 16, 2015 at 7:58 AM, Ben Bolker <bbolker at gmail.com> wrote:
> 
>> On Mon, Nov 16, 2015 at 5:56 AM, Henry Travers
>> <henry.travers at zoo.ox.ac.uk> wrote:
>>> I have what I hope is a relatively straightforward question about how to
>> interpret the results of a mixed effects model of the form:
>>>
>>> fm1 <- lmer(Reaction ~ Days + (Days | Subject))
>>>
>>> I am running an experiment such that I am most interested in the
>> (equivalent of the) effect of Days for each Subject, rather than say fitted
>> values. I understand how to derive the point estimates for this effect, but
>> I am struggling to see how to calculate confidence intervals for these
>> estimates that take account of both the standard error in the parameter
>> estimate for Days and the uncertainty in the corresponding slope estimates
>> for each Subject.
>>>
>>> I would be very grateful if someone could point me in the right
>> direction or to a suitable reference.
>>
>>   This is a surprisingly difficult question to answer.
>>   There have been extended discussions in the mailing list (which I
>> don't have time to dig for now) about whether it's OK to add the
>> conditional variances of the conditional modes to the variances of the
>> fixed-effect predictions, and the circumstances under which this would
>> be (in)accurate/(anti)conservative.  The other alternative is to use
>> bootMer + predict to get confidence intervals ...
>>
>>   This should probably be added to the FAQ ...
>>
>>>
>>> ------------------------------------------------
>>> Henry Travers, PhD
>>> Research Associate
>>>
>>> Interdisciplinary Centre for Conservation Science
>>> Department of Zoology
>>> University of Oxford
>>>
>>>
>>>
>>>
>>>
>>>
>>>         [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> R-sig-mixed-models at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>