[R-lang] Re: p-values from pvals.fnc

Sat Jul 30 15:29:15 PDT 2011

Jon's proposal appears to work if we do not care about the intercept  
of X. Also, we need to estimate the parameters for the full model  
first to compute the X-offset for the zero intercept-slope correlation  
parameter. Specifically, if a and b are intercept and slope, then the  
offset for X to obtain a zero intercept-slope correlation parameter is  
cov(a,b)/ var(b). (I saw this derivation, which also yields a minimum  
variance estimate for the intercept, in a 2007 GLMM handout by Ulrich  
Halekoh.)  To illustrate this with the sleepstudy data.

> (fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
Linear mixed model fit by REML
Formula: Reaction ~ Days + (Days | Subject)
    Data: sleepstudy
   AIC  BIC logLik deviance REMLdev
  1756 1775 -871.8     1752    1744
Random effects:
  Groups   Name        Variance Std.Dev. Corr
  Subject  (Intercept) 612.092  24.7405
           Days         35.072   5.9221  0.066
  Residual             654.941  25.5918
Number of obs: 180, groups: Subject, 18

Fixed effects:
             Estimate Std. Error t value
(Intercept)  251.405      6.825   36.84
Days          10.467      1.546    6.77

# Compute the offset and add it to Days

> sleepstudy$DaysA <- sleepstudy$Days + (0.066*24.7405*5.9221)/35.072
> (fm3 <- lmer(Reaction ~ DaysA + (DaysA|Subject), sleepstudy))
Linear mixed model fit by REML
Formula: Reaction ~ DaysA + (DaysA | Subject)
    Data: sleepstudy
   AIC  BIC logLik deviance REMLdev
  1756 1775 -871.8     1752    1744
Random effects:
  Groups   Name        Variance Std.Dev. Corr
  Subject  (Intercept) 609.460  24.6872
           DaysA        35.072   5.9221  0.000
  Residual             654.941  25.5918
Number of obs: 180, groups: Subject, 18

Fixed effects:
             Estimate Std. Error t value
(Intercept)  248.519      6.896   36.04
DaysA         10.467      1.546    6.77

Now to the model with two independent variance components for which  
the MCMC p-values can be computed:
> (fm4 <- lmer(Reaction ~ DaysA + (1|Subject) + (0+DaysA|Subject), sleepstudy))
Linear mixed model fit by REML
Formula: Reaction ~ DaysA + (1 | Subject) + (0 + DaysA | Subject)
    Data: sleepstudy
   AIC  BIC logLik deviance REMLdev
  1754 1770 -871.8     1752    1744
Random effects:
  Groups   Name        Variance Std.Dev.
  Subject  (Intercept) 609.346  24.6849
  Subject  DaysA        35.066   5.9217
  Residual             654.951  25.5920
Number of obs: 180, groups: Subject, 18

Fixed effects:
             Estimate Std. Error t value
(Intercept)  248.519      6.896   36.04
DaysA         10.467      1.546    6.77

Not sure how well this generalizes but it appears that Jon's proposal  
could be implemented at least for the simple  
varying-intercept/varying-slope model.

Reinhold Kliegl
http://read.psych.uni-potsdam.de/pmr2

Quoting Jonathan Baron <baron@psych.upenn.edu>:

> I have a question about this, which has been bugging me for a while.
> It is very probably very naive, but here is a chance to ask it.
>
> pvals.fnc() does not deal with random slopes because it complains
> about "random correlation parameters".  I'm not quite sure what this
> means, but I think it may mean that it wants to take into account the
> correlation of the slopes and intercepts, but it can't.  Consider the
> following:
>
> l1 <- lmer(Y ~ X + (1+X|S)
> l2 <- lmer(Y ~ X + (1|S) + (0+X|S)
> l3 <- lmer(Y' ~ X' + (1|S) + (0+X'|S)
>
> pvals.fnc() does not work with l1.  It will work with l2, but I think
> l2 imposes the assumption of no correlation of slopes and intercepts.
> This may be approximately true, and maybe you can tell by looking at
> the output of l1, where it gives correlations in the random effects
> part.  (It also gives them for the fixed effects, but I think that is
> not the issue, although I'm not sure.)
>
> If the assumption of no correlation isn't true, can't you transform X
> and Y into X' and Y' (or maybe just transform one of them) by adding
> or subtracting a constant?  I have tried various ways of doing this,
> but none works for me except trial and error.  In many cases,
> transformation by a constant would not change the substantive
> inferences that we want to draw.
>
> Jon
>
> On 07/30/11 20:53, Jakke Tamminen wrote:
>> Roger: Thanks for the information, I guess I have a lot of reading to do!
>>
>> Alex and Roger: Looks like my version of lme4 is pretty old, 0.99875-6. If
>> the more recent versions don't give you p-values for models with random
>> slopes, should I be looking at them (the p-values) at all, or rely on the
>> t-statistic (and probably update my packages!)?
>>
>> Jakke
>>
>> On 30 July 2011 18:45, Alex Fine <afine@bcs.rochester.edu> wrote:
>>
>> > Jakke,
>> >
>> > I'm probably missing something, so I'm not replying-all.  How do you even
>> > get pvals.fnc() to work with a model that has random slopes?  I have the
>> > most up-to-date version of the languageR package and it won't take models
>> > that have anything other than random intercepts.
>> >
>> > thanks,
>> > Alex
>
> --
> Jonathan Baron, Professor of Psychology, University of Pennsylvania
> Home page: http://www.sas.upenn.edu/~baron
>