Mixed-effect model selection with aictab

if (!require("pacman")) install.packages("pacman")
pacman::p_load(s20x, lme4, AICcmodavg, MASS)

# Data prep
data(Orthodont,package="nlme")
head(Orthodont)

## Grouped Data: distance ~ age | Subject
##   distance age Subject  Sex
## 1     26.0   8     M01 Male
## 2     25.0  10     M01 Male
## 3     29.0  12     M01 Male
## 4     31.0  14     M01 Male
## 5     21.5   8     M02 Male
## 6     22.5  10     M02 Male

# Model definition
null <- lmer(distance ~ 1 + (1|Subject), data=Orthodont, REML=FALSE)
age <- lmer(distance ~ age + (1|Subject), data=Orthodont, REML=FALSE)
sex <- lmer(distance ~ Sex + (1|Subject), data=Orthodont, REML=FALSE)
add <- lmer(distance ~ age + Sex + (1|Subject), data=Orthodont, REML=FALSE)
# Model selection
 cand.mod.names <- c("null", "age", "sex", "add") 
cand.mods <- list( ) 
      
       # This function fills the list by model names
          for(i in 1:length(cand.mod.names)) {
            cand.mods[[i]] <- get(cand.mod.names[i]) }
      
      # Function aictab does the AICc-based model comparison
       print(aictab(cand.set = cand.mods, 
                    modnames = cand.mod.names))

## 
## Model selection based on AICc:
## 
##      K   AICc Delta_AICc AICcWt Cum.Wt      LL
## add  5 445.44       0.00   0.96   0.96 -217.43
## age  4 451.78       6.33   0.04   1.00 -221.69
## sex  4 515.35      69.90   0.00   1.00 -253.48
## null 3 521.72      76.28   0.00   1.00 -257.75