[R-lang] Plotting CI from Mixed Effects Model

[Roger Levy]

Hi David,

David Reitter wrote:
> Roger, Corey,
> 
> On 30 Sep 2007, at 00:47, Roger Levy wrote:
>> So with that explanation, the recommendation I'd give you for what you
>> want to do is: rerun lmer four different ways, each time resetting the
>> contrasts for place and voice such that you cover the four logical
>> possibilities, and report the CI for each cell on the basis of the CIs
>> (either classical or HPD) obtained when that cell is the default.
> 
> Do you see a problem with explicitly getting the model to estimate  
> the constrasts for all the factor combinations without parameters for  
> the main effects an the intercept, using
> 
> model <- glmmPQL(data=d, fixed=c ~ (a*b) - a - b - 1, random=...,  
> family=...)
> (This is with the MASS / nlme libraries, syntax from the top of my  
> head.)
> 
> and then using `intervals.lme',
> 
> intervals(model, level=0.95)
> 
> Would this produce the same results?

This is a good question.  I believe this would in general give the same 
results as long as you are absolutely sure to subtract out all the 
lower-order terms.   For the example I gave, it's fine:


 > lm.a.1 <- lm(y ~ x.a - 1)
 > summary(lm.a.1)

Call:
lm(formula = y ~ x.a - 1)

Coefficients:
      Estimate Std. Error t value Pr(>|t|)
x.aa    9.991      1.040   9.603 2.78e-14 ***
x.ab   12.789      0.658  19.436  < 2e-16 ***


Here's another example where it works:


 > n <- 50
 > s <- 5
 > a <- factor(ifelse(runif(n) > 0.5,"a1","a2"))
 > b <- factor(ifelse(runif(n) > 0.5,"b1","b2"))
 > ab <- ifelse(as.numeric(a) + as.numeric(b) == 4, 1,0)
 > model.mat <- cbind(rep(1,n),a,b,ab)
 > bet <- c(10,3,5,2)
 > pred <- model.mat %*% bet
 > y <- pred + rnorm(n,0,s)
 > plot(pred,y)
 >
 > lm.ab <- lm(y ~ a*b - a - b - 1)
 > summary(lm.ab)

Call:
lm(formula = y ~ a * b - a - b - 1)

Coefficients:
         Estimate Std. Error t value Pr(>|t|)
aa1:bb1   18.242      1.363   13.39   <2e-16 ***
aa2:bb1   19.797      1.493   13.26   <2e-16 ***
aa1:bb2   20.976      1.493   14.05   <2e-16 ***
aa2:bb2   27.069      1.113   24.33   <2e-16 ***


These coefficients shouldn't be correlated with each other. 
(unfortunately I don't see how to extract the covariance matrix for the 
model parameters from lm.  Do you know how?)

> 
>> This means that if you wanted to construct a simultaneous
>> confidence "interval" on multiple parameter of the model, it would
>> actually be a confidence *region* in the form of an ellipsoid. The
>> symmetric confidence intervals we look at around the point estimate
>> of a single parameter of the model have simply marginalized out all
>> the other model parameters.
> 
> In the above approach, we would basically only estimate parameters  
> for a single factor - which is simply the product of factors a and b,  
> e.g. all combinations of their respective levels. That way, there  
> shouldn't be interactions between multiple parameters.

Yes, that's right!

Roger

-- 

Roger Levy                      Email: rlevy at ucsd.edu
Assistant Professor             Phone: 858-534-7219
Department of Linguistics       Fax:   858-534-4789
UC San Diego                    Web:   http://ling.ucsd.edu/~rlevy