Mixed Models: Random Effects

Hello there, friend! Welcome to the most random Shiny App you've ever visited.

Purpose:

The purpose of this Shiny app is to visualize the changes that happen to your mixed-effects models. Specifically, to visualize the residuals and the standard error of the fixed effects when you add random effects to a mixed-effects model.

Audience:

This Shiny app is for you if you have ever fit a mixed-effects model and have had difficulty imagining what is happening to the explained variance of your model as you add different parameters. It can be used to understand your own models but also as a teaching tool to help others understand the theoretical principles behind adding random effect parameters to your mixed-effects model.

Stability:

You will see that as you click around on the random effect parameters, they will be added or removed from the model, which in turn has an effect on the residuals and standard error.

A guide to our beautiful Shiny App:

The Simulate Button: This button will simulate data based on our generative (correct) model.
The Random Effects: These checkboxes allow you to select the random effects you would like to add to the model you are estimating.
The Residual Plot Tab: A demonstration of the residuals from the generative (correct) model (black dots) compared to the residuals from the model you are estimating (red dots).
The Forest Plot Tab: A demonstration of the fixed effects (pink squares) and their confidence intervals (blue lines). This plot also contains a vertical line at 0 to see if your fixed effects are significant.

A suggested order of operations:

Simulate data based on the generative model
Choose which random effects you would like to include in your model
Pick a plot to view
Change around the random effects you chose originally
See the dynamic magic before your eyes!

Estimated residuals from current model: ⠄ Estimated residuals from the generative (correct) model:⠄

The Residual Plot: A demonstration of the residuals from the generative (correct) model (black dots) compared to the residuals from the model you are estimating (red dots). The x-axis is broken up into chunks representing a subset of the participants in the study. Within each chunk, the residuals associated with each item completed can be found. The y-axis, with 0 in the middle, demonstrates the size of the residual.

Note: In mixed-effects models, you would like to reduce the amount of residuals (i.e., have them closest to 0) because you would like to explain as much variance as possible (i.e., have as little unexplained variance as possible).

Note some more: In mixed-effects models, you want the residuals to be completely independent from one another, as this is an assumption that is made when running this model. In other words, you do not want to see any pattern in the residual plot.

Note even more: In mixed-effects models, it is possible to overfit your data such that you become overconfident in the model you have estimated.

The Forest Plot: A demonstration of the fixed effects (pink squares) and their confidence intervals (blue lines). This plot also contains a vertical line at 0 to see if your fixed effects are significant.

Note: As the model becomes complex, the confidence intervals become larger. This is because when models are too simple, we become overly confident about the precision of our estimates (i.e., have smaller confidence intervals).

This app was created during a TquanT seminar by Kelsey MacKay (Leuven), Hannah Rós Sigurðardóttir (Amsterdam), Johanna Xemaire (Tübingen) and Cristina Mendonça (Lisbon).

shinyServer(function(input, output) { library("simgen") library("ggplot2") library("sweetalertR") # install with #devtools::install_github("timelyportfolio/sweetalertR") library("forestplot") library("lme4") observeEvent(input$sim, ignoreNULL = TRUE, { # #####BIG RED SIMULATE BUTTON####### # ################################## # ## define generative parameters: this will be sent to mkDf() as the argument 'mcr.params' v <- c(int = 0, eff = 0, # intercept + slope (fixed effects) err = 1, # residual error miss = 0, pMin = 0, pMax = 0, # no missing data t00 = 10, t11 = 10, rsub = 0, # by-subj random slopes and intercepts (unit) w00 = 1, w11 = 0, ritm = 0) # by-item random intercept (unit) ## only thing you should never change are: rsub and ritm. Must be zero otherwise you can't ## independently scale the item and slope variances. ## simulate some data (mkDf is a function in the 'simgen' package my_data <- mkDf(nsubj = 100, nitem = 10, # small N for visualization mcr.params = v, # generative parameters wsbi = TRUE, # Within Subjects Between Item (WSBI) design verbose = TRUE) # return random effects for subjects/items # ## the 'verbose = TRUE' option returns additional generative information in the form of a list, ## with the following elements: # ## # ## my_data$dat : data # ## my_data$subj_re : subject random effects # ## my_data$item_re : item random effects # ## my_data$err : model residuals # ## create a different 'condition' variable my_data$dat$cond <- factor(ifelse(my_data$dat$Cond == -.5, "A", "B")) # # ######THE FULL MODEL (FROM WHICH WE SIMULATE THE DATA)############# # ################################################################## # ## fitting the (maximal random effects) model mod <- lmer(Resp ~ Cond + (Cond || SubjID) + (1 | ItemID), data = my_data$dat, REML = FALSE) summary(mod) real.res <- resid(mod) # To use for residual plot function # # ################## Plot function ######################### # residual plot res.plot <- function(real.res=real.res, model, no.items=4, no.persons=4){ validate( need(input$randomSettings != "", HTML('ERROR: Please set at least one random factor.'))) c<-no.persons*no.items real.res<-real.res[1:c] est.res <- resid(model) est.res=est.res[1:c] items<-rep(1:no.items, no.persons) it.s <- paste(rep(1:no.persons, each=no.items), items, sep=".") labels<-paste("Person", 1:no.persons, sep=" ") breaks <- as.numeric(paste(1:no.persons, 1, sep=".")) qplot(it.s, real.res, xlab="Items per person", ylab="Residuals")+geom_point(aes(it.s, est.res), color="red")+ scale_x_discrete("breaks"=breaks, labels=labels)+ylim(-3,3)+ggtitle("Residuals") } # forest plot library("forestplot") forest <- function(model){ validate( need(input$randomSettings != "", HTML('ERROR: Please set at least one random factor.'))) int.est<-summary(model)$coefficients[1,1] slope.est <- summary(model)$coefficients[2,1] int.se <- summary(model)$coefficients[1,2] slope.se <- summary(model)$coefficients[2,2] forestplot(c("Intercept", "Slope"), mean=c(int.est, slope.est), lower=c(int.est-int.se,slope.est-slope.se), upper=c(int.est+int.se,slope.est+slope.se), col=fpColors(box="pink", lines="blue", text="black", axes="black"), xlab="Estimated values",boxsize=0.05, graphwidth = unit(x = 5, units = "inches")) } #######SLOPE AND INTERCEPT COMBINATIONS################### values <- reactiveValues() obsB <- observe({ test <- paste(unlist(input$randomSettings), collapse=' ') if (test == ""){ # No random intercepts, model is estimated to avoid a crash values$equation <- HTML('ERROR: Please set at least one random factor.') values$model <- lmer(Resp ~ Cond + (1 | SubjID), data = my_data$dat, REML = FALSE) }else if(test == "SUBI"){ #1##random subject intercept only values$model <- lmer(Resp ~ Cond + (1 | SubjID), data = my_data$dat, REML = FALSE) values$equation <- "$$Y_{si} = \\beta_0 + S_{0s} + \\beta_1X_i + e_{si}$$" } else if (test == "SUBS"){ #2##random subject slope only values$model <- lmer(Resp ~ Cond + (0 + Cond || SubjID), data = my_data$dat, REML = FALSE) values$equation <- "$$Y_{si} = \\beta_0 + (\\beta_1 + S_{1s})X_i + e_{si}$$" } else if (test == "ITI"){ #3#random item intercept only values$model <- lmer(Resp ~ Cond + (1 | ItemID), data = my_data$dat, REML = FALSE) values$equation <- "$$Y_{si} = \\beta_0 + I_{0i} + \\beta_1X_i + e_{si}$$" } else if (test == "ITS"){ #4#random item slope only values$model <- lmer(Resp ~ Cond + (0 + Cond | ItemID), data = my_data$dat, REML = FALSE) values$equation <- "$$Y_{si} = \\beta_0 + (\\beta_1 + I_{1s})X_i + e_{si}$$" } else if (test == "SUBI SUBS") { #5#subject intercept and slope only values$model <- lmer(Resp ~ Cond + (Cond || SubjID), data = my_data$dat, REML = FALSE) values$equation <- "$$Y_{si} = \\beta_0 + S_{0s} + (\\beta_1 + S_{1s})X_i + e_{si}$$" } else if (test == "ITI ITS") { #6#item intercept and slope only values$model <- lmer(Resp ~ Cond + (1 + Cond | ItemID), data = my_data$dat, REML = FALSE) values$equation <- "$$Y_{si} = \\beta_0 + I_{0i} + (\\beta_1 + I_{1s})X_i + e_{si}$$" } else if (test == "SUBI ITI") { #7#subject intercept and item intercept values$model <- lmer(Resp ~ Cond + (1 | SubjID) + (1 | ItemID), data = my_data$dat, REML = FALSE) values$equation <- "$$Y_{si} = \\beta_0 + S_{0s} + I_{0i} + \\beta_1X_i + e_{si}$$" } else if (test == "SUBS ITS") { #8#subject slope and item slope values$model <- lmer(Resp ~ Cond + (0 + Cond || SubjID) + (0 + Cond | ItemID), data = my_data$dat, REML = FALSE) values$equation <- "$$Y_{si} = \\beta_0 + (\\beta_1 + S_{1s} + I_{1s})X_i + e_{si}$$" } else if (test == "SUBS ITI") { #9#item intercept and subject slope values$model <- lmer(Resp ~ Cond + (0 + Cond || SubjID) + (1 | ItemID), data = my_data$dat, REML = FALSE) values$equation <- "$$Y_{si} = \\beta_0 + I_{0i} + (\\beta_1 + S_{1s})X_i + e_{si}$$" } else if (test == "SUBI ITS") { #10#subject intercept and item slope values$model <- lmer(Resp ~ Cond + (1 | SubjID) + (0+Cond | ItemID), data = my_data$dat, REML = FALSE) values$equation <- "$$Y_{si} = \\beta_0 + S_{0s} + (\\beta_1 + I_{1s})X_i + e_{si}$$" } else if (test == "SUBI SUBS ITI") { #11#subject intercept and slope and item intercept values$model <- lmer(Resp ~ Cond + (Cond || SubjID) + (1 | ItemID), data = my_data$dat, REML = FALSE) values$equation <- "$$Y_{si} = \\beta_0 + S_{0s} + I_{0i} +(\\beta_1 + S_{1s})X_i + e_{si}$$" } else if (test == "SUBI SUBS ITS") { #12#subject intercept and slope and item slope values$model <- lmer(Resp ~ Cond + (Cond || SubjID) + (0+Cond | ItemID), data = my_data$dat, REML = FALSE) values$equation <- "$$Y_{si} = \\beta_0 + S_{0s} + (\\beta_1 + S_{1s} + I_{1s})X_i + e_{si}$$" } else if (test == "SUBI ITI ITS") { #13#subject intercept and item intercept and slope values$model <- lmer(Resp ~ Cond + (1 | SubjID) + (1 + Cond | ItemID), data = my_data$dat, REML = FALSE) values$equation <- "$$Y_{si} = \\beta_0 + S_{0s} + I_{0i} + (\\beta_1 + I_{1s})X_i + e_{si}$$" } else if (test == "SUBS ITI ITS") { #14#subject slope and item intercept and slope values$model <- lmer(Resp ~ Cond + (0+Cond || SubjID) + (1 + Cond| ItemID), data = my_data$dat, REML = FALSE) values$equation <- "$$Y_{si} = \\beta_0 + I_{0i} + (\\beta_1 + S_{1s} + I_{1s})X_i + e_{si}$$" } else if (test == "SUBI SUBS ITI ITS") { #15#subject intercept and slope and item intercept and slope - ALL IN BABY values$model <- lmer(Resp ~ Cond + (Cond || SubjID) + (Cond | ItemID), data = my_data$dat, REML = FALSE) values$equation <- "$$Y_{si} = \\beta_0 + S_{0s} + I_{0i} + (\\beta_1 + S_{1s} + I_{1s})X_i + e_{si}$$"} if(any(grepl("failed to converge", values$model@optinfo$conv$lme4$messages))) { values$warning <- HTML('ERROR: Model failed to converge! Please resimulate data.') }else if(any(grepl("Hessian is numerically singular", values$model@optinfo$conv$lme4$messages))) { values$warning <- HTML('ERROR: Model Hessian is single! Please resimulate data.') }else if(any(grepl("nearly unidentifiable", values$model@optinfo$conv$lme4$messages))){ values$warning <- HTML('ERROR: Model is nearly unidentifiable! Please resimulate data.') }else if(!any(grepl("failed to converge", values$model@optinfo$conv$lme4$messages))) { values$warning <- "" }}) output$text <- renderUI({ withMathJax( values$equation ) }) output$warning <- renderUI({ values$warning }) output$plot.res <- renderPlot({ res.plot(real.res=real.res, model=values$model) }) output$plot.forest <- renderPlot({ forest(values$model) }) }) }) shinyApp(ui=shinyUI, server=shinyServer)

library(sweetalertR) shinyUI(fluidPage( includeScript("../../../Matomo-tquant.js"), titlePanel("Mixed Models: Random Effects"), sidebarLayout( sidebarPanel( br(), HTML(""), HTML(' <button id="sim" type="button" class="btn btn-danger action-button btn-lg"> Simulate </button> '), HTML(""), sweetalert(selector = "#sim", text = "Your simulation was successful!", title = "Yay :D", type = "success"), br(), withMathJax(), "We used this super cool model to generate the data: $$Y_{si} = \\beta_0 + S_{0s} + I_{0i} + (\\beta_1 + S_{1s})X_i + e_{si}$$", p("And this is the model you are currently estimating:"), uiOutput("text"), br(), checkboxGroupInput(inputId = "randomSettings", label = "Define Random Effects", choices=list("Random Subject Intercept (\$S_{0s}\$)" = "SUBI", "Random Subject Slope (\$S_{1s}\$)" = "SUBS", "Random Item Intercept (\$I_{0i}\$)" = "ITI", "Random Item Slope (\$I_{1i}\$)" = "ITS"), selected="SUBI"), br(), uiOutput("warning"), width = 4 ), mainPanel( tabsetPanel( tabPanel("Welcome Page", HTML(" Hello there, friend! Welcome to the most random Shiny App you've ever visited. Purpose: The purpose of this Shiny app is to visualize the changes that happen to your mixed-effects models. Specifically, to visualize the residuals and the standard error of the fixed effects when you add random effects to a mixed-effects model. Audience: This Shiny app is for you if you have ever fit a mixed-effects model and have had difficulty imagining what is happening to the explained variance of your model as you add different parameters. It can be used to understand your own models but also as a teaching tool to help others understand the theoretical principles behind adding random effect parameters to your mixed-effects model. Stability: You will see that as you click around on the random effect parameters, they will be added or removed from the model, which in turn has an effect on the residuals and standard error. A guide to our beautiful Shiny App: <ul> <li>The Simulate Button: This button will simulate data based on our generative (correct) model.</il> <li>The Random Effects: These checkboxes allow you to select the random effects you would like to add to the model you are estimating.</li> <li>The Residual Plot Tab: A demonstration of the residuals from the generative (correct) model (black dots) compared to the residuals from the model you are estimating (red dots).</li> <li>The Forest Plot Tab: A demonstration of the fixed effects (pink squares) and their confidence intervals (blue lines). This plot also contains a vertical line at 0 to see if your fixed effects are significant.</li> </ul> A suggested order of operations: <ol> <li>Simulate data based on the generative model</li> <li>Choose which random effects you would like to include in your model</li> <li>Pick a plot to view</li> <li>Change around the random effects you chose originally</li> <li>See the dynamic magic before your eyes!</li></ol>")), tabPanel("Residual Plot", plotOutput("plot.res"), HTML(' Estimated residuals from current model: ⠄ Estimated residuals from the generative (correct) model:⠄ '), HTML("The Residual Plot: A demonstration of the residuals from the generative (correct) model (black dots) compared to the residuals from the model you are estimating (red dots). The x-axis is broken up into chunks representing a subset of the participants in the study. Within each chunk, the residuals associated with each item completed can be found. The y-axis, with 0 in the middle, demonstrates the size of the residual. "), br(), HTML("Note: In mixed-effects models, you would like to reduce the amount of residuals (i.e., have them closest to 0) because you would like to explain as much variance as possible (i.e., have as little unexplained variance as possible)."), HTML("Note some more: In mixed-effects models, you want the residuals to be completely independent from one another, as this is an assumption that is made when running this model. In other words, you do not want to see any pattern in the residual plot."), HTML("Note even more: In mixed-effects models, it is possible to overfit your data such that you become overconfident in the model you have estimated.")), tabPanel("Forest Plot", plotOutput("plot.forest"), HTML("The Forest Plot: A demonstration of the fixed effects (pink squares) and their confidence intervals (blue lines). This plot also contains a vertical line at 0 to see if your fixed effects are significant."), br(), HTML("Note: As the model becomes complex, the confidence intervals become larger. This is because when models are too simple, we become overly confident about the precision of our estimates (i.e., have smaller confidence intervals).")), tabPanel("The Team", HTML(""), br(),br(), img(src="pic.png", height = 334, width = 400), br(),br(), HTML("This app was created during a TquanT seminar by Kelsey MacKay (Leuven), Hannah Rós Sigur&#240ardóttir (Amsterdam), Johanna Xemaire (Tübingen) and Cristina Mendonça (Lisbon)."), HTML("")) ) ))) )