I have a built a mixed effect model using the lmer()
function from the lme4
package. The lme4
package does not output the p-value of the coefficients for some good philosophical reason. However, I still need p-values to report in my publication. I know there are multiple ways to calculate p-values using the model created by lmer()
, e.g. here.
My problem is: I want to extract p-value using the tidy()
function from the broom
package. Here, I really want to stick with tidy()
because I want to maintain the following pipeline:
data_frame %>% group_by(grouping variables) %>% do(tidy(fitted_model))
One option would be to create a custom function and append it to the pipeline. However, the man page of the broom
package (http://rpackages.ianhowson.com/cran/broom/man/lme4_tidiers.html) says:
"p.value P-value computed from t-statistic (may be missing/NA)".
By this I am assuming a function to calculate p-value from the t-value given by lmer model has already been implemented in broom. So, I am reluctant to reinvent the wheel.
The problem is I don't get the column with name p.value at all. I was expecting a column named p.value with NAs as the worst case scenario.
Code:
library(lme4)
library(broom)
lme <- lmer(Reaction ~ Days + (1 + Days | Subject), sleepstudy)
tidy(lme)
tidy(lme, effects = "fixed")
Output:
> tidy(lme)
term estimate std.error statistic group
1 (Intercept) 251.40510485 6.824557 36.838306 fixed
2 Days 10.46728596 1.545789 6.771485 fixed
3 sd_(Intercept).Subject 24.74045195 NA NA Subject
4 sd_Days.Subject 5.92213312 NA NA Subject
5 cor_(Intercept).Days.Subject 0.06555113 NA NA Subject
6 sd_Observation.Residual 25.59181564 NA NA Residual
> tidy(lme, effects = "fixed")
term estimate std.error statistic
1 (Intercept) 251.40510 6.824557 36.838306
2 Days 10.46729 1.545789 6.771485
You will need package lmerTest
to obtain p-values. tidy
will not work on the lme
object and you will need to append it to your format.
attach(mtcars)
lme <- lmer(mpg ~ cyl + (1 + cyl | carb), mtcars)
summary(lme)