[R-sig-ME] Fwd: same old question - lme4 and p-values
Robert Kushler
kushler at oakland.edu
Sat Apr 5 17:50:18 CEST 2008
Unfortunately, the latest changes to lme4 seem to have broken languageR.
Many of the plotting functions are now built in to lme4, but the p-value
obsessed among us would love to have the other functions in languageR
back in action.
Harold, any plans for an update? Will Doug let you do it? :-)
Rob Kushler
Andy Fugard wrote:
> I made (more) sense of mixed effects models when I went back to standard
> regression and ANOVAs with the philosophy: Everything is a Comparison.
>
> So for instance noting different ways of getting the magical numbers
> that result from doing regression:
>
> ---8<--------------------------------------------------------------------
>
> > m.0 = lm(Fertility ~ 1, data = swiss)
> > m.full = lm(Fertility ~ ., data = swiss)
> > summary(m.full)
>
> ...
>
> Residual standard error: 7.17 on 41 degrees of freedom
> Multiple R-squared: 0.707, Adjusted R-squared: 0.671
> F-statistic: 19.8 on 5 and 41 DF, p-value: 5.59e-10
> ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
>
> > anova(m.0,m.full)
> Analysis of Variance Table
>
> Model 1: Fertility ~ 1
> Model 2: Fertility ~ Agriculture + Examination + Education + Catholic +
> Infant.Mortality
> Res.Df RSS Df Sum of Sq F Pr(>F)
> 1 46 7178
> 2 41 2105 5 5073 19.8 5.6e-10 ***
> ^^ ^ ^^^^^^^^^^^^
> ----------------------------------------------------------------->8------
>
> Noting what happens when you compare nested models (incidentally, a
> confusing term when used the context of "multilevel" models):
>
> ---8<--------------------------------------------------------------------
>
> > m.full = lm(Fertility ~ . , data = swiss)
> > summary(m.full)
>
> ...
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 66.9152 10.7060 6.25 1.9e-07 ***
> Agriculture -0.1721 0.0703 -2.45 0.0187 *
> Examination -0.2580 0.2539 -1.02 0.3155
> Education -0.8709 0.1830 -4.76 2.4e-05 ***
> Catholic 0.1041 0.0353 2.95 0.0052 **
> Infant.Mortality 1.0770 0.3817 2.82 0.0073 **
> ^^^^^^
> ...
>
> >
> > m1 = update(full.model, ~. -Infant.Mortality)
> > anova(m.full,m1)
> Analysis of Variance Table
>
> Model 1: Fertility ~ Agriculture + Examination + Education + Catholic +
> Infant.Mortality
> Model 2: Fertility ~ Agriculture + Examination + Education + Catholic
> Res.Df RSS Df Sum of Sq F Pr(>F)
> 1 41 2105
> 2 42 2514 -1 -409 7.96 0.0073 **
> ^^^^^^
>
> ----------------------------------------------------------------->8------
>
> Playing around with this sort of thing made it a lot easier to
> understand why statisticians get annoyed (e.g.,
> <http://www.stats.ox.ac.uk/pub/MASS3/Exegeses.pdf> for a good example)
> with us non-statisticians.
>
> On a more practical note, the languageR package has an especially useful
> pvals.fnc function for when you really do need a test of whether a slope
> is significantly different to zero, e.g. when you've got a load of
> categorical predictors and want to see where the difference is.
>
> Hope that wasn't too far off topic.
>
> Cheers,
>
> Andy
>
>
>
> Hank Stevens wrote:
>> Google:
>> p-values lmer wiki
>>
>> On Apr 4, 2008, at 9:33 AM, Douglas Bates wrote:
>>
>>> ---------- Forwarded message ----------
>>> From: Douglas Bates <bates at stat.wisc.edu>
>>> Date: Fri, Apr 4, 2008 at 7:54 AM
>>> Subject: Re: same old question - lme4 and p-values
>>> To: andreas.nord at zooekol.lu.se
>>>
>>>
>>>
>>> On Fri, Apr 4, 2008 at 5:24 AM, <andreas.nord at zooekol.lu.se> wrote:
>>>> Dear Prof. Bates,
>>>> I've recently switched to using R for my analyses, and I find the
>>> lme4 package to be extremely helpful. I have read your explanation
>>> (posted on the mailing list) of why you choose not to display
>>> p-values. Unfortunately, most of the journals I publish in require
>>> that I include p-values, which is why I have to find a way of
>>> calculating them from the lmer output. However, not being a trained
>>> statistician I have some difficulties following your recommendations
>>> given in the explanatory text. In other words, after having fitted my
>>> model, I am not at all sure on what to do in order to obtain p-values
>>> (or similar).
>>>
>>>> I am sorry to have to bother you with a question I know you have
>>> already answered many times, but perhaps you would be so kind as to
>>> give me some hints on how to proceed.
>>>
>>> I understand your situation. Statisticians have created the "every
>>> question of scientific interest must be answered by a p-value" monster
>>> and now it turns on us. Nevertheless I am reluctant to give advice on
>>> p-values in lme4 because apparently I don't know how to do it
>>> correctly.
>>>
>>> May I send a copy of this reply to the
>>> R-SIG-Mixed-Models at r-project.org mailing list? ("SIG" == "Special
>>> Interest Group")? (I ask your permission to send the copy because I am
>>> quoting your original question.) Some who subscribe to that mailing
>>> list may have the courage to wade into this swamp and offer their
>>> advice.
>>>
>>> _______________________________________________
>>> R-sig-mixed-models at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>
>> Dr. Hank Stevens, Associate Professor
>> 338 Pearson Hall
>> Botany Department
>> Miami University
>> Oxford, OH 45056
>>
>> Office: (513) 529-4206
>> Lab: (513) 529-4262
>> FAX: (513) 529-4243
>> http://www.cas.muohio.edu/~stevenmh/
>> http://www.cas.muohio.edu/ecology
>> http://www.muohio.edu/botany/
>>
>> "If the stars should appear one night in a thousand years, how would men
>> believe and adore." -Ralph Waldo Emerson, writer and philosopher
>> (1803-1882)
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>
>
More information about the R-sig-mixed-models
mailing list