Re: [R] Equivalent of intervals() in lmer

From: Michael Kubovy <kubovy_at_virginia.edu>
Date: Wed, 23 Apr 2008 11:14:23 -0400

Dear Friends,

(1) There may be a solution for those (e.g., experimental psychologists) who are *not at all* interested in generalizing the absolute level of the response variable (say, reaction time, rt) to other subjects, but *only* to generalize the effect of the manipulated variables within subjects (because that's what the theory speaks to) to other subjects. First, center the predictor(s) so as to remove any spurious correlation between mean and intercept, and then write the random effect w/o an intercept. Now the model will not address whether rt is different from 0, it will only estimate the slope and whether it's different from 0:

require(lme4)
require(gmodels)
data(sleepstudy)
ss <- sleepstudy

ss$days <- with(ss, Days - mean(Days))
(fm1 <- lmer(Reaction ~ days + (days|Subject), ss))
(fm3 <- lmer(Reaction ~ days + (-1 + days|Subject), ss))
ci(fm1)
ci(fm3) # CIs for days are about 14% smaller

(2) When the predictor is not continuous, this approach doesn't work. A solution for a designed experiment is to plot CIs for differences, i.e., 5%LSDs (rather than plot the CIs on cell means). Those CIs are smaller and address the question of interest. Here is a one-way ANOVA:

recall <- c(10, 13, 13, 6, 8, 8, 11, 14, 14, 22, 23, 25, 16, 18, 20, 15, 17, 17, 1, 1, 4, 12, 15, 17, 9, 12, 12, 8, 9, 12) fr <- data.frame(rcl = recall, time = factor(rep(c(1, 2, 5), 10)), subj = factor(rep(1:10, each = 3)))
(fr.lmer <- lmer(rcl ~ time -1 +(1 | subj), fr)) mm <- unique(model.matrix(~ time -1, fr)) cm <- mm[1, ] - mm[3, ]
cm1 <- mm[1, ] - mm[2, ]
estimable(fr.lmer, cm = cm, conf.into = 0.95) estimable(fr.lmer, cm = cm1, conf.into = 0.95)

plot mean \pm 2* 0.366 and call it a 5%LSD (uncorrected for multiple comparisons)

In the case of a more-than-one-way ANOVA with interaction (say 2x2), choose which simple effects are of interest, get SEs for those; and then plot the four points with the CIs. These are not quite right for the difference between simple effects, but I don't know what to do about that.



Professor Michael Kubovy
University of Virginia
Department of Psychology
USPS:     P.O.Box 400400    Charlottesville, VA 22904-4400
Parcels:    Room 102        Gilmer Hall
         McCormick Road    Charlottesville, VA 22903
Office:    B011    +1-434-982-4729
Lab:        B019    +1-434-982-4751
Fax:        +1-434-982-4766

WWW: http://www.people.virginia.edu/~mk9y/

On Apr 21, 2008, at 9:05 AM, Dieter Menne wrote:

> Douglas Bates <bates <at> stat.wisc.edu> writes:
>
>> If you want to examine the three means then you should fit the
>> model as
>> lmer(rcl ~ time - 1 + (1 | subj), fr)
>
> True, but for the notorious "error bars" in plots that reviewers
> always request
> the 0.35 is probable more relevant than the 1.87. Which I think is
> justified in
> this case, but in most non-orthogonal designs with three or more
> factors, where
> we have a mixture of between/withing subject, there is no clear
> solution. What
> to do when required to produce "error-bars" that reasonably mirror p-
> values?
>
> It's easier with British Journals in the medical field that often
> have
> statistical professionals as reviewers, but many American Journals
> with their
> amateur physician/statisticians (why no t-test on raw data?) drive
> me nuts.
>
> Dieter
>
> #-------------
> library(lme4)
> recall <- c(10, 13, 13, 6, 8, 8, 11, 14, 14, 22, 23, 25, 16, 18, 20,
> 15, 17, 17, 1, 1, 4, 12, 15, 17, 9, 12, 12, 8, 9, 12)
> fr <- data.frame(rcl = recall, time = factor(rep(c(1, 2, 5), 10)),
> subj = factor(rep(1:10, each = 3)))
> fr.lmer <- lmer(rcl ~ time -1 +(1 | subj), fr)
> summary(fr.lmer)
> fr.lmer <- lmer(rcl ~ time +(1 | subj), fr)
> summary(fr.lmer)
> ------------------
> Fixed effects:
> Estimate Std. Error t value
> time1 11.000 1.879 5.853
> time2 13.000 1.879 6.918
> time5 14.200 1.879 7.556
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 11.0000 1.8793 5.853
> time2 2.0000 0.3507 5.703
> time5 3.2000 0.3507 9.125

        [[alternative HTML version deleted]]



R-help_at_r-project.org 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. Received on Wed 23 Apr 2008 - 16:41:15 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Wed 23 Apr 2008 - 17:30:31 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive