Comparision of lmms and lms based on visual comparison, as well as deviance of each model

You can specify the parameters of the dataset you want to generate. For example you you set the fixed intercept to 8, wich means that the mean of Y will be 8. Random effects are represented by their standardeviation. If you choose to select a correlation of 0 for the random intercept and the random slope, you see them as uncorrelated. Dont input any negative standarddeviation Q&A is on it. Probably

How many subjects?

How many items?

Which fixed intercept do you want to use?

Which fixed slope do you want to use?

What standarddeviation you want to use for the distribution of the random intercept effects (subjects)?

What standarddeviation you want to use for the distribution of the random slope effects (subjects)?

What correlation between intercept and slope you want to use (subjects)?

What standarddeviation you want to use for the distribution of the error terms?

Generated dataset

Model comparision Chisq-Test Deviance

Is the difference in deviance significant? (p-value)

show with app

data_gen <- function(nsubj, nitem, beta.0, beta.1, int.sd, slope.sd, r.corr, int.sd.it, slope.sd.it, r.corr.it, error.sd){
  
  subj<-data.frame(subjid=1:nsubj,
                   gender=sample(c("M","F"), nsubj,replace=TRUE))
  item<-data.frame(itemid=1:nitem,
                   cond=rep(c(-0.5, 0.5), # this is basically effect coding
                            each=nitem/2))
  trials<-expand.grid(subjid=1:nsubj,
                      itemid=1:nitem)
  dat1<-merge(subj,trials, "subjid")
  dat2<-merge(dat1,item,"itemid")
  
  vcovmx<-matrix(c(int.sd^2,
                   r.corr*int.sd*slope.sd,
                   r.corr*int.sd*slope.sd,
                   slope.sd^2), nrow=2, ncol=2)
  vcovitem <- matrix(c(int.sd.it^2,
                       r.corr.it*int.sd.it*slope.sd.it,
                       r.corr.it*int.sd.it*slope.sd.it,
                       slope.sd.it^2), nrow=2, ncol=2)
  subjrandom <- mvrnorm(nsubj, c(0, 0), vcovmx)
  itemrandom <- mvrnorm(nitem, c(0, 0), vcovitem)
  subj$rint<-subjrandom[,1]
  subj$rslo<-subjrandom[,2]
  item$rint.it <- itemrandom[, 1]
  item$rslo.it <- itemrandom[, 2]
  
  #put dataset together
  dat1<-merge(subj, trials, "subjid")
  dat2<-merge(dat1, item, "itemid")
  dat2$err<-rnorm(nsubj * nitem, sd=error.sd )
  
  #equation
  dat2$Y <- beta.0 + dat2$rint + dat2$rint.it + (beta.1 + dat2$rslo + dat2$rslo.it)*dat2$cond + dat2$err
  
  dat2  
}

#### server.R
## implementation of comparison of different linear mixed models
## using shiny, on a fixed design
##
## Team:
##    Agnes Altmanninger, Graz
##    Alexandra Martinez Rodriguez, Madrid
##    Dominik Frommeyer, Tuebingen
##    Eline Belmans, Leuven
##    Filipa Almeida, Lisbon
##    Katarzyna Bzymek, Leuven
##
## last modified, 08.03.2016, dF

## to do list:
#
# Loading of datasets (!) (obsolete)
#
# Generating dataset  (!)
#
# Previewing datasets (!)
#
# Specify Model parameters ()
#   - just lm   (!)
# model-type
#       glmm or lmm        (obs)
#       family, in case of glmm (obs)
#       link, in case of glmm   (obs)
#   - fixed effects from dataframe      (!)
#   - dependent variable from dataframe (obs)
#   - random effects -"-                (!)
#   - Display model
#       - Having 2 to 3 outputs available   (!)
#         - summary,                        (!)
#         - prediction (xyplot?),           ()
#         - modeltest                       (?)
#

library(shiny)
library(lme4)
library(MASS)
library(ggplot2)
source("data_gen.R")
#install.packages("devtools")
#devtools::install_github("dalejbarr/simgen")
#library(simgen)

shinyServer(function(input, output){
  ## getting data
  dat <- reactiveValues()
  observeEvent(input$go_comp, {
#      dat$file <- mkDf(input$nsubj, input$nitem, input$check_within,
#                       c(input$beta = 3,   ## fixed intercept (beta_0)
#                         eff = 6,   ## fixed effect of condition (beta 1)
#                         err = 10,  ## residual SD
#                         t00 = 5,   ## by-subject random intercept (SD)
#                         t11 = 3,   ## by-subject random slope (SD)
#                         rsub = .5, ## by-subject random correlation (-1 to 1)
#                         w00 = 5,   ## by-item random intercept
#                         w11 = 3,   ## by-item random slope
#                         ritm = .5, ## by-item random correlation
#                         ## the params below aren't really necessary to alter
#                         miss = 0,  ## proportion of data that are missing
#                         pMin = 0,  ## minimum proportion of obs missing
#                         pMax = 0)  ## maximum proportion of obs missing))
       dat$file <- data_gen(input$nsubj,
                       input$nitem,
                       input$beta0,
                       input$beta1,
                       abs(input$intsd),
                       abs(input$slopesd),
                       input$rcorr, 
                       0, #input$intsdit,
                       0, #input$slopesdit,
                       0, #input$rcorrit,
                       input$errorsd)
       # print(trigger$model)
      output$go_model_choice <- renderUI({
        actionButton("go_model_choice", "Display selected Model")
      })
     })
  #read.table("paldata.txt", header=TRUE)
  # AV = Response
  # UV = Age*GoogleCoocFreq 
  # random = (1 + Age|Pair)
  trigger <- reactiveValues()
  trigger$model <- 1

  output$preview_data <- renderPrint({
    print(
      head(dat$file
           #[order(dat$file$subjid,
            #                  dat$file$itemid),]
               )
    , digits=2)
  })
  output$summary <- renderPrint({
    str(dat$file)
  })
  
  
  
  # (x/y)%%1==0 checks if there are fractions left
  observeEvent(input$go_model_choice, {
    trigger$model <- trigger$model + 1
    
    dat$mod1 <- lm(Y ~ 1, dat$file)
    dat$mod2 <- lm(Y ~ 1 + cond, dat$file)
    dat$mod3 <- lmer(Y ~ 1 + cond + (1 | subjid), 
                     dat$file, REML=input$choose_method)
    dat$mod4 <- lmer(Y ~ 1 + cond + (1 + cond | subjid),
                     dat$file, REML=input$choose_method)
    dat$mod5 <- lmer(Y ~ (1 | subjid),
                     dat$file, REML=input$choose_method)
#     dat$mod5 <- lmer(Y ~ 1 + cond + (1 | itemid),
#                      dat$file, REML=input$choose_method)
    dat$mod6 <- lmer(Y ~ 1 + cond + (0 + cond | subjid),
                     dat$file, REML=input$choose_method)
#     dat$mod6 <- lmer(Y ~ 1 + cond + (1 + cond | itemid),
#                      dat$file, REML=input$choose_method)
#     dat$mod7 <- lmer(Y ~ 1 + cond + (1 | subjid) + (1 | itemid),
#                      dat$file, REML=input$choose_method)
#     dat$mod8 <- lmer(Y ~ 1 + cond + (1 + cond | itemid) + (1 | subjid),
#                      dat$file, REML=input$choose_method)
#     dat$mod9 <- lmer(Y ~ 1 + cond + (1 + cond | subjid) + (1 | itemid),
#                      dat$file, REML=input$choose_method)
    #     }
#     dat$mod10 <- lmer(Y ~ 1 + cond + (1 + cond | subjid) + (1 + cond | itemid),

    dat$mod0 <- dat[[input$choose_model]]
    dat$mod0.predict <- predict(dat[[input$choose_model]])
    # (x/y)%%1==0
    if((trigger$model/2) == floor(trigger$model/2)){
      dat$mod.panel1 <- dat$mod0
      output$summary_model1 <- renderPrint({
        print(summary(isolate(dat$mod0)), digits=1)
      })
      
#       output$model1_converg <- renderPrint({
#         dat$mod0$converge
#       })
      
      output$pred_model1 <- renderPlot({
        # predict(lmer2, data.frame(subj_id=rep(1:2, each=10), item_id=rep(1:10, times=2), cond=-0.5))
        # predicts for the first two participants for each item the values, in condition -0.5
        ggplot(isolate(dat$file), aes(cond, Y, group=interaction(subjid), col=subjid)) + 
          geom_line(aes(y=isolate(dat$mod0.predict)), size=0.8) +
          geom_point(alpha = 0.3) + 
          geom_hline(yintercept=0, linetype="dashed") +
          theme_bw()
      })
    }
    else{
      dat$mod.panel2 <- dat$mod0
#       output$model2_converge <- renderPrint({
#         dat$mod0$converge
#       })
      output$summary_model2 <- renderPrint({
        print(summary(dat$mod0), digits=1)
      })
      
      output$pred_model2 <- renderPlot({
        # predict(lmer2, data.frame(subj_id=rep(1:2, each=10), item_id=rep(1:10, times=2), cond=-0.5))
        # predicts for the first two participants for each item the values, in condition -0.5
        ggplot(isolate(dat$file), aes(cond, Y, group=interaction(subjid), col=subjid)) + 
          geom_line(aes(y=isolate(dat$mod0.predict)), size=0.8) +
          geom_point(alpha = 0.3) + 
          geom_hline(yintercept=0, linetype="dashed") +
          theme_bw()
      })
    }
    
    output$test_model1 <- renderPrint({
      if(is.null(dat$mod.panel1)){
        return("Select a model")
      }
      dat$mod1.deviance <- deviance(dat$mod.panel1)
      dat$mod1.deviance
      #logLik(dat$mod.panel1, REML=input$choose_method)
    })
    output$test_model2 <- renderPrint({
      if(is.null(dat$mod.panel2)){
        return("Select a model")
      }
      dat$mod2.deviance <- deviance(dat$mod.panel2)
      dat$mod2.deviance
      #logLik(dat$mod.panel2, REML=input$choose_method)
    })
    
    output$custom.chisq <- renderPrint({
      if(is.null(dat$mod.panel2) | is.null(dat$mod.panel1)){
        "Select two models"
      }
      else{
        pval <- pchisq(abs(dat$mod1.deviance - dat$mod2.deviance), abs(attr(logLik(dat$mod.panel1), "df") - attr(logLik(dat$mod.panel2), "df")), lower.tail = FALSE)
        pval
      }
      
   })
    
#     output$pred_ranef_model1 <- renderPlot({
#       if(is.null(dat$mod.panel1)){
#         return (NULL)
#       }
#       else{
#         if(class(dat$mod.panel1)=="lm"){
#           if(is.na(coef(dat$mod.panel1)["(Intercept)"])){
#             plot(rep(0,
#                      times=input$nsubj*input$nitem) ~ dat$file$rint,
#                  ylab="Estimated random intercept/Subj",
#                  xlab="Generated random intercept/Subj")
#           }
#           plot(rep(as.vector(coef(dat$mod.panel1)["(Intercept)"]),
#                    times=isolate(input$nsubj*input$nitem)) ~ dat$file$rint,
#                ylab="Estimated random intercept/Subj",
#                xlab="Generated random intercept/Subj")
#         }
#       else{
#         if(is.null(ranef(dat$mod.panel1)$subjid["(Intercept)"])){
#           plot(rep(as.vector(t(coef(dat$mod.panel1)$subjid["(Intercept)"]))
#                    , times=isolate(input$nsubj*input$nitem)) ~ dat$file$rint,
#                ylab="Estimated random intercept/Subj",
#                xlab="Generated random intercept/Subj")
#         }
#         else{
#           plot(rep(as.vector(t(ranef(dat$mod.panel1)$subjid["(Intercept)"]))
#                    , times=isolate(input$nitem)) ~ dat$file$rint,
#                ylab="Estimated random intercept/Subj",
#                xlab="Generated random intercept/Subj")
#         }
#       }
#       }
#     })
#     output$pred_ranef_model2 <- renderPlot({
#       if(is.null(dat$mod.panel2)){
#         return (NULL)
#       }
#       else{
#         if(class(dat$mod.panel2)=="lm"){
#           if(is.na(coef(dat$mod.panel2)["(Intercept)"])){
#             plot(rep(0,
#                      times=input$nsubj*input$nitem) ~ dat$file$rint,
#                  ylab="Estimated random intercept/Subj",
#                  xlab="Generated random intercept/Subj")
#           }
#           else{
#             plot(rep(as.vector(coef(dat$mod.panel2)["(Intercept)"]),
#                      times=isolate(input$nsubj*input$nitem)) ~ dat$file$rint,
#                  ylab="Estimated random intercept/Subj",
#                  xlab="Generated random intercept/Subj")
#             }
#         }
#         else{
#           if(is.null(ranef(dat$mod.panel2)$subjid["(Intercept)"])){
#             plot(rep(as.vector(t(coef(dat$mod.panel2)$subjid["(Intercept)"]))
#                      , times=isolate(input$nsubj*input$nitem)) ~ dat$file$rint,
#                  ylab="Estimated random intercept/Subj",
#                  xlab="Generated random intercept/Subj")
#           }
#           else{
#             plot(rep(as.vector(t(ranef(dat$mod.panel2)$subjid["(Intercept)"]))
#                      , times=isolate(input$nitem)) ~ dat$file$rint,
#                  ylab="Estimated random intercept/Subj",
#                  xlab="Generated random intercept/Subj")
#           }
#         }
#       }
#     })
# #    print(as.vector(t(coef(dat$mod.panel1)$subjid[, "cond"])))
#     output$pred_ranef_slope.1 <- renderPlot({
#       if(is.null(dat$mod.panel1)){
#         return (NULL)
#       }
#       else{
#         if(class(dat$mod.panel1)=="lm"){
#           if(is.na(coef(dat$mod.panel1)["cond"])){
#             plot(rep(0,
#                      times=input$nsubj*input$nitem) ~ dat$file$rslo,
#                  ylab="Estimated random slope/Subj",
#                  xlab="Generated random slope/Subj")
#           }
#           else{
#             plot(rep(as.vector(coef(dat$mod.panel1)["cond"]),
#                      times=isolate(input$nsubj*input$nitem)) ~ dat$file$rslo,
#                  ylab="Estimated random slope/Subj",
#                  xlab="Generated random slope/Subj")
#           }
#         }
#         else{
#           if(is.null(ranef(dat$mod.panel1)$subjid[,"cond"])){
#             plot(rep(as.vector(t(coef(dat$mod.panel1)$subjid[, "cond"])),
#                      times=isolate(input$nsubj*input$nitem)) ~ dat$file$rslo,
#                  ylab="Estimated random slope/Subj",
#                  xlab="Generated random slope/Subj")
#           }
#           else{
#             plot(rep(as.vector(t(ranef(dat$mod.panel1)$subjid[, "cond"])),
#                      times=isolate(input$nitem)) ~ dat$file$rslo,
#                  ylab="Estimated random slope/Subj",
#                  xlab="Generated random slope/Subj")
#           }
#         }
#       }
#     })
#     output$pred_ranef_slope.2 <- renderPlot({
#       if(is.null(dat$mod.panel2)){
#         return (NULL)
#       }
#       else{
#         if(class(dat$mod.panel2)=="lm"){
#           if(is.na(coef(dat$mod.panel2)["cond"])){
#             plot(rep(0,
#                      times=input$nsubj*input$nitem) ~ dat$file$rslo,
#                  ylab="Estimated random slope/Subj",
#                  xlab="Generated random slope/Subj")
#           }
#           else{
#             plot(rep(as.vector(coef(dat$mod.panel2)["cond"]),
#                      times=isolate(input$nsubj*input$nitem)) ~ dat$file$rslo,
#                  ylab="Estimated random slope/Subj",
#                  xlab="Generated random slope/Subj")
#           }
#         }
#         else{
#           if(is.null(ranef(dat$mod.panel2)$subjid["cond"])){
#             plot(rep(as.vector(t(coef(dat$mod.panel2)$subjid["cond"])),
#                      times=isolate(input$nsubj*input$nitem)) ~ dat$file$rslo,
#                  ylab="Estimated random slope/Subj",
#                  xlab="Generated random slope/Subj")
#           }
#           else{
#             plot(rep(as.vector(t(ranef(dat$mod.panel2)$subjid["cond"])),
#                      times=isolate(input$nitem)) ~ dat$file$rslo,
#                  ylab="Estimated random slope/Subj",
#                  xlab="Generated random slope/Subj")
#           }
#         }
#       }
#     })
  })
})

#### ui.R
## implementation of comparison of different linear mixed models
## using shiny, on a fixed design
##
## Team:
##    Agnes Altmanninger, Graz
##    Alexandra Martinez Rodriguez, Madrid
##    Dominik Frommeyer, Tuebingen
##    Eline Belmans, Leuven
##    Filipa Almeida, Lisbon
##    Katarzyna Bzymek, Leuven
##
## last modified, 08.03.2016, dF

## to do list:
#
# Loading of datasets (!) (obsolete)
#
# Generating dataset  ()
#
# Previewing datasets (!) 
#
# Specify Model parameters ()
#   - just lm   (!)
# model-type
#       glmm or lmm        (obs)
#       family, in case of glmm (obs)
#       link, in case of glmm   (obs)
#   - fixed effects from dataframe      (!)
#   - dependent variable from dataframe (obs)
#   - random effects -"-                (!)
#   - Display model
#       - Having 2 to 3 outputs available   (!)
#         - summary,                        (!)
#         - prediction (xyplot?),           (?)
#         - modeltest                       (?)
#

library(shiny)

shinyUI(fluidPage(
  includeScript("../../../Matomo-tquant.js"),
  titlePanel("Comparision of lmms and lms based on visual comparison, as well as deviance of each model"),
  p("You can specify the parameters of the dataset you want to generate. For example you you set the fixed intercept to 8, wich means that the mean of Y will be 8. Random effects are represented by their standardeviation. If you choose to select a correlation of 0 for the random intercept and the random slope, you see them as uncorrelated. Dont input any negative standarddeviation Q&A is on it. Probably"),
#     column(4, p(fileInput("choose_file", "Choose dataset"),
#                 checkboxInput("choose_header", "Header?"),
#                 radioButtons("choose_sep", "Give seperation value",
#                              c(Comma=",",
#                                Semicolon=";",
#                                Tab="\t",
#                                Space=" "),
#                              ','),
#                 selectInput("choose_dec", "Decimalsymbol", 
#                             choices=c(".", ","), selected="."),
#                 actionButton("go_data", "Upload selected file"))),
      fluidRow(column(6, sliderInput("nsubj", 
                  label = c("How many subjects?"), 
                  value = 10,
                  min=10, max=50)),
      column(6, sliderInput("nitem", 
                  label = c("How many items?"), 
                  value = 6,
                  step=2, min=6, max=20))),
      fluidRow(column(6, numericInput("beta0", 
                   label = c("Which fixed intercept do you want to use?"), 
                   value = 1)),
      column(6, numericInput("beta1", 
                   label = c("Which fixed slope do you want to use?"), 
                   value = 1))),
      fluidRow(column(4, numericInput("intsd", 
                   label = c("What standarddeviation you want to use for the distribution of the random intercept effects (subjects)?"), 
                   value = 1, min=0)),
      column(4, numericInput("slopesd", 
                   label = c("What standarddeviation you want to use for the distribution of the random slope effects (subjects)?"), 
                   value = 1, min=0)),
      column(4, sliderInput("rcorr", 
                  label = c("What correlation between intercept and slope you want to use (subjects)?"), 
                  value = 0.3, min=-1, max=1, step=0.05))),
#       fluidRow(column(4, numericInput("intsdit",
#                    "What standarddeviation do you want to us for the distribution of the random slope (items)",
#                    value=1, min=0)),
#       column(4, numericInput("slopesdit", 
#                    "What standarddeviation you want to use for the distribution of the random slope effects (items)?",
#                    value=1, min=0)),
#       column(4, sliderInput("rcorrit",
#                   "What correlation between intercept and slope you want to use (items)?",
#                   value=0.3, min=-1, max=1, step=0.05))),
      fluidRow(column(4, numericInput("errorsd", 
                   label = c("What standarddeviation you want to use for the distribution of the error terms?"), 
                   value = 1, min=0)),
      column(8, #checkboxInput("check_within", "Is cond within subject?"),
      actionButton("go_comp", "Generate data"))),
   fluidRow(column(6, verbatimTextOutput("preview_data")),
            column(6, h2("Generated dataset"))),
#    verbatimTextOutput("summary"),
  sidebarLayout(
    sidebarPanel(width=3,
      ## specify the models for each case, then have radio Buttons for those 
      ## models
      "Choose the model",
      radioButtons("choose_model", "",
                   choices=c("Fixed Intercept"="mod1",
                             "Fixed Slope|Intercept"="mod2",
                             "Fixed Slope|Intercept + Random Intercept/Subj"="mod3",
                             "Fixed Slope|Intercept + Random Slope|Intercept/Subj"="mod4",
                             "Random Intercept/Subj"="mod5",
                             "Fixed Slope|Intercept + Random Slope/Subj"="mod6"
#                              "Fixed Slope|Intercept + Random Slope|Intercept/Item"="mod6",
#                              "Fixed Slope|Intercept + Random Intercept/Subj|Item"="mod7",
#                              "Fixed Slope|Intercept + Random Slope|Intercept/Item + Intercept/Subj"="mod8",
#                              "Fixed Slope|Intercept + Random Slope|Intercept/Subj + Intercept/Item"="mod9",
#                              "Fixed Slope|Intercept + Random Slope|Intercept/Subj|Item"="mod10"
                             )),

      checkboxInput("choose_method", "Use REML?"),
      uiOutput("go_model_choice")
#       radioButtons("choose_type", "Choose model type",
#                    choices=c("lmer", "glmer")),
#       uiOutput("gen_family"),
#       radioButtons("choose_cdng", "Choose coding of effects",
#                    choices=c("dummy", "effect")),
#       radioButtons("choose_Y", "Choose dependent variable"),
#       uiOutput("gen_UV")
    ),
    mainPanel(width=9,
      fluidRow(
        column(6, p(verbatimTextOutput("summary_model1"),
                    plotOutput("pred_model1")
                    )),
        column(6, p(verbatimTextOutput("summary_model2"),
                    plotOutput("pred_model2")
                    ))),
      p("Model comparision Chisq-Test",
        "Deviance",
        verbatimTextOutput("test_model1"), verbatimTextOutput("test_model2")
        ,
         "Is the difference in deviance significant? (p-value)",
         verbatimTextOutput("custom.chisq")
        )
#      ,
#         "Random intercepts",
#       fluidRow(column(6, plotOutput("pred_ranef_model1")
#                       ),
#                column(6, plotOutput("pred_ranef_model2")
#                       )),
#       "Random Slopes",
#       fluidRow(column(6, plotOutput("pred_ranef_slope.1") 
#                       ),
#                column(6, plotOutput("pred_ranef_slope.2")
#                       ))
      )
    )
))