[R-sig-ME] A consultation about DF of the result of lmer
Fox, John
j|ox @end|ng |rom mcm@@ter@c@
Wed May 8 18:04:31 CEST 2019
Dear Rumeng He,
First, when you ask a question on the R-sig-ME list, it's polite to copy further messages to the list, as I've done with this response to your latest message.
Second, you don't say what model you fit to the sleepstudy data, so I'll answer by ESP:
Here I fit two models with lmer() to the sleepstudy data, one treating Days as a numeric predictor and one as a factor. You'll see that I get 1 df for the term in the first case and 9 in the second, as one would expect, which is what Ben Bolker and I both suggested:
--------------- snip ------------------
> library(car)
Loading required package: carData
> library(lme4)
Loading required package: Matrix
Registered S3 methods overwritten by 'lme4':
method from
cooks.distance.influence.merMod car
influence.merMod car
dfbeta.influence.merMod car
dfbetas.influence.merMod car
> # example from ?lmer
> (fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy))
Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (Days | Subject)
Data: sleepstudy
REML criterion at convergence: 1743.628
Random effects:
Groups Name Std.Dev. Corr
Subject (Intercept) 24.737
Days 5.923 0.07
Residual 25.592
Number of obs: 180, groups: Subject, 18
Fixed Effects:
(Intercept) Days
251.41 10.47
> Anova(fm1)
Analysis of Deviance Table (Type II Wald chisquare tests)
Response: Reaction
Chisq Df Pr(>Chisq)
Days 45.843 1 1.281e-11 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> # making Days a factor
> sleepstudy$DaysF <- as.factor(sleepstudy$Days)
> levels(sleepstudy$DaysF)
[1] "0" "1" "2" "3" "4" "5" "6" "7" "8" "9"
> # can only fit random-intercept model
> (fm1F <- lmer(Reaction ~ DaysF + (1 | Subject), sleepstudy))
Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ DaysF + (1 | Subject)
Data: sleepstudy
REML criterion at convergence: 1729.493
Random effects:
Groups Name Std.Dev.
Subject (Intercept) 37.09
Residual 31.43
Number of obs: 180, groups: Subject, 18
Fixed Effects:
(Intercept) DaysF1 DaysF2 DaysF3 DaysF4 DaysF5 DaysF6 DaysF7 DaysF8 DaysF9
256.652 7.844 8.710 26.340 31.998 51.867 55.526 62.099 79.978 94.199
> Anova(fm1F)
Analysis of Deviance Table (Type II Wald chisquare tests)
Response: Reaction
Chisq Df Pr(>Chisq)
DaysF 168.32 9 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> Anova(fm1F, test.statistic="F")
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: Reaction
F Df Df.res Pr(>F)
DaysF 18.703 9 153 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------- snip ---------------
So what's the problem?
John
> -----Original Message-----
> From: 何如梦 [mailto:18754808835 using 163.com]
> Sent: Wednesday, May 8, 2019 2:14 AM
> To: Fox, John <jfox using mcmaster.ca>
> Subject: Re:RE: [R-sig-ME] A consultation about DF of the result of lmer
>
> Dear John,
> Thanks for your reply. I tried to translate random effect as factors
> but it didn't work. (The dataset of sleepstudy is in R)
>
>
> Then I also made analyse of "aov". The Df is right. So I think there
> would be other reasons for this wrong.
>
> I also found other pepole met the same question (as the follow website)
> but for too many years ago. I want to know whethere I can correct
> https://stat.ethz.ch/pipermail/r-help/2008-October/176084.html
>
>
> At 2019-05-08 01:49:39, "Fox, John" <jfox using mcmaster.ca> wrote:
> >Dear Rumeng He,
> >
> >For what it's worth, my guess is the same as Ben's. If we're wrong,
> then you'll have to send more information about what you did, ideally
> including a reproducible example of the problem.
> >
> >Best,
> > John
> >
> >--------------------------------------
> >John Fox, Professor Emeritus
> >McMaster University
> >Hamilton, Ontario, Canada
> >Web: socialsciences.mcmaster.ca/jfox/
> >
> >
> >
> >> -----Original Message-----
> >> From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces using r-
> >> project.org] On Behalf Of Ben Bolker
> >> Sent: Tuesday, May 7, 2019 12:31 PM
> >> To: r-sig-mixed-models using r-project.org
> >> Subject: Re: [R-sig-ME] A consultation about DF of the result of lmer
> >>
> >> It's very hard to say without more information (try e.g. here
> >> <https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-
> >> reproducible-example>
> >> or look at other posts in the list archive
> >> <https://stat.ethz.ch/pipermail/r-sig-mixed-models/>
> >>
> >> A first guess is that you have a fixed-effect predictor that's
> >> supposed to be a factor with 6 levels, but is actually being
> >> interpreted as a numeric variable. If that's the case, then either
> >> changing it within the data set
> >>
> >> mydata$pred1 <- factor(mydata$pred1)
> >>
> >> or doing it on the fly in your model
> >>
> >> lmer(response ~ factor(pred1) + ...)
> >>
> >> should fix the problem.
> >>
> >>
> >> On 2019-05-05 9:11 p.m., 何如梦 wrote:
> >> > Dear,
> >> > I am a graduated student who's topic is ecology. Recently, I am
> >> studying how to establish a liner mixed model to exclude the error
> >> caused by the difference of site by using lme4 in R. I found when I
> >> calculated p value by using the function of "Anova" of car package
> >> the Df is 1. But in fact the data set has 6 levels. Then I also
> >> operate according to the code of reference PDF of lme4 in P52(lmer)
> >> without any change. When run "anova(fm1, fm2)" I found the Df of fm1
> >> also is 1. As I see if the Df is wrong the p value would be wrong
> >> either. I want to know how to correct my wrong. I will very
> appreciate your reply.
> >> >
> >> >
> >> > Best regards,
> >> > Rumeng He
> >> > [[alternative HTML version deleted]]
> >> >
> >> > _______________________________________________
> >> > R-sig-mixed-models using r-project.org mailing list
> >> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >> >
> >>
> >> _______________________________________________
> >> R-sig-mixed-models using 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