17

I have a mixed effects model and I would like to see the R²- and p-value. I thought this is acessible by summary() but it's not. Or at least I don't recognize it.

> summary(fit1.lme <- lmer(log(log(Amplification)) ~ poly(Voltage, 3) + (1 | Serial_number), data = bdf))
Linear mixed model fit by REML ['lmerMod']
Formula: log(log(Amplification)) ~ poly(Voltage, 3) + (1 | Serial_number)
   Data: bdf

REML criterion at convergence: -253237.6

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-14.8183  -0.4863  -0.0681   0.2941   9.3292 

Random effects:
 Groups        Name        Variance Std.Dev.
 Serial_number (Intercept) 0.008435 0.09184 
 Residual                  0.001985 0.04456 
Number of obs: 76914, groups:  Serial_number, 1270

Fixed effects:
                    Estimate Std. Error t value
(Intercept)         0.826745   0.002582     320
poly(Voltage, 3)1 286.978430   0.045248    6342
poly(Voltage, 3)2 -74.061993   0.045846   -1615
poly(Voltage, 3)3  39.605454   0.045505     870

Correlation of Fixed Effects:
            (Intr) p(V,3)1 p(V,3)2
ply(Vlt,3)1 0.001                 
ply(Vlt,3)2 0.002  0.021          
ply(Vlt,3)3 0.001  0.032   0.028  
Ben
  • 1,432
  • 4
  • 20
  • 43
  • 2
    (1) `?pvalues`, http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#testing-hypotheses; (2) http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#how-do-i-compute-a-coefficient-of-determination-r2-or-an-analogue-for-glmms – Ben Bolker Jul 26 '17 at 12:42
  • I'm still not used to the fact that I can get help from my code :) Thanks! – Ben Jul 26 '17 at 12:46

3 Answers3

31

For the R², you can use r.squaredGLMM(fit1.lme) from the MuMIn package. It will returns the marginal and the conditional R².

For the p-value, you can find them by using summary with the lmerTest package.

For more information on p-values with mixed models: http://mindingthebrain.blogspot.ch/2014/02/three-ways-to-get-parameter-specific-p.html

Azor Ahai -him-
  • 123
  • 1
  • 7
abichat
  • 2,317
  • 2
  • 21
  • 39
13

I add a very small demos with hierarchial modeling for ozone layer where the modeling acknowledges that it varies by month. You can find comparisons below. I could find the R squared term only in MuMIn package.

MuMIn package

> data(airquality)

> MuMIn::r.squaredGLMM(lme4::lmer(data=airquality, Ozone ~ 1 + (1|Month)))
     R2m       R2c
[1,]   0 0.2390012
> summary(lm(data=airquality, Ozone ~ 1 + (1|Month)))$r.squared
[1] 0

where we compare the linear regression and the mixed effect model aka hierarchial regression model.

Linear regression

> summary(lm(data=airquality, Ozone ~ 1 + (1|Month)))

Call:
lm(formula = Ozone ~ 1 + (1 | Month), data = airquality)

Residuals:
   Min     1Q Median     3Q    Max 
-41.13 -24.13 -10.63  21.12 125.87 

Coefficients: (1 not defined because of singularities)
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)     42.129      3.063   13.76   <2e-16 ***
1 | MonthTRUE       NA         NA      NA       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 32.99 on 115 degrees of freedom
  (37 observations deleted due to missingness)

lmer4

> summary(lme4::lmer(data=airquality, Ozone ~ 1 + (1|Month)))
Linear mixed model fit by REML ['lmerMod']
Formula: Ozone ~ 1 + (1 | Month)
   Data: airquality

REML criterion at convergence: 1116.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.7084 -0.6269 -0.2669  0.4121  3.7507 

Random effects:
 Groups   Name        Variance Std.Dev.
 Month    (Intercept) 270.6    16.45   
 Residual             861.6    29.35   
Number of obs: 116, groups:  Month, 5

Fixed effects:
            Estimate Std. Error t value
(Intercept)   41.093      7.922   5.187

lmerTest

library(lmerTest)

> lmerTest::lmer(data=airquality, Ozone ~ 1 + (1|Month))
Linear mixed model fit by REML ['lmerModLmerTest']
Formula: Ozone ~ 1 + (1 | Month)
   Data: airquality
REML criterion at convergence: 1116.544
Random effects:
 Groups   Name        Std.Dev.
 Month    (Intercept) 16.45   
 Residual             29.35   
Number of obs: 116, groups:  Month, 5
Fixed Effects:
(Intercept)  
      41.09  
> summary(lmerTest::lmer(data=airquality, Ozone ~ 1 + (1|Month)))
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Ozone ~ 1 + (1 | Month)
   Data: airquality

REML criterion at convergence: 1116.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.7084 -0.6269 -0.2669  0.4121  3.7507 

Random effects:
 Groups   Name        Variance Std.Dev.
 Month    (Intercept) 270.6    16.45   
 Residual             861.6    29.35   
Number of obs: 116, groups:  Month, 5

Fixed effects:
            Estimate Std. Error     df t value Pr(>|t|)   
(Intercept)   41.093      7.922  4.096   5.187  0.00616 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hhh
  • 50,788
  • 62
  • 179
  • 282
5

You can try packages sjPlot or sjstats. The first package help create APA-style tables from lme4 analyses, and the second package for extracting fit statistics.

You just need to simply write code:

tab_model(fit1.lme)

It will output an APA table including estimates slope, intercept, CI, p-value, variances, residuals, observation number, ICC, marginal and conditional R square etc.

Looks like this: enter image description here

Issa Chi
  • 61
  • 1
  • 3
  • If I remembered correct, this plot is from sjPlot, and sjstat can help extract p value with the code: p_value(). I personally prefer use sjPlot, because it is more clear and return more statistical values. Hope it helps. – Issa Chi Nov 06 '21 at 21:48