lme X lmer results

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

lme X lmer results

John Maindonald
Surely there is a correct denominator degrees of freedom if the design
is balanced, as Ronaldo's design seems to be. Assuming that he has
specified the design correctly to lme() and that lme() is getting the df
right, the difference is between 2 df and 878 df.  If the t-statistic  
for the
second level of Xvar had been 3.0 rather than 1.1, the difference
would be between a t-statistic equal to 0.095 and 1e-6.  In a design
where there are 10 observations on each experimental unit, and all
comparisons are at the level of experimental units or above, df for
all comparisons will be inflated by a factor of at least 9.

Rather than giving df that for the comparison(s) of interest may be
highly inflated, I'd prefer to give no degrees of freedom at all, & to
encourage users to work out df for themselves if at all possible.
If they are not able to do this, then mcmcsamp() is a good alternative,
and may be the way to go in any case.  This has the further advantage
of allowing assessments in cases where the relevant distribution is
hard to get at. I'd think a warning in order that the df are upper  
bounds,
and may be grossly inflated.

Incidentally, does mcmcsamp() do its calculations pretty well
independently of the lmer results?

John Maindonald.

On 29 Dec 2005, at 10:00 PM, [hidden email] wrote:

> From: Douglas Bates <[hidden email]>
> Date: 29 December 2005 5:59:07 AM
> To: "Ronaldo Reis-Jr." <[hidden email]>
> Cc: R-Help <[hidden email]>
> Subject: Re: [R] lme X lmer results
>
>
> On 12/26/05, Ronaldo Reis-Jr. <[hidden email]> wrote:
>> Hi,
>>
>> this is not a new doubt, but is a doubt that I cant find a good  
>> response.
>>
>> Look this output:
>>
>>> m.lme <- lme(Yvar~Xvar,random=~1|Plot1/Plot2/Plot3)
>>
>>> anova(m.lme)
>>             numDF denDF  F-value p-value
>> (Intercept)     1   860 210.2457  <.0001
>> Xvar            1     2   1.2352  0.3821
>>> summary(m.lme)
>> Linear mixed-effects model fit by REML
>>  Data: NULL
>>       AIC      BIC    logLik
>>   5416.59 5445.256 -2702.295
>>
>> Random effects:
>>  Formula: ~1 | Plot1
>>         (Intercept)
>> StdDev: 0.000745924
>>
>>  Formula: ~1 | Plot2 %in% Plot1
>>         (Intercept)
>> StdDev: 0.000158718
>>
>>  Formula: ~1 | Plot3 %in% Plot2 %in% Plot1
>>         (Intercept) Residual
>> StdDev: 0.000196583 5.216954
>>
>> Fixed effects: Yvar ~ Xvar
>>                    Value Std.Error  DF  t-value p-value
>> (Intercept)    2.3545454 0.2487091 860 9.467066  0.0000
>> XvarFactor2    0.3909091 0.3517278   2 1.111397  0.3821
>>
>> Number of Observations: 880
>> Number of Groups:
>>                          Plot1               Plot2 %in% Plot1
>>                              4                              8
>>    Plot3 %in% Plot2 %in% Plot1
>>                             20
>>
>> This is the correct result, de correct denDF for Xvar.
>>
>> I make this using lmer.
>>
>>> m.lmer <- lmer(Yvar~Xvar+(1|Plot1)+(1|Plot1:Plot2)+(1|Plot3))
>>> anova(m.lmer)
>> Analysis of Variance Table
>>            Df Sum Sq Mean Sq  Denom F value Pr(>F)
>> Xvar  1  33.62   33.62 878.00  1.2352 0.2667
>>> summary(m.lmer)
>> Linear mixed-effects model fit by REML
>> Formula: Yvar ~ Xvar + (1 | Plot1) + (1 | Plot1:Plot2) + (1 | Plot3)
>>      AIC     BIC    logLik MLdeviance REMLdeviance
>>  5416.59 5445.27 -2702.295   5402.698      5404.59
>> Random effects:
>>  Groups        Name        Variance   Std.Dev.
>>  Plot3         (Intercept) 1.3608e-08 0.00011665
>>  Plot1:Plot2   (Intercept) 1.3608e-08 0.00011665
>>  Plot1         (Intercept) 1.3608e-08 0.00011665
>>  Residual                  2.7217e+01 5.21695390
>> # of obs: 880, groups: Plot3, 20; Plot1:Plot2, 8; Plot1, 4
>>
>> Fixed effects:
>>                 Estimate Std. Error  DF t value Pr(>|t|)
>> (Intercept)      2.35455    0.24871 878  9.4671   <2e-16 ***
>> XvarFactor2      0.39091    0.35173 878  1.1114   0.2667
>> ---
>> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>>
>> Look the wrong P value, I know that it is wrong because the DF  
>> used. But, In
>> this case, the result is not correct. Dont have any difference of  
>> the result
>> using random effects with lmer and using a simple analyses with lm.
>
> You are assuming that there is a correct value of the denominator
> degrees of freedom.  I don't believe there is.  The statistic that is
> quoted there doesn't have exactly an F distribution so there is no
> correct degrees of freedom.
>
> One thing you can do with lmer is to form a Markov Chain Monte Carlo
> sample from the posterior distribution of the parameters so you can
> check to see whether the value of zero is in the middle of the
> distribution of XvarFactor2 or not.
>
> It would be possible for me to recreate in lmer the rules used in lme
> for calculating denominator degrees of freedom associated with terms
> of the random effects.  However, the class of models fit by lmer is
> larger than the class of models fit by lme (at least as far as the
> structure of the random-effects terms goes).  In particular lmer
> allows for random effects associated with crossed or partially crossed
> grouping factors and the rules for denominator degrees of freedom in
> lme only apply cleanly to nested grouping factors.  I would prefer to
> have a set of rules that would apply to the general case.
>
> Right now I would prefer to devote my time to other aspects of lmer -
> in particular I am still working on code for generalized linear mixed
> models using a supernodal Cholesky factorization.  I am willing to put
> this aside and code up the rules for denominator degrees of freedom
> with nested grouping factors BUT first I want someone to show me an
> example demonstrating that there really is a problem.  The example
> must show that the p-value calculated in the anova table or the
> parameter estimates table for lmer is seriously wrong compared to an
> empirical p-value - obtained from simulation under the null
> distribution or through MCMC sampling or something like that.  Saying
> that "Software XYZ says there are n denominator d.f. and lmer says
> there are m" does NOT count as an example.  I will readily concede
> that the denominator degrees of freedom reported by lmer are wrong but
> so are the degrees of freedom reported by Software XYZ because there
> is no right answer (in general - in a few simple balanced designs
> there may be a right answer).
>
>>
>>> m.lm <- lm(Yvar~Xvar)
>>>
>>> anova(m.lm)
>> Analysis of Variance Table
>>
>> Response: Nadultos
>>             Df  Sum Sq Mean Sq F value Pr(>F)
>> Xvar         1    33.6    33.6  1.2352 0.2667
>> Residuals  878 23896.2    27.2
>>>
>>> summary(m.lm)
>>
>> Call:
>> lm(formula = Yvar ~ Xvar)
>>
>> Residuals:
>>     Min      1Q  Median      3Q     Max
>> -2.7455 -2.3545 -1.7455  0.2545 69.6455
>>
>> Coefficients:
>>                Estimate Std. Error t value Pr(>|t|)
>> (Intercept)      2.3545     0.2487   9.467   <2e-16 ***
>> XvarFactor2      0.3909     0.3517   1.111    0.267
>> ---
>> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>>
>> Residual standard error: 5.217 on 878 degrees of freedom
>> Multiple R-Squared: 0.001405,   Adjusted R-squared: 0.0002675
>> F-statistic: 1.235 on 1 and 878 DF,  p-value: 0.2667
>>
>> I read the rnews about this use of the full DF in lmer, but I dont  
>> undestand
>> this use with a gaussian error, I undestand this with glm data.
>>
>> I need more explanations, please.
>>
>> Thanks
>> Ronaldo

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: lme X lmer results

Douglas Bates
On 12/29/05, John Maindonald <[hidden email]> wrote:

> Surely there is a correct denominator degrees of freedom if the design
> is balanced, as Ronaldo's design seems to be. Assuming that he has
> specified the design correctly to lme() and that lme() is getting the df
> right, the difference is between 2 df and 878 df.  If the t-statistic
> for the
> second level of Xvar had been 3.0 rather than 1.1, the difference
> would be between a t-statistic equal to 0.095 and 1e-6.  In a design
> where there are 10 observations on each experimental unit, and all
> comparisons are at the level of experimental units or above, df for
> all comparisons will be inflated by a factor of at least 9.

I don't want to be obtuse and argumentative but I still am not
convinced that there is a correct denominator degrees of freedom for
_this_ F statistic.  I may be wrong about this but I think you are
referring to an F statistic based on a denominator from a different
error stratum, which is not what is being quoted.  (Those are not
given because they don't generalize to unbalanced designs.)

This is why I would like to see someone undertake a simulation study
to compare various approaches to inference for the fixed effects terms
in a mixed model, using realistic (i.e. unbalanced) examples.

It seems peculiar to me that the F statistics are being created from
the ratios of mean squares for different terms to the _same_ mean
square (actually a penalized sum of squares divided by the degrees of
freedom) and the adjustment suggested to take into account the
presence of the random effects is to change the denominator degrees of
freedom.  I think the rationale for this is an attempt to generalized
another approach (the use of error strata) even though it is not being
used here.

> Rather than giving df that for the comparison(s) of interest may be
> highly inflated, I'd prefer to give no degrees of freedom at all, & to
> encourage users to work out df for themselves if at all possible.
> If they are not able to do this, then mcmcsamp() is a good alternative,
> and may be the way to go in any case.  This has the further advantage
> of allowing assessments in cases where the relevant distribution is
> hard to get at. I'd think a warning in order that the df are upper
> bounds,
> and may be grossly inflated.

As I said, I am willing to change this if it is shown to be grossly
inaccurate but please show me.

> Incidentally, does mcmcsamp() do its calculations pretty well
> independently of the lmer results?

mcmcsamp starts from the parameter estimates when creating the chain
but that is the extent to which it depends on the lmer results.

> John Maindonald.
>
> On 29 Dec 2005, at 10:00 PM, [hidden email] wrote:
>
> > From: Douglas Bates <[hidden email]>
> > Date: 29 December 2005 5:59:07 AM
> > To: "Ronaldo Reis-Jr." <[hidden email]>
> > Cc: R-Help <[hidden email]>
> > Subject: Re: [R] lme X lmer results
> >
> >
> > On 12/26/05, Ronaldo Reis-Jr. <[hidden email]> wrote:
> >> Hi,
> >>
> >> this is not a new doubt, but is a doubt that I cant find a good
> >> response.
> >>
> >> Look this output:
> >>
> >>> m.lme <- lme(Yvar~Xvar,random=~1|Plot1/Plot2/Plot3)
> >>
> >>> anova(m.lme)
> >>             numDF denDF  F-value p-value
> >> (Intercept)     1   860 210.2457  <.0001
> >> Xvar            1     2   1.2352  0.3821
> >>> summary(m.lme)
> >> Linear mixed-effects model fit by REML
> >>  Data: NULL
> >>       AIC      BIC    logLik
> >>   5416.59 5445.256 -2702.295
> >>
> >> Random effects:
> >>  Formula: ~1 | Plot1
> >>         (Intercept)
> >> StdDev: 0.000745924
> >>
> >>  Formula: ~1 | Plot2 %in% Plot1
> >>         (Intercept)
> >> StdDev: 0.000158718
> >>
> >>  Formula: ~1 | Plot3 %in% Plot2 %in% Plot1
> >>         (Intercept) Residual
> >> StdDev: 0.000196583 5.216954
> >>
> >> Fixed effects: Yvar ~ Xvar
> >>                    Value Std.Error  DF  t-value p-value
> >> (Intercept)    2.3545454 0.2487091 860 9.467066  0.0000
> >> XvarFactor2    0.3909091 0.3517278   2 1.111397  0.3821
> >>
> >> Number of Observations: 880
> >> Number of Groups:
> >>                          Plot1               Plot2 %in% Plot1
> >>                              4                              8
> >>    Plot3 %in% Plot2 %in% Plot1
> >>                             20
> >>
> >> This is the correct result, de correct denDF for Xvar.
> >>
> >> I make this using lmer.
> >>
> >>> m.lmer <- lmer(Yvar~Xvar+(1|Plot1)+(1|Plot1:Plot2)+(1|Plot3))
> >>> anova(m.lmer)
> >> Analysis of Variance Table
> >>            Df Sum Sq Mean Sq  Denom F value Pr(>F)
> >> Xvar  1  33.62   33.62 878.00  1.2352 0.2667
> >>> summary(m.lmer)
> >> Linear mixed-effects model fit by REML
> >> Formula: Yvar ~ Xvar + (1 | Plot1) + (1 | Plot1:Plot2) + (1 | Plot3)
> >>      AIC     BIC    logLik MLdeviance REMLdeviance
> >>  5416.59 5445.27 -2702.295   5402.698      5404.59
> >> Random effects:
> >>  Groups        Name        Variance   Std.Dev.
> >>  Plot3         (Intercept) 1.3608e-08 0.00011665
> >>  Plot1:Plot2   (Intercept) 1.3608e-08 0.00011665
> >>  Plot1         (Intercept) 1.3608e-08 0.00011665
> >>  Residual                  2.7217e+01 5.21695390
> >> # of obs: 880, groups: Plot3, 20; Plot1:Plot2, 8; Plot1, 4
> >>
> >> Fixed effects:
> >>                 Estimate Std. Error  DF t value Pr(>|t|)
> >> (Intercept)      2.35455    0.24871 878  9.4671   <2e-16 ***
> >> XvarFactor2      0.39091    0.35173 878  1.1114   0.2667
> >> ---
> >> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> >>
> >> Look the wrong P value, I know that it is wrong because the DF
> >> used. But, In
> >> this case, the result is not correct. Dont have any difference of
> >> the result
> >> using random effects with lmer and using a simple analyses with lm.
> >
> > You are assuming that there is a correct value of the denominator
> > degrees of freedom.  I don't believe there is.  The statistic that is
> > quoted there doesn't have exactly an F distribution so there is no
> > correct degrees of freedom.
> >
> > One thing you can do with lmer is to form a Markov Chain Monte Carlo
> > sample from the posterior distribution of the parameters so you can
> > check to see whether the value of zero is in the middle of the
> > distribution of XvarFactor2 or not.
> >
> > It would be possible for me to recreate in lmer the rules used in lme
> > for calculating denominator degrees of freedom associated with terms
> > of the random effects.  However, the class of models fit by lmer is
> > larger than the class of models fit by lme (at least as far as the
> > structure of the random-effects terms goes).  In particular lmer
> > allows for random effects associated with crossed or partially crossed
> > grouping factors and the rules for denominator degrees of freedom in
> > lme only apply cleanly to nested grouping factors.  I would prefer to
> > have a set of rules that would apply to the general case.
> >
> > Right now I would prefer to devote my time to other aspects of lmer -
> > in particular I am still working on code for generalized linear mixed
> > models using a supernodal Cholesky factorization.  I am willing to put
> > this aside and code up the rules for denominator degrees of freedom
> > with nested grouping factors BUT first I want someone to show me an
> > example demonstrating that there really is a problem.  The example
> > must show that the p-value calculated in the anova table or the
> > parameter estimates table for lmer is seriously wrong compared to an
> > empirical p-value - obtained from simulation under the null
> > distribution or through MCMC sampling or something like that.  Saying
> > that "Software XYZ says there are n denominator d.f. and lmer says
> > there are m" does NOT count as an example.  I will readily concede
> > that the denominator degrees of freedom reported by lmer are wrong but
> > so are the degrees of freedom reported by Software XYZ because there
> > is no right answer (in general - in a few simple balanced designs
> > there may be a right answer).
> >
> >>
> >>> m.lm <- lm(Yvar~Xvar)
> >>>
> >>> anova(m.lm)
> >> Analysis of Variance Table
> >>
> >> Response: Nadultos
> >>             Df  Sum Sq Mean Sq F value Pr(>F)
> >> Xvar         1    33.6    33.6  1.2352 0.2667
> >> Residuals  878 23896.2    27.2
> >>>
> >>> summary(m.lm)
> >>
> >> Call:
> >> lm(formula = Yvar ~ Xvar)
> >>
> >> Residuals:
> >>     Min      1Q  Median      3Q     Max
> >> -2.7455 -2.3545 -1.7455  0.2545 69.6455
> >>
> >> Coefficients:
> >>                Estimate Std. Error t value Pr(>|t|)
> >> (Intercept)      2.3545     0.2487   9.467   <2e-16 ***
> >> XvarFactor2      0.3909     0.3517   1.111    0.267
> >> ---
> >> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> >>
> >> Residual standard error: 5.217 on 878 degrees of freedom
> >> Multiple R-Squared: 0.001405,   Adjusted R-squared: 0.0002675
> >> F-statistic: 1.235 on 1 and 878 DF,  p-value: 0.2667
> >>
> >> I read the rnews about this use of the full DF in lmer, but I dont
> >> undestand
> >> this use with a gaussian error, I undestand this with glm data.
> >>
> >> I need more explanations, please.
> >>
> >> Thanks
> >> Ronaldo
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: lme X lmer results

David Atkins
In reply to this post by John Maindonald

Message: 18
Date: Fri, 30 Dec 2005 12:51:59 -0600
From: Douglas Bates <[hidden email]>
Subject: Re: [R] lme X lmer results
To: John Maindonald <[hidden email]>
Cc: [hidden email]
Message-ID:
        <[hidden email]>
Content-Type: text/plain; charset=ISO-8859-1

On 12/29/05, John Maindonald <[hidden email]> wrote:

 >> Surely there is a correct denominator degrees of freedom if the design
 >> is balanced, as Ronaldo's design seems to be. Assuming that he has
 >> specified the design correctly to lme() and that lme() is getting the df
 >> right, the difference is between 2 df and 878 df.  If the t-statistic
 >> for the
 >> second level of Xvar had been 3.0 rather than 1.1, the difference
 >> would be between a t-statistic equal to 0.095 and 1e-6.  In a design
 >> where there are 10 observations on each experimental unit, and all
 >> comparisons are at the level of experimental units or above, df for
 >> all comparisons will be inflated by a factor of at least 9.

Doug Bates commented:

I don't want to be obtuse and argumentative but I still am not
convinced that there is a correct denominator degrees of freedom for
_this_ F statistic.  I may be wrong about this but I think you are
referring to an F statistic based on a denominator from a different
error stratum, which is not what is being quoted.  (Those are not
given because they don't generalize to unbalanced designs.)

This is why I would like to see someone undertake a simulation study
to compare various approaches to inference for the fixed effects terms
in a mixed model, using realistic (i.e. unbalanced) examples.

Doug--

Here is a paper that focused on the various alternatives to denominator degrees
of freedom in SAS and does report some simulation results:

http://www2.sas.com/proceedings/sugi26/p262-26.pdf

Not sure whether it argues convincingly one way or the other in the present
discussion.

cheers, Dave

--
Dave Atkins, PhD
[hidden email]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: lme X lmer results

Peter Dalgaard
Dave Atkins <[hidden email]> writes:

> Message: 18
> Date: Fri, 30 Dec 2005 12:51:59 -0600
> From: Douglas Bates <[hidden email]>
> Subject: Re: [R] lme X lmer results
> To: John Maindonald <[hidden email]>
> Cc: [hidden email]
> Message-ID:
> <[hidden email]>
> Content-Type: text/plain; charset=ISO-8859-1
>
> On 12/29/05, John Maindonald <[hidden email]> wrote:
>
>  >> Surely there is a correct denominator degrees of freedom if the design
>  >> is balanced, as Ronaldo's design seems to be. Assuming that he has
>  >> specified the design correctly to lme() and that lme() is getting the df
>  >> right, the difference is between 2 df and 878 df.  If the t-statistic
>  >> for the
>  >> second level of Xvar had been 3.0 rather than 1.1, the difference
>  >> would be between a t-statistic equal to 0.095 and 1e-6.  In a design
>  >> where there are 10 observations on each experimental unit, and all
>  >> comparisons are at the level of experimental units or above, df for
>  >> all comparisons will be inflated by a factor of at least 9.
>
> Doug Bates commented:
>
> I don't want to be obtuse and argumentative but I still am not
> convinced that there is a correct denominator degrees of freedom for
> _this_ F statistic.  I may be wrong about this but I think you are
> referring to an F statistic based on a denominator from a different
> error stratum, which is not what is being quoted.  (Those are not
> given because they don't generalize to unbalanced designs.)
>
> This is why I would like to see someone undertake a simulation study
> to compare various approaches to inference for the fixed effects terms
> in a mixed model, using realistic (i.e. unbalanced) examples.
>
> Doug--
>
> Here is a paper that focused on the various alternatives to denominator degrees
> of freedom in SAS and does report some simulation results:
>
> http://www2.sas.com/proceedings/sugi26/p262-26.pdf
>
> Not sure whether it argues convincingly one way or the other in the present
> discussion.
>
> cheers, Dave

It's certainly informative. I was not quite aware that what SAS calls
"Satterthwaite" could perform so poorly relative to "KenwardRoger"
(which to my mind is much closer in structure to a generalized
Satterthwaite approximation).

In either case I think it might be worth emphasizing that the issue is
not really whether you are testing at alpha=.050 or at alpha=.093 --
the corrections are likely to be so highly dependent on higher order
moments of the normal distribution to be wildly off in practice. The
important thing is that you get a low denominator DF value when you
are far from "Asymptopia" and any good statistician should know better
than to trust such tests.

The nasty problem with the "containment" type methods is that they can
get things completely wrong and, e.g., give DFs on the order of the
number of observations where in truth (insofar as it exists) it should
be closer to the number of subjects in the study.

I don't think anyone, including Doug, would be opposed to including a
Kenward-Roger style DF calculation in lmer. It "just" has to be worked
out how to convert the calculations to work with the sparse-matrix,
penalized least squares techniques that it uses, and Doug himself has
his mind elsewhere.

--
   O__  ---- Peter Dalgaard             Ă˜ster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics     PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark          Ph:  (+45) 35327918
~~~~~~~~~~ - ([hidden email])                  FAX: (+45) 35327907

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: lme X lmer results

John Maindonald
In reply to this post by Douglas Bates
Douglas -

As I understand Ronaldo's experiment, there are 4 plots, 8 subplots  
within
each of those 4 plots, and 20 subsubplots within each of the 8 subplots.
Within each subsbubplot there is an average of 1.375 observational
units.  We do not however need to know the distribution of observational
units per plot, for the discussion that follows.

A treatment Xvar, with 2 levels, has been applied at the plot level.
So 2 plots get the level 1 treatment, and 2 plots get the level 2  
treatment.
A simple-minded (and insightful) way to start thinking about the  
analysis
is to calculate averages for each of the 4 plots, and base the  
analysis on
those averages.

There are 220 observational units (SD 5.217), 32 plot 1 units (SD  
0.001965),
8 plot 2 units (SD 0.001587), together with an SD of 0.000741 at the  
plot 1,
that contribute to variation at the plot 1 level.  So the SD for an  
inidividual
plot 1 unit is
sd1 = sqrt(5.217^2/220+.001965^2/32+.001587^2/8+.000741^2) = 0.35173

Thus the estimated SED for comparing the two levels of Xvar is
sqrt(s^2/2+s^2/2) = s = 0.35173
Not accidentally, this is exactly the SED that is given by both lme  
and lmer.
Because the estimate is dominated by variation at the level of  
observational
units this is also, essentially, the value given by lm.

If we knew that the variance component at the plot 1 level was indeed  
very
small, giving the variance estimate 2df would be very conservative.  The
problem is that we do not know this.  We can either calculate sd1^2  
directly
from the plot 1 means, or we can go through the more convoluted process
of estimating the various variance components, then dividing by the  
relevant
numbers and adding.  Either way, both the estimate of sd1 and the  
estimate
of the component of variance at the plot 1 level is a 2df  
calculation, with a
CI that relies on assuming that sd1^2 has something like a chi-squared
distribution.  For moderately unbalanced designs the df calculation  
has to
be fudged a bit.

The following function has defaults that reproduce the essential  
features
or Ronaldo's analysis, though with a treatment effect that is noise.  
For
simplicity, there are just two levels of variation:

"mlsim" <-
   function(s1=0.00074, s2=5.22, n1=4, n2=220, nsamples=1000){
     ## Generate data with n1 plots, each with 220 subplots
     ## Data are such that estimated variance component at plot
     ## level is s1^2; at subplot level is s2^2
     ## A treatment with 2 levels is applied at the plot level;
     ## with n1/2 plots per treatment.
     ## n1 should be an even number
     nu1 <- n1-2       ## df for treatment comparison
     nu2 <- n1*(n2-1)  ## df for estimate of s2
     xp <- expand.grid(plot1=1:n1, plot2=1:n2)
     eff1 <- rnorm(n1)
     trt <- factor((xp$plot1+1) %% 2 + 1)
     sd1 <- sqrt(s1^2+s2^2/n2)
     eff1 <- residuals(lm(eff1~factor((2:(n1+1)) %% 2 + 1)))
     eff1 <- eff1*sd1/sd(eff1)
     eff2 <- rnorm(n1*n2)
     eff2 <- resid(lm(eff2~factor(xp$plot1)))
     eff2 <- eff2*s2/sqrt(sum(eff2^2)/nu2)
     y <- eff1[xp$plot1] + eff2
     df <- cbind(xp, trt=trt, y=y)
     df$plot1 <- factor(df$plot1)
     df$plot2 <- factor(df$plot2)
     asum <- summary(aov(y~trt +Error(plot1), data=df))
     bms2 <- asum[[1]][[1]][2,3]
     ws2 <- asum[[2]][[1]][1,3]
     eff1 <- eff1*sd1/sqrt(bms2/n2)
     df$y <- eff1[xp$plot1] + eff2
     if(nsamples>1){
       df.lmer <- lmer(y~trt + (1|plot1), data=df)
       print(summary(df.lmer))
       df.mc <- mcmcsamp(df.lmer, n=1000)
       par(mfrow=c(2,2), mar=c(3.1,3.1,2.1,1.1), mgp=c(2,.5,0))
       on.exit(par(oldpar))
       qqplot(qchisq(ppoints(1000),nu2)*s2^2/nu2, exp(df.mc[,3]))
       mtext(side=3,line=0.25, paste("s2 =", s2, " (nu2=", nu2, ")",  
sep=""))
       abline(0,1,col=2)
       qqplot(qchisq(ppoints(1000), nu1)*sd1^2/nu1,
              exp(df.mc[,3])/n2+exp(df.mc[,4]))
       abline(0,1,col=2)
       mtext(side=3,line=1.25, paste("s1=",s1, ";  Var of gp means=",
                      paste(round(sd1^2,2)), ";  df=", nu1, sep=""))
       mtext(side=3,line=.25, paste("(",signif(s1^2, 4), " (from s1); ",
                      signif(s2^2/n2,4)," ((from s2)", ")", sep=""))
       invisible(df.mc)
     }
     else invisible(df)
   }

Observe that the sample distribution of sd1^2 (the estimated variance
of the means at the plot 1 level) changes quite dramatically from one
run to the next, withe the direction of convexity in the qqplot changing
around the line y=x.  The distribution is not chi-squared 2, but it  
is not
consistently anything else either.  The distribution is very likely,  
I am
guessing, sensitive to the choice of prior.  Possibly a choice where
log(sigma2) is locally uniform, where sigma2^2 is the variance at
the plot 2 level, would yield a distribution that is closer to chi-
squared 2.
The plots do however indicate that the chi-squared 2 distribution is in
the right ballpark, as judged by the results from multiple runs of
mcmcsamp().

With mlsim(s1=2.5, s2=5, n1=40, n2=4), the two components contribute
equally to sd1^2, and sd1^2 usually tracks quite closely to a
chi-squared distribution.

John Maindonald.


On 31 Dec 2005, at 5:51 AM, Douglas Bates wrote:

> On 12/29/05, John Maindonald <[hidden email]> wrote:
>> Surely there is a correct denominator degrees of freedom if the  
>> design
>> is balanced, as Ronaldo's design seems to be. Assuming that he has
>> specified the design correctly to lme() and that lme() is getting  
>> the df
>> right, the difference is between 2 df and 878 df.  If the t-statistic
>> for the
>> second level of Xvar had been 3.0 rather than 1.1, the difference
>> would be between a t-statistic equal to 0.095 and 1e-6.  In a design
>> where there are 10 observations on each experimental unit, and all
>> comparisons are at the level of experimental units or above, df for
>> all comparisons will be inflated by a factor of at least 9.
>
> I don't want to be obtuse and argumentative but I still am not
> convinced that there is a correct denominator degrees of freedom for
> _this_ F statistic.  I may be wrong about this but I think you are
> referring to an F statistic based on a denominator from a different
> error stratum, which is not what is being quoted.  (Those are not
> given because they don't generalize to unbalanced designs.)
>
> This is why I would like to see someone undertake a simulation study
> to compare various approaches to inference for the fixed effects terms
> in a mixed model, using realistic (i.e. unbalanced) examples.
>
> It seems peculiar to me that the F statistics are being created from
> the ratios of mean squares for different terms to the _same_ mean
> square (actually a penalized sum of squares divided by the degrees of
> freedom) and the adjustment suggested to take into account the
> presence of the random effects is to change the denominator degrees of
> freedom.  I think the rationale for this is an attempt to generalized
> another approach (the use of error strata) even though it is not being
> used here.
>
>> Rather than giving df that for the comparison(s) of interest may be
>> highly inflated, I'd prefer to give no degrees of freedom at all,  
>> & to
>> encourage users to work out df for themselves if at all possible.
>> If they are not able to do this, then mcmcsamp() is a good  
>> alternative,
>> and may be the way to go in any case.  This has the further advantage
>> of allowing assessments in cases where the relevant distribution is
>> hard to get at. I'd think a warning in order that the df are upper
>> bounds,
>> and may be grossly inflated.
>
> As I said, I am willing to change this if it is shown to be grossly
> inaccurate but please show me.
>
>> Incidentally, does mcmcsamp() do its calculations pretty well
>> independently of the lmer results?
>
> mcmcsamp starts from the parameter estimates when creating the chain
> but that is the extent to which it depends on the lmer results.
>
>> John Maindonald.
>>
>> On 29 Dec 2005, at 10:00 PM, [hidden email] wrote:
>>
>>> From: Douglas Bates <[hidden email]>
>>> Date: 29 December 2005 5:59:07 AM
>>> To: "Ronaldo Reis-Jr." <[hidden email]>
>>> Cc: R-Help <[hidden email]>
>>> Subject: Re: [R] lme X lmer results
>>>
>>>
>>> On 12/26/05, Ronaldo Reis-Jr. <[hidden email]> wrote:
>>>> Hi,
>>>>
>>>> this is not a new doubt, but is a doubt that I cant find a good
>>>> response.
>>>>
>>>> Look this output:
>>>>
>>>>> m.lme <- lme(Yvar~Xvar,random=~1|Plot1/Plot2/Plot3)
>>>>
>>>>> anova(m.lme)
>>>>             numDF denDF  F-value p-value
>>>> (Intercept)     1   860 210.2457  <.0001
>>>> Xvar            1     2   1.2352  0.3821
>>>>> summary(m.lme)
>>>> Linear mixed-effects model fit by REML
>>>>  Data: NULL
>>>>       AIC      BIC    logLik
>>>>   5416.59 5445.256 -2702.295
>>>>
>>>> Random effects:
>>>>  Formula: ~1 | Plot1
>>>>         (Intercept)
>>>> StdDev: 0.000745924
>>>>
>>>>  Formula: ~1 | Plot2 %in% Plot1
>>>>         (Intercept)
>>>> StdDev: 0.000158718
>>>>
>>>>  Formula: ~1 | Plot3 %in% Plot2 %in% Plot1
>>>>         (Intercept) Residual
>>>> StdDev: 0.000196583 5.216954
>>>>
>>>> Fixed effects: Yvar ~ Xvar
>>>>                    Value Std.Error  DF  t-value p-value
>>>> (Intercept)    2.3545454 0.2487091 860 9.467066  0.0000
>>>> XvarFactor2    0.3909091 0.3517278   2 1.111397  0.3821
>>>>
>>>> Number of Observations: 880
>>>> Number of Groups:
>>>>                          Plot1               Plot2 %in% Plot1
>>>>                              4                              8
>>>>    Plot3 %in% Plot2 %in% Plot1
>>>>                             20
>>>>
>>>> This is the correct result, de correct denDF for Xvar.
>>>>
>>>> I make this using lmer.
>>>>
>>>>> m.lmer <- lmer(Yvar~Xvar+(1|Plot1)+(1|Plot1:Plot2)+(1|Plot3))
>>>>> anova(m.lmer)
>>>> Analysis of Variance Table
>>>>            Df Sum Sq Mean Sq  Denom F value Pr(>F)
>>>> Xvar  1  33.62   33.62 878.00  1.2352 0.2667
>>>>> summary(m.lmer)
>>>> Linear mixed-effects model fit by REML
>>>> Formula: Yvar ~ Xvar + (1 | Plot1) + (1 | Plot1:Plot2) + (1 |  
>>>> Plot3)
>>>>      AIC     BIC    logLik MLdeviance REMLdeviance
>>>>  5416.59 5445.27 -2702.295   5402.698      5404.59
>>>> Random effects:
>>>>  Groups        Name        Variance   Std.Dev.
>>>>  Plot3         (Intercept) 1.3608e-08 0.00011665
>>>>  Plot1:Plot2   (Intercept) 1.3608e-08 0.00011665
>>>>  Plot1         (Intercept) 1.3608e-08 0.00011665
>>>>  Residual                  2.7217e+01 5.21695390
>>>> # of obs: 880, groups: Plot3, 20; Plot1:Plot2, 8; Plot1, 4
>>>>
>>>> Fixed effects:
>>>>                 Estimate Std. Error  DF t value Pr(>|t|)
>>>> (Intercept)      2.35455    0.24871 878  9.4671   <2e-16 ***
>>>> XvarFactor2      0.39091    0.35173 878  1.1114   0.2667
>>>> ---
>>>> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>>>>
>>>> Look the wrong P value, I know that it is wrong because the DF
>>>> used. But, In
>>>> this case, the result is not correct. Dont have any difference of
>>>> the result
>>>> using random effects with lmer and using a simple analyses with lm.
>>>
>>> You are assuming that there is a correct value of the denominator
>>> degrees of freedom.  I don't believe there is.  The statistic  
>>> that is
>>> quoted there doesn't have exactly an F distribution so there is no
>>> correct degrees of freedom.
>>>
>>> One thing you can do with lmer is to form a Markov Chain Monte Carlo
>>> sample from the posterior distribution of the parameters so you can
>>> check to see whether the value of zero is in the middle of the
>>> distribution of XvarFactor2 or not.
>>>
>>> It would be possible for me to recreate in lmer the rules used in  
>>> lme
>>> for calculating denominator degrees of freedom associated with terms
>>> of the random effects.  However, the class of models fit by lmer is
>>> larger than the class of models fit by lme (at least as far as the
>>> structure of the random-effects terms goes).  In particular lmer
>>> allows for random effects associated with crossed or partially  
>>> crossed
>>> grouping factors and the rules for denominator degrees of freedom in
>>> lme only apply cleanly to nested grouping factors.  I would  
>>> prefer to
>>> have a set of rules that would apply to the general case.
>>>
>>> Right now I would prefer to devote my time to other aspects of  
>>> lmer -
>>> in particular I am still working on code for generalized linear  
>>> mixed
>>> models using a supernodal Cholesky factorization.  I am willing  
>>> to put
>>> this aside and code up the rules for denominator degrees of freedom
>>> with nested grouping factors BUT first I want someone to show me an
>>> example demonstrating that there really is a problem.  The example
>>> must show that the p-value calculated in the anova table or the
>>> parameter estimates table for lmer is seriously wrong compared to an
>>> empirical p-value - obtained from simulation under the null
>>> distribution or through MCMC sampling or something like that.  
>>> Saying
>>> that "Software XYZ says there are n denominator d.f. and lmer says
>>> there are m" does NOT count as an example.  I will readily concede
>>> that the denominator degrees of freedom reported by lmer are  
>>> wrong but
>>> so are the degrees of freedom reported by Software XYZ because there
>>> is no right answer (in general - in a few simple balanced designs
>>> there may be a right answer).
>>>
>>>>
>>>>> m.lm <- lm(Yvar~Xvar)
>>>>>
>>>>> anova(m.lm)
>>>> Analysis of Variance Table
>>>>
>>>> Response: Nadultos
>>>>             Df  Sum Sq Mean Sq F value Pr(>F)
>>>> Xvar         1    33.6    33.6  1.2352 0.2667
>>>> Residuals  878 23896.2    27.2
>>>>>
>>>>> summary(m.lm)
>>>>
>>>> Call:
>>>> lm(formula = Yvar ~ Xvar)
>>>>
>>>> Residuals:
>>>>     Min      1Q  Median      3Q     Max
>>>> -2.7455 -2.3545 -1.7455  0.2545 69.6455
>>>>
>>>> Coefficients:
>>>>                Estimate Std. Error t value Pr(>|t|)
>>>> (Intercept)      2.3545     0.2487   9.467   <2e-16 ***
>>>> XvarFactor2      0.3909     0.3517   1.111    0.267
>>>> ---
>>>> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>>>>
>>>> Residual standard error: 5.217 on 878 degrees of freedom
>>>> Multiple R-Squared: 0.001405,   Adjusted R-squared: 0.0002675
>>>> F-statistic: 1.235 on 1 and 878 DF,  p-value: 0.2667
>>>>
>>>> I read the rnews about this use of the full DF in lmer, but I dont
>>>> undestand
>>>> this use with a gaussian error, I undestand this with glm data.
>>>>
>>>> I need more explanations, please.
>>>>
>>>> Thanks
>>>> Ronaldo
>>
>> ______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide! http://www.R-project.org/posting- 
>> guide.html
>>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html