There are four related functions for mixed effects analyses: mixed_model, mixed_anova, mixed_model_slopes, and mixed_anova_slopes.

mixed_model(data, Y_value, Fixed_Factor, Random_Factor, ...)

Arguments

data

a data table object, e.g. data.frame or tibble.

Y_value

name of column containing quantitative (dependent) variable, provided within "quotes".

Fixed_Factor

name(s) of categorical fixed factors (independent variables) provided as a vector if more than one or within "quotes".

Random_Factor

name(s) of random factors to allow random intercepts; to be provided as a vector when more than one or within "quotes".

...

any additional arguments to pass on to lmer if required.

Value

This function returns an S4 object of class "lmerModLmerTest".

Details

This function uses lmer to fit a linear mixed effect model and provides the model object, which could be used for post-hoc comparisons. The model object is converted to class lmerModLmerTest object by as_lmerModLmerTest.

It requires a data table, one dependent variable (Y_value), one or more independent variables (Fixed_Factor), and at least one random factor (Random_Factor). These should match names of variables in the long-format data table exactly. This function is related to mixed_anova. Output of this function can be used with posthoc_Pairwise, posthoc_Levelwise and posthoc_vsRef, or with emmeans.

More than one fixed factors can be provided as a vector (e.g. c("A", "B")). A full model with interaction term is fitted. This means when Y_value = Y, Fixed_factor = c("A", "B"), Random_factor = "R" are entered as arguments, these are passed on as Y ~ A*B + (1|R) (which is equivalent to Y ~ A + B + A:B + (1|R)). For simplicity, only random intercepts are fitted ((1|R)).

Also see mixed_anova_slopes and mixed_model_slopes for similar functions where variable slopes and intercept models are fit.

Examples

#one fixed factor and random factor
mixed_model(data = data_doubling_time, 
Y_value = "Doubling_time", 
Fixed_Factor = "Student", 
Random_Factor = "Experiment")
#> boundary (singular) fit: see help('isSingular')
#> Linear mixed model fit by REML ['lmerModLmerTest']
#> Formula: Doubling_time ~ Student + (1 | Experiment)
#>    Data: data_doubling_time
#> REML criterion at convergence: 95.2499
#> Random effects:
#>  Groups     Name        Std.Dev.
#>  Experiment (Intercept) 0.000   
#>  Residual               1.989   
#> Number of obs: 30, groups:  Experiment, 3
#> Fixed Effects:
#> (Intercept)     StudentB     StudentC     StudentD     StudentE     StudentF  
#>     19.9619      -0.3713      -0.5929      -1.1728      -0.6286       0.7416  
#>    StudentG     StudentH     StudentI     StudentJ  
#>      0.4885      -0.8083       0.6775       1.2115  
#> optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings 

#two fixed factors as a vector, one random factor
mixed_model(data = data_cholesterol, 
Y_value = "Cholesterol", 
Fixed_Factor = c("Treatment", "Hospital"), 
Random_Factor = "Subject")
#> Linear mixed model fit by REML ['lmerModLmerTest']
#> Formula: Cholesterol ~ Treatment * Hospital + (1 | Subject)
#>    Data: data_cholesterol
#> REML criterion at convergence: 376.5107
#> Random effects:
#>  Groups   Name        Std.Dev.
#>  Subject  (Intercept) 51.756  
#>  Residual              6.524  
#> Number of obs: 50, groups:  Subject, 25
#> Fixed Effects:
#>                         (Intercept)                 TreatmentBefore_drug  
#>                             147.011                               26.165  
#>                      HospitalHosp_b                       HospitalHosp_c  
#>                              -4.587                              -42.004  
#>                      HospitalHosp_d                       HospitalHosp_e  
#>                              15.770                               15.432  
#> TreatmentBefore_drug:HospitalHosp_b  TreatmentBefore_drug:HospitalHosp_c  
#>                              -8.075                               -7.001  
#> TreatmentBefore_drug:HospitalHosp_d  TreatmentBefore_drug:HospitalHosp_e  
#>                             -27.202                              -25.803  

#save model
model <- mixed_model(data = data_doubling_time, 
Y_value =  "Doubling_time", 
Fixed_Factor = "Student", 
Random_Factor = "Experiment")
#> boundary (singular) fit: see help('isSingular')

#get model summary
summary(model)
#> Linear mixed model fit by REML. t-tests use Satterthwaite's method [
#> lmerModLmerTest]
#> Formula: Doubling_time ~ Student + (1 | Experiment)
#>    Data: data_doubling_time
#> 
#> REML criterion at convergence: 95.2
#> 
#> Scaled residuals: 
#>      Min       1Q   Median       3Q      Max 
#> -1.94563 -0.40680 -0.04097  0.42601  1.45893 
#> 
#> Random effects:
#>  Groups     Name        Variance Std.Dev.
#>  Experiment (Intercept) 0.000    0.000   
#>  Residual               3.956    1.989   
#> Number of obs: 30, groups:  Experiment, 3
#> 
#> Fixed effects:
#>             Estimate Std. Error      df t value Pr(>|t|)    
#> (Intercept)  19.9619     1.1484 20.0000  17.383 1.54e-13 ***
#> StudentB     -0.3713     1.6241 20.0000  -0.229    0.821    
#> StudentC     -0.5929     1.6241 20.0000  -0.365    0.719    
#> StudentD     -1.1728     1.6241 20.0000  -0.722    0.479    
#> StudentE     -0.6286     1.6241 20.0000  -0.387    0.703    
#> StudentF      0.7416     1.6241 20.0000   0.457    0.653    
#> StudentG      0.4885     1.6241 20.0000   0.301    0.767    
#> StudentH     -0.8083     1.6241 20.0000  -0.498    0.624    
#> StudentI      0.6775     1.6241 20.0000   0.417    0.681    
#> StudentJ      1.2115     1.6241 20.0000   0.746    0.464    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Correlation of Fixed Effects:
#>          (Intr) StdntB StdntC StdntD StdntE StdntF StdntG StdntH StdntI
#> StudentB -0.707                                                        
#> StudentC -0.707  0.500                                                 
#> StudentD -0.707  0.500  0.500                                          
#> StudentE -0.707  0.500  0.500  0.500                                   
#> StudentF -0.707  0.500  0.500  0.500  0.500                            
#> StudentG -0.707  0.500  0.500  0.500  0.500  0.500                     
#> StudentH -0.707  0.500  0.500  0.500  0.500  0.500  0.500              
#> StudentI -0.707  0.500  0.500  0.500  0.500  0.500  0.500  0.500       
#> StudentJ -0.707  0.500  0.500  0.500  0.500  0.500  0.500  0.500  0.500
#> optimizer (nloptwrap) convergence code: 0 (OK)
#> boundary (singular) fit: see help('isSingular')
#>