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

[Levy, Roger]

Hi Reinhold & co.,

I'm not sure whether this is exactly what Jon was suggesting, but it seems pretty clever!  Do you think you could point us to the handout and/or summarize this derivation?

Best

Roger


On Jul 30, 2011, at 4:29 PM, <kliegl@uni-potsdam.de>
 <kliegl@uni-potsdam.de> wrote:

> 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
>> 
> 
> 
> 
> 

--

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