library(shiny)
library(DT)
library(dplyr)
library(markdown)
library(simgen)
library(stargazer)
library(lme4)
#Shiny server function
shinyServer(function(input,output,session){
dat <- eventReactive(input$go,{
prams <- c(int = input$int,
eff = input$eff,
err = input$err,
t00 = input$t00,
t11 = input$t11,
rsub = input$rsub,
w00 = input$w00,
w11 = input$w11,
ritm = input$ritm,
miss = 0,
pMin = 0,
pMax = 0
)
is_between_item <- input$C2
nsubj <- input$nsubj
nitem <- input$nitem
data <- mkDf(nsubj, nitem, is_between_item, prams)
## make simulated data
data <- mkDf(nsubj, nitem, is_between_item, prams)
return(data)
})
output$tbl = DT::renderDataTable(dat(),
server = TRUE,
options = list(dom = 'C<"clear">lfrtip'))
output$downloadData <- downloadHandler(
filename = paste0("simulated_data_", Sys.Date(),".csv"),
content = function(file) {
write.csv(dat(),file,row.names=F)
}
)
#The models
#Model with by-subject random intercept + slope and by-item random intercept
mod <- eventReactive(input$go,{
# if ("Model with by-subject random intercept + slope and by-item random intercept" %in% input$model_types){
inut_data <- dat()
mf <- (Resp ~ Cond + (1 + Cond | SubjID) + (1 | ItemID))
mod <- lmer(mf, inut_data, REML=FALSE)
return(mod)
# }
# else{
# return(NULL)
# }
})
output$mod_1 <- renderPrint({
summary(mod())
})
#Model with by-subject random intercept and by-item random intercept
mod_no_slope <- eventReactive(input$go,{
# if ("Model with by-subject random intercept and by-item random intercept" %in% input$model_types){
inut_data <- dat()
mf_no_slope <- (Resp ~ Cond + (1| SubjID) + (1 | ItemID))
mod <- lmer(mf_no_slope, inut_data, REML=FALSE)
return(mod)
# }
#
# else{
# return(NULL)
# }
})
output$mod_2 <- renderPrint({
summary(mod_no_slope())
})
#Model with only fixed effects
mod_no_random <- eventReactive(input$go,{
# if ("Model with only fixed effects" %in% input$model_types){
inut_data <- dat()
mf_no_random <- (Resp ~ Cond)
mod <- lm(mf_no_random, inut_data)
return(mod)
# }
# else{
# return(NULL)
# }
})
output$mod_3 <- renderPrint({
summary(mod_no_random())
})
#Model with by-item intercept and slope
mod_overfitted <- eventReactive(input$go,{
# if ("Model with by-item intercept and slope" %in% input$model_types){
inut_data <- dat()
mf_overfitted <- (Resp ~ Cond + (1 + Cond | SubjID) + (1 + Cond | ItemID))
mod <- lmer(mf_overfitted, inut_data, REML = FALSE)
return(mod)
# }
#
# else{
# return(NULL)
# }
})
output$mod_4 <- renderPrint({
summary(mod_overfitted())
})
# output$mod_11 <- renderUI({
# (HTML(stargazer(mod_1(),mod_2(),mod_3(),mod_4(),type = "html")))
# })
#
# output$mod_22 <- renderUI({
# (HTML(stargazer(mod_3(),mod_4(),type = "html")))
# })
#
##Stargazer's tables
fixed_mod_no_random <- eventReactive(input$go,{
mod_no_random_int <- round(((summary(mod_no_random()))$coefficients[1,1:3]),3)
mod_no_random_slp <- round(((summary(mod_no_random()))$coefficients[2,1:3]),3)
fixed_mod_no_random <- data.frame(mod_no_random_int, mod_no_random_slp)
colnames(fixed_mod_no_random) <- c("Intercept", "Slope")
return(fixed_mod_no_random)
})
fixed_mod_no_slope <- eventReactive(input$go,{
mod_no_slope_int <- round(((summary(mod_no_slope()))$coefficients[1,1:3]),3)
mod_no_slope_slp <- round(((summary(mod_no_slope()))$coefficients[2,1:3]),3)
fixed_mod_no_slope <- data.frame(mod_no_slope_int, mod_no_slope_slp)
colnames(fixed_mod_no_slope) <- c("Intercept", "Slope")
return(fixed_mod_no_slope)
})
random_mod_no_slope <- eventReactive(input$go,{
vc_mod_no_slope <- VarCorr(mod_no_slope())
s_mod_ns <- sqrt(diag(vc_mod_no_slope [["SubjID"]])) # by-subject random intercept, SD
i_mod_ns <- sqrt(diag(vc_mod_no_slope [["ItemID"]])) # by-item random intercept, SD
#attr(vc_mod_no_slope, "sc") # residuals
mod_ns_re <- data.frame(Groups=c("SubjID", "SubjID", "ItemID", "ItemID", "Residuals"),
Name = c("Intercept", "Cond", "Intercept", "Cond","-"),
Std.Dev = c(round((s_mod_ns),3), "-", round((i_mod_ns),3), "-", round((attr(vc_mod_no_slope, "sc")),3) ) )
return(mod_ns_re)
})
fixed_mod <- eventReactive(input$go,{
mod_int <- round(((summary(mod()))$coefficients[1,1:3]),3)
mod_slp <- round(((summary(mod()))$coefficients[2,1:3]),3)
fixed_mod <- data.frame(mod_int, mod_slp)
colnames(fixed_mod) <- c("Intercept", "Slope")
return(fixed_mod)
})
random_mod <- eventReactive(input$go,{
vc_mod <- VarCorr(mod())
s_mod <- sqrt(diag(vc_mod[["SubjID"]])) # by-subject random intercept and slope, SD
i_mod <- sqrt(diag(vc_mod[["ItemID"]])) # by-item random intercept, SD
#attr(vc_mod, "sc") # residuals
#attr(vc_mod[["SubjID"]], "correlation")[1,2] # correlation between by-subject intercept and slope
mod_re <- data.frame(Groups=c("SubjID", "SubjID", "ItemID", "ItemID", "Residuals"),
Name = c("Intercept", "Cond", "Intercept", "Cond","-"),
Std.Dev = c(round((s_mod[1]),3), round((s_mod[2]),3), round((i_mod[1]),3), "-", round ((attr(vc_mod, "sc")),3) ),
Correlation=c( round((attr(vc_mod[["SubjID"]], "correlation")[1,2]), 3), "-", "-", "-", "-"))
return(mod_re)
})
fixed_mod_overfitted <- eventReactive(input$go,{
mod_overfitted_int <- round(((summary(mod_overfitted()))$coefficients[1,1:3]),3)
mod_overfitted_slp <- round(((summary(mod_overfitted()))$coefficients[2,1:3]),3)
fixed_mod_overfitted <- data.frame(mod_overfitted_int, mod_overfitted_slp)
colnames(fixed_mod_overfitted) <- c("Intercept", "Slope")
return(fixed_mod_overfitted)
})
random_mod_overfitted <- eventReactive(input$go,{
vc_mod_overfitted <- VarCorr(mod_overfitted())
s_mod_o <- sqrt(diag(vc_mod_overfitted[["SubjID"]])) # by-subject random intercept and slope, SD
i_mod_o <- sqrt(diag(vc_mod_overfitted[["ItemID"]])) # by-item random intercept and slope, SD
#attr(vc_mod_overfitted, "sc") # residuals
#attr(vc_mod_overfitted[["SubjID"]], "correlation")[1,2] # correlation between by-subject intercept and slope
#attr(vc_mod_overfitted[["ItemID"]], "correlation")[1,2] # correlation between by-item intercept and slope
mod_o_re <- data.frame(Groups=c("SubjID", "SubjID", "ItemID", "ItemID", "Residuals"),
Name = c("Intercept", "Cond", "Intercept", "Cond","-"),
Std.Dev = c(round((s_mod_o[1]),3), round((s_mod_o[2]),3), round((i_mod_o[1]),3), round((i_mod_o[2]),3), round((attr(vc_mod_overfitted, "sc")),3) ),
Correlation=c(round((attr(vc_mod_overfitted[["SubjID"]], "correlation")[1,2]),3), "-", round((attr(vc_mod_overfitted[["ItemID"]], "correlation")[1,2]),3), "-", "-"))
return(mod_o_re)
})
output$fixed_mod_no_random_tab <- renderUI({
(HTML(stargazer(fixed_mod_no_random(),title = "Model with no random efects - fixed effects ",
notes="Model syntax: (Resp ~ Cond)", nobs = FALSE,
font.size = "Huge",
column.sep.width = "10pt",summary = FALSE,
min.max= FALSE, type = "html")))
})
output$fixed_mod_no_slope_tab <- renderUI({
(HTML(stargazer(fixed_mod_no_slope(),title="Model with by-subject random intercept and by-item random intercept - fixed effects",
notes="Model syntax: (Resp ~ Cond + (1| SubjID) + (1 | ItemID))", nobs = FALSE,
font.size = "Huge",
column.sep.width = "10pt",summary = FALSE,
min.max= FALSE, type = "html")))
})
output$random_mod_no_slope_tab <- renderUI({
(HTML(stargazer(random_mod_no_slope(),title="Model with by-subject random intercept and by-item random intercept - random effects",
notes="Model syntax: (Resp ~ Cond + (1| SubjID) + (1 | ItemID))", nobs = FALSE,
font.size = "Huge",
column.sep.width = "10pt",summary = FALSE,
min.max= FALSE, type = "html")))
})
output$fixed_mod_tab <- renderUI({
(HTML(stargazer(fixed_mod(),title="Model with by-subject random intercept + slope and by-item random intercept - fixed effects",
notes="Model syntax: (Resp ~ Cond + (1 + Cond | SubjID) + (1 | ItemID))", nobs = FALSE,
font.size = "Huge",
column.sep.width = "10pt",summary = FALSE,
min.max= FALSE, type = "html")))
})
output$random_mod_tab <- renderUI({
(HTML(stargazer(random_mod(),title="Model with by-subject random intercept + slope and by-item random intercept - random effects",
notes="Model syntax: (Resp ~ Cond + (1 + Cond | SubjID) + (1 | ItemID))", nobs = FALSE,
font.size = "Huge",
column.sep.width = "10pt", summary = FALSE,
min.max= FALSE, type = "html")))
})
output$fixed_mod_overfitted_tab <- renderUI({
(HTML(stargazer(fixed_mod_overfitted(),title="Model with by-item and by-subject intercept and slope - fixed effects",
notes="Model syntax: (Resp ~ Cond + (1 + Cond | SubjID) + (1 + Cond | ItemID))", nobs = FALSE,
font.size = "Huge",
column.sep.width = "10pt",summary = FALSE,
min.max= FALSE, type = "html")))
})
output$random_mod_overfitted_tab <- renderUI({
(HTML(stargazer(random_mod_overfitted(),title="Model with by-item and by-subject intercept and slope - random effects",
notes="Model syntax: (Resp ~ Cond + (1 + Cond | SubjID) + (1 + Cond | ItemID))", nobs = FALSE,
font.size = "Huge",
column.sep.width = "10pt",
min.max= FALSE, summary = FALSE, type = "html")))
})
#Plots
plots_combined <- eventReactive(input$go,{
fitted_overfitted <- fitted(mod_overfitted())
residuals_overfitted <- residuals(mod_overfitted(),"pearson")
fitted_no_slope <- fitted(mod_no_slope())
residuals_no_slope <- residuals(mod_no_slope(),"pearson")
fitted_no_random <- fitted(mod_no_random())
residuals_no_random <- residuals(mod_no_random(),"pearson")
fitted_ <- fitted(mod())
residuals_ <- residuals(mod(),"pearson")
all_fitted <- c(fitted_overfitted, fitted_no_slope, fitted_no_random, fitted_)
all_res <- c(residuals_overfitted, residuals_no_slope, residuals_no_random, residuals_)
names <- c(rep("crossed random effects", length(fitted_overfitted)),
rep("by-item by-subject random intercept", length(fitted_no_slope)),
rep("only fixed effects", length(fitted_no_random)),
rep("by-item by-subject random intercept and by-subject random slope", length(fitted_)))
names <- factor(names, levels = c("crossed random effects", "by-item by-subject random intercept",
"only fixed effects",
"by-item by-subject random intercept and by-subject random slope"), ordered=TRUE)
input <- data.frame(all_fitted, all_res, names)
?geom_smooth
library(ggplot2)
p <- ggplot(data = input, aes(x=all_fitted, y=all_res, colour=names, alpha=0.3)) +
geom_point(size = 2) +
geom_smooth(data = input[input$names != "only fixed effects", ]) +
#geom_jitter() +
geom_abline(col="gold") +
facet_wrap(~ names)+
ggtitle("Fitted vs residual plot")+
theme(plot.title = element_text(lineheight=.7, face="bold", size = 16))+
theme(axis.title.x = element_text(face="bold", size=12))+
theme(axis.title.y = element_text(face="bold", size=12))+
theme(legend.position="none")+
scale_y_continuous(name="Residual values (Pearson residuals)")+
scale_x_continuous(name="Fitted values")
return(p)
})
output$plot <- renderPlot({plots_combined()}, height = 900)
##Model comparison
model_comparison <- eventReactive(input$go,{
aic_data <- AIC(mod_no_random(), mod(), mod_no_slope(), mod_overfitted())
aic_dataframe <- data.frame(Models = c("Fixed Effects Only", "By-S and By-I Random Intercept",
"By-S and By-I Random Intercept, and By-S Random Slope" ,"Crossed Random Effects"),
Df = c(aic_data[,1]),
AIC = c(aic_data[,2]))
return(aic_dataframe)
})
output$AIC <- renderUI({
(HTML(stargazer(model_comparison(),title="Comparing the models",
notes="Model comparison is done by Aikake Information Criteria (AIC)", nobs = FALSE,
font.size = "Huge",
column.sep.width = "10pt",
min.max= FALSE, summary = FALSE, type = "html")))
})
#Select a model tab
observeEvent(input$select_model, {
if (input$select_model == ""){
output$the_answer <- renderText( print("") )
}
if (input$select_model == "fixed"){
output$the_answer <- renderText( print("Try again! This is the simplest model that most people will use, but it does not take into accout individual differences in reaction times
or individual differences in the effect of alcohol in reaction times!") )
}
if (input$select_model == "rand_intercept"){
output$the_answer <- renderText( print("OK, you are almost there, but something is still missing! Now you have included overall individual differences in reaction times
between subjects and between items. But will everyone be affected by alcohol in the same way?") )
}
if (input$select_model == "rand_2" & input$C2 == FALSE){
output$the_answer <- renderText( print("Correct! You used different stimuli for the after alcohol and before alcohol condition, so you only needed a
by-item random intercept! And, as all participants completed both conditions you needed both by-subject random slope and intercept.
Well done!") )
}
if (input$select_model == "rand_2" & input$C2 == TRUE){
output$the_answer <- renderText( print("Hmmm, almost there! All participants completed both conditions, so you correctly included the by-subject random slope and intercept.
And you correctly identified that you need a by-item random intercept. However, you used the same items in both the before and after alcohol conditions,
so there is still something missing!") )
}
if (input$select_model == "rand_3" & input$C2 == FALSE){
output$the_answer <- renderText( print("Oooops, the model you selected has too many random effects given the design you selected! You used different stimuli for the before
and after alcohol conditions, so you only need a by-item random intercept!") )
}
if (input$select_model == "rand_3" & input$C2 == TRUE){
output$the_answer <- renderText( print("Absolutely correct! Well done!") )
}
})
})
### Names
### Reny Baykova, University of Glasgow,
### Anna-Elisabeth Schnell, Katholieke Universiteit Leuven
### Mait Metelitsa, University of Tartu
### Jakob Scheunemann, Carl von Ossietzky Universitaet Oldenburg
### Kadi Tulver, University of Tartu
### Nathasja Franke, University of Leuven
library(shiny)
shinyUI(fluidPage(
includeScript("../../../Matomo-tquant.js"),
titlePanel("Comparing model specifications"),
sidebarLayout(
sidebarPanel("Data simulation parameters",
#uiOutput("model_types_"),
#checkboxInput("C1", label="Select all the model specification types", value = FALSE),
sliderInput("nsubj", "Set the number of subjects", min = 10, max = 100, value = 50, step = 2),
sliderInput("nitem", "Set the number of items", min = 10, max = 100, value = 50, step = 2),
sliderInput("int", "Set the fixed intercept", min = 1, max = 500, value = 250),
sliderInput("eff", "Set the fixed effect of condition", min = 1, max = 100, value = 50),
sliderInput("err", "Set the residual SD", min = 1, max = 100, value = 50),
sliderInput("t00", "Set the by-subject random intercept (SD)", min = 1, max = 100, value = 50),
sliderInput("t11", "Set the by-subject random slope (SD)", min = 1, max = 40, value = 20),
sliderInput("rsub", "Set the by-subject random correlation (-1 to 1)", min = -1, max = 1, value = 0,step = 0.1),
sliderInput("w00", "Set the by-item random intercept", min = 1, max = 20, value = 10),
sliderInput("w11", "Set the by-item random slope", min = 1, max = 20, value = 10),
sliderInput("ritm", "Set the by-item random correlation", min = -1, max = 1, value = 0,step = 0.1),
checkboxInput("C2", label="Check this box to use within-item condition (defaults to between-item)", value = FALSE),
actionButton("go","See the results"),
downloadButton('downloadData','Save my file!')
),
mainPanel(
tags$style(type="text/css",
".shiny-output-error { visibility: hidden; }",
".shiny-output-error:before { visibility: hidden; }"),
tabsetPanel(
tabPanel(title = "Description of our app",
textOutput("text"),
fluidRow(includeMarkdown("documentation.md"))),
tabPanel("Explore the data",
DT::dataTableOutput('tbl')),
tabPanel("Comparison of simulation parameters and models",
#uiOutput(outputId = "mod_22"),
#uiOutput(outputId = "mod_11"),
#verbatimTextOutput("mod_1"),
#verbatimTextOutput("mod_2"),
#verbatimTextOutput("mod_3"),
#verbatimTextOutput("mod_4"),
uiOutput(outputId = "fixed_mod_no_random_tab"),
p(""),
p(""),
uiOutput(outputId = "fixed_mod_no_slope_tab"),
p(""),
p(""),
uiOutput(outputId = "random_mod_no_slope_tab"),
p(""),
p(""),
uiOutput(outputId = "fixed_mod_tab"),
p(""),
p(""),
uiOutput(outputId = "random_mod_tab"),
p(""),
p(""),
uiOutput(outputId = "fixed_mod_overfitted_tab"),
p(""),
p(""),
uiOutput(outputId = "random_mod_overfitted_tab"),
p(""),
p(""),
uiOutput(outputId = "AIC")
),
tabPanel("Diagnostic plots for fitted models",
plotOutput("plot")),
tabPanel("Select a model",
sidebarLayout(
sidebarPanel(
radioButtons("select_model", "Select a model:", c("None selected" = "",
"Fixed effects only" = "fixed",
"By-subject and by-item random intercept" = "rand_intercept",
"By-subject and by-item random intercept, and by-subject random slope"= "rand_2",
"By-subject and by-item random intercept and slope (crossed random effects)" = "rand_3"))
),
mainPanel(
textOutput(outputId = "the_answer")
) ))
)))))