Introduction to LearnVizLMM

Why without data?

• To expand the types of examples used in statistics education.
• To inform study design (e.g., model and variable selection prior to data collection or model fitting).

Installation

To install this package, you can use CRAN (the central R package repository).

library(LearnVizLMM)

Background

In many settings, multiple samples are collected from the same individual or group (e.g., achievement scores of students from different schools and classrooms within schools; impact of different diets on the weights of individuals sampled once a month for six months). This type of data is often referred to as multilevel and can be analyzed using LMMs. LMMs can be fit in R using the following packages and functions.

library(lme4)
lme4::lmer(formula = outcome ~ fixed_effects + random_effects)

library(nlme)
nlme::lme(fixed = outcome ~ fixed_effects, random = random_effects)

The outcome and fixed effects are similar to that of a linear model.

lmer(formula = Y ~ 1 + X1 + X2 + random_effects)
lme(fixed = Y ~ 1 + X1 + X2, random = random_effects)

The random effects are specified with respect to groups or grouping factors (GFs), which are factor variables under which observations are nested or grouped. Using different notations, users can specify random intercepts, random slopes, and more to be estimated by the model.

# Random Intercept
lmer(formula = outcome ~ fixed_effects + (1|GF))
lme(fixed = outcome ~ fixed_effects, random = ~1|GF)
lme(fixed = outcome ~ fixed_effects, random = list(GF=~1))
# Random Intercept & Slope
lmer(formula = outcome ~ fixed_effects + (1 + X1|GF))
lme(fixed = outcome ~ fixed_effects, random = ~1 + X1|GF)
lme(fixed = outcome ~ fixed_effects, random = list(GF=~1 + X1))

output_type = “latex” (default)

extract_equation(model = "lmer(score ~ 1 + age + (1|subject))")

$$
\begin{aligned}
  \operatorname{score}_{ij} &= \beta_{0} + \beta_{1}(\operatorname{age}) \\
  &+ u_{0i} \\
 &+ \epsilon_{ij} \\
 u_{0i} &\sim N \left(0, \tau^2_{u_{0}} \right), \text{ for subject i = 1,} \dots \text{,a} \\
 \epsilon_{ij} &\sim N \left(0, \sigma^2 \right), \text{ for Observation j = 1,} \dots \text{,n}  
\end{aligned}
$$

output_type = “string”

extract_equation(model = "lmer(score ~ 1 + age + (1|subject))",
                 output_type = "string")

[1] "$$\n\\begin{aligned}\n  \\operatorname{score}_{ij} &= \\beta_{0} + \\beta_{1}(\\operatorname{age}) \\\\\n  &+ u_{0i} \\\\\n &+ \\epsilon_{ij} \\\\\n u_{0i} &\\sim N \\left(0, \\tau^2_{u_{0}} \\right), \\text{ for subject i = 1,} \\dots \\text{,a} \\\\\n \\epsilon_{ij} &\\sim N \\left(0, \\sigma^2 \\right), \\text{ for Observation j = 1,} \\dots \\text{,n}  \n\\end{aligned}\n$$\n"

output_type = “none”

extract_equation(model = "lmer(score ~ 1 + age + (1|subject))",
                 output_type = "none")

One Group, Random Intercept and Slopes

By default, all predictor variables are assumed to be numeric.

extract_variables(model = "lme(weight~1+sex+time+I(time^2), random=~1+time+I(time^2)|ID)")

Two Groups, Random Intercepts and Slope

To include one or more categorical predictors, use cat_vars and cat_vars_nlevels. In addition, fixed and random intercepts are added to the data frame unless explicitly excluded (e.g., “age-1” and “0+age”).

extract_variables(model = "lme(Score~type+age*sex,random=list(School=pdDiag(~type),Class=~1))",
                  cat_vars = c("type", "sex"),
                  cat_vars_nlevels = c(3, 2))

One group

extract_structure(n_gf = 1, 
                  gf_names = "Subject")

To specify the number of Subjects, use gf_nlevels.

extract_structure(model = "lme(Weight ~ Time, random=~Time|Subject, data)",
                  gf_nlevels = 47)

To remove the labels for the levels on the left side of the image, use label_levels.

extract_structure(model = "lme(Weight ~ Time, random=~Time|Subject, data)",
                  gf_nlevels = "m", 
                  label_levels = "no")

Two groups

When there are two groups, they can either be crossed (e.g., every level of Worker is observed for every level of Machine)

extract_structure(model = "lmer(Strength ~ 1 + (1|Machine) + (1|Worker))",
                  gf_nlevels = c(10, 20))

or nested (e.g., class is nested within school).

extract_structure(model = "lme(score~type, random=list(school=pdDiag(~1+type),class=~1))")

These figures can also be created using the following code.

extract_structure(n_gf = 2, 
                  gf_description = "crossed",
                  gf_names = c("Machine", "Worker"),
                  gf_nlevels = c(10, 20))

extract_structure(n_gf = 2, 
                  gf_description = "nested",
                  gf_names = c("school", "class"))

Three groups

When there are three groups, they can also be crossed or nested.

extract_structure(n_gf = 3,
                  gf_description = "crossed",
                  gf_names = c("District", "School", "Class"),
                  gf_nlevels = c(8, 15, 5))

The index used for the highest-level group or the group that entered the model or gf_names first can be left as "i" (see above) or redefined using gf3_index.

extract_structure(n_gf = 3,
                  gf_description = "nested",
                  gf_names = c("District", "School", "Class"),
                  gf_nlevels = c(8, 15, 5),
                  gf3_index = "q")

There may also be a combination of nested and crossed groups. For example, one group may be crossed with two nested groups (i.e. crossed with nested).

extract_structure(n_gf = 3,
                  gf_description = "crossed with nested",
                  gf_names = c("GF1", "GF2", "GF3"))

Another example is when two crossed groups are nested within another group (i.e. crossed within nested).

extract_structure(n_gf = 3,
                  gf_description = "crossed within nested",
                  gf_names = c("GF1", "GF2", "GF3"))

Export Type

If you wish to edit the output beyond what is offered by extract_structure(), set export_type = "text". After editing the text, the figure can be created using grViz() in the DiagrammeR package.

diagram_text <- extract_structure(n_gf = 1, 
                                  gf_names = "Subject",
                                  export_type = "text")

DiagrammeR::grViz(diagram_text)