shinyApp(
ui = tagList(
navbarPage(
theme = "united", # <--- To use a theme, uncomment this
"Probabilistic Knowledge Structures",
tabPanel("Estimation",
includeScript("../../../Matomo-tquant.js"),
sidebarPanel(
tags$h4("Model specification"), br(),
tags$h5("Change initial values:"),
checkboxInput("initCareless", "Careless errors", value = FALSE),
checkboxInput("initLuckyguess", "Lucky guesses", value = FALSE),
checkboxInput("initResponses", "Response frequencies", value = FALSE),
selectInput("analysis", "Analysis types:",
choices = c("Bayesian analysis" = "bayesian",
"Classical analysis" = "classic"),
selected = "bayesian"
)
),
mainPanel(
tabsetPanel(
tabPanel("Analysis",
h4("Initialize your model:"),
uiOutput("initializeCareless"),
uiOutput("initializeLuckyguess"),
uiOutput("initializeResponses"),
uiOutput("Blim")
),
tabPanel("Additional information",
h4("Convergence diagnostic for Bayesian parameter estimation"),
p("This panel provides information about the quality of the obtained MCMC samples,
by means of Gelman and Rubins (1992) scale reduction factor and traceplots for each
parameter of interest. Note that a scale reduction factor close to
1 indicates good convergence. In addition, a good overlap between the chains in the traceplot
indicates that the chains have converged from their initial starting values to
their stationary distribution. If convergence cannot be assumed, try out a larger number
of samples, a longer burn-in period or more thinning. Details about Bayesian cognitive modelling
and MCMC methods can be found, e.g., in Lee and Wagenmakers (2013)."),
dataTableOutput("convergenceTable"),
plotOutput("bayesTrace")
)
)
)
)
)
),
server = function(input, output) {
## source all necessary scripts
source("functions.R")
source("PlotFunctions.R")
# initialize parameter values
output$initializeResponses <- renderUI({
if(input$initResponses){
fluidPage(
h5("Insert the number of observed response patterns:"),
fluidRow(
column(3,
numericInput("NR1", "0", 0, 1e5, value = 10)),
column(3,
numericInput("NR2", "bc", 0, 1e5, value = 10)),
column(3,
numericInput("NR3", "bd", 0, 1e5, value = 10)),
column(3,
numericInput("NR4", "abc", 0, 1e5, value = 10))
),
fluidRow(
column(3,
numericInput("NR5", "acd", 0, 1e5, value = 10)),
column(3,
numericInput("NR6", "abd", 0, 1e5, value = 10)),
column(3,
numericInput("NR7", "abcd", 0, 1e5, value = 10))
)
)
}
})
output$initializeCareless <- renderUI({
if(input$initCareless){
fluidPage(
h5("Define the probability that a person makes a careless error on the following items:"),
fluidRow(
column(3,
numericInput("beta1", "Item a", 0, 1, value = 0.1, step = 0.1)),
column(3,
numericInput("beta2", "Item b", 0, 1, value = 0.1, step = 0.1)),
column(3,
numericInput("beta3", "Item c", 0, 1, value = 0.1, step = 0.1)),
column(3,
numericInput("beta4", "Item d", 0, 1, value = 0.1,step = 0.1))
)
)
}
})
output$initializeLuckyguess <- renderUI({
if(input$initLuckyguess){
fluidPage(
h5("Define the probability that a person makes a lucky guess on the following items:"),
fluidRow(
column(3,
numericInput("eta1", "Item a", 0, 1, value = 0.1, step = 0.1)),
column(3,
numericInput("eta2", "Item b", 0, 1, value = 0.1, step = 0.1)),
column(3,
numericInput("eta3", "Item c", 0, 1, value = 0.1, step = 0.1)),
column(3,
numericInput("eta4", "Item d", 0, 1, value = 0.1, step = 0.1))
)
)
}
})
# Frequentist Analysis
goFrequentistAnalysis <- eventReactive(input$runFrequentistAnalysis, {
if(input$initLuckyguess){
etaInput <- c(input$eta1, input$eta2, input$eta3, input$eta4)
} else {
etaInput <- rep(0.1, 4)
}
if(input$initCareless){
betaInput <- c(input$beta1, input$beta2, input$beta3, input$beta4)
} else {
betaInput <- rep(0.1, 4)
}
if(input$initResponses){
NRInput <- c(input$NR1, input$NR2, input$NR3, input$NR4, input$NR5, input$NR6, input$NR7)
} else {
NRInput <- rep(10, 7)
}
list(method = input$estimationMethod,
eta = etaInput,
beta = betaInput,
NR = NRInput)
})
output$freqBlimPlot <- renderPlot({
print(goFrequentistAnalysis())
blim1 <- frequentistParameterEstimation(N.R = goFrequentistAnalysis()$NR,
estimation.method = goFrequentistAnalysis()$method,
beta.parameters = goFrequentistAnalysis()$beta,
eta.parameters = goFrequentistAnalysis()$eta)
print(blim1)
blimPlot(blim1$beta, blim1$eta, blim1$P.K)
})
# Bayesian Analysis
goBayesianAnalysis <- eventReactive(input$runBayesianAnalysis, {
if(input$initLuckyguess){
eta.percentages <- c(input$eta1, input$eta2, input$eta3, input$eta4)
} else {
eta.percentages <- rep(0.5, 4)
}
if(input$initCareless){
beta.percentages <- c(input$beta1, input$beta2, input$beta3, input$beta4)
} else {
beta.percentages <- rep(0.5, 4)
}
if(input$initResponses){
NR.Input <- c(input$NR1, input$NR2, input$NR3, input$NR4, input$NR5, input$NR6, input$NR7)
} else {
NR.Input <- rep(10, 7)
}
list(out = bayesianParameterEstimation(N.R = NR.Input,
beta.percentages = beta.percentages,
eta.percentages = eta.percentages,
n.iter = input$niter,
n.thin = input$nthin,
n.burnin = input$nburnin),
item.number = as.numeric(input$plotPosterior))
})
output$bayesSummaryPlotPK <- renderPlot({
plotBayesSummary(goBayesianAnalysis()$out$PK$info, param = "PK")}, height= 450, width = 600)
output$bayesSummaryPlotBeta <- renderPlot({
plotBayesSummary(goBayesianAnalysis()$out$beta$info, param = "beta")}, height= 250, width = 600)
output$bayesSummaryPlotEta <- renderPlot({
plotBayesSummary(goBayesianAnalysis()$out$eta$info, param = "eta")}, height = 250, width = 600)
output$bayesPriorPosteriorPlot <- renderPlot({
item.number <- goBayesianAnalysis()$item.number
print(item.number)
bayesHist(goBayesianAnalysis()$out, item.number = item.number)
}, height = 1000, width = 600 )
output$bayesBlimPlot <- renderPlot({
blimPlot(goBayesianAnalysis()$out$beta$info$median,
goBayesianAnalysis()$out$eta$info$median,
goBayesianAnalysis()$out$PK$info$median)
})
# Bayesian convergence diagnostics
output$convergenceTable <- renderDataTable({
convergenceDiagnostics(goBayesianAnalysis()$out$samples)
})
output$bayesTrace <- renderPlot({
bayesTrace(goBayesianAnalysis()$out$samples)
}, height = 1000, width = 600)
# Output
output$Blim <- renderUI({
if(input$analysis == "classic"){
fluidPage(
selectInput("estimationMethod", label = "Choose a parameter estimation method:",
choices = c("Maximum likelihood method (ML)" = "ML",
"Minimum discrepancy method (MD)" = "MD",
"Mixed (MDML)" = "MDML"),
selected = "ML"),
actionButton("runFrequentistAnalysis", "Estimate model", class = "btn-primary"),
plotOutput("freqBlimPlot")
)
} else {
fluidPage(
numericInput("niter", "Specify number of iterations", 15e4, min = 6000, max = 1e6),
numericInput("nburnin", "Specify the burn-in period", 5000, min = 0, max = 15000),
numericInput("nthin", "Specify the thinning", 10, min = 0, max = 200),
checkboxGroupInput("plotPosterior", "Plot the prior and posterior distribution for the following items:",
c("Item a" = 1,
"Item b" = 2,
"Item c" = 3,
"Item d" = 4),
selected = c(1:4),
inline = TRUE),
actionButton("runBayesianAnalysis", "Estimate model", class = "btn-primary"),
plotOutput("bayesBlimPlot"),
plotOutput("bayesSummaryPlotPK", height = "600px"),
plotOutput("bayesSummaryPlotBeta"), plotOutput("bayesSummaryPlotEta"),
plotOutput("bayesPriorPosteriorPlot")
)
}
})
# end of the app
}
)
# load packages
library(pks)
library(R2jags)
library(ggplot2)
# general settings
# N.R <- rep(10, 7)
# names(N.R) <- as.pattern(K)
# Frequentist parameter estimation
frequentistParameterEstimation <- function(N.R, estimation.method, beta.parameters = rep(0.1, 4),
eta.parameters = rep(0.1, 4)){
# estimates the frequentist model for blim
K <- as.binmat(c("0000", "0110", "0101", "1110", "1011", "1101", "1111"))
rownames(K) <- as.pattern(K)
names(N.R) <- as.pattern(K)
classic.blim <- blim(K, N.R, method = estimation.method, beta = beta.parameters,
eta = eta.parameters)
}
# Bayesian parameter estimation
bayesianParameterEstimation <- function(N.R, beta.percentages, eta.percentages, n.iter, n.thin, n.burnin){
# estimates the bayesian graphical model for blim
# general settings
K <- as.binmat(c("0000", "0110", "0101", "1110", "1011", "1101", "1111"))
rownames(K) <- as.pattern(K)
nitems <- ncol(K)
nstates <- nrow(K)
names(N.R) <- as.pattern(K)
R <- as.binmat(N.R)
npatterns <- 7
n.observ <- sum(N.R)
# beta and eta priors
beta.a <- beta.percentages*10
beta.b <- 10- beta.a
eta.a <- eta.percentages*10
eta.b <- 10- eta.a
# set up Bayesian graphical model
data <- list("N.R", "n.observ", "R", "K", "nitems", "nstates", "npatterns",
"beta.a", "beta.b", "eta.a", "eta.b")
myinits <- list(list(beta = rep(0.1, nitems), eta = rep(0.1, nitems),
P.K = rep(1/nstates, nstates)),
list(beta = rep(0.1, nitems), eta = rep(0.1, nitems),
P.K = rep(1/nstates, nstates)))
parameters <- c("P.K", "eta", "beta")
# all posterior samples
samples <- jags(data, inits = myinits, parameters,
model.file = "blim_model.txt", n.chains = 2, n.iter = n.iter,
n.thin = n.thin, n.burnin=n.burnin)
# information about beta parameters
beta <- list(post.samples = samples$BUGSoutput$sims.list[["beta"]],
prior.samples = sapply(1:4, function(i) rbeta(1e5, beta.a[i], 10-beta.a[i])),
info = estimateMedian(samples$BUGSoutput$sims.list[["beta"]]))
# information about eta parameters
eta <- list(post.samples = samples$BUGSoutput$sims.list[["eta"]],
prior.samples = sapply(1:4, function(i) rbeta(1e5, eta.a[i], 10-eta.a[i])),
info = estimateMedian(samples$BUGSoutput$sims.list[["eta"]]))
PK <- list(post.samples = samples$BUGSoutput$sims.list[["P.K"]],
prior.samples = sapply(1:7, function(i) rbeta(1e5, 1, 6)),
info = estimateMedian(samples$BUGSoutput$sims.list[["P.K"]]))
# output
output <- list(samples = samples,
beta = beta,
eta = eta,
PK = PK)
return(output)
}
# convergence diagnostics
convergenceDiagnostics <- function(samples, item.number = c(1:4)){
# calculates for each item Gelmans R for beta and eta parameter
convergenceTable <- data.frame(Item = NA, Scale.reduction.factor = NA, Upper.CI = NA)
for(i in seq_along(item.number)){
item.name.beta <- paste0("beta[", i, "]")
item.name.eta <- paste0("eta[", i, "]")
# transform samples into mcmc list
beta.list <- list(mcmc(samples$BUGSoutput$sims.array[, 1, item.name.beta]), # chain 1
mcmc(samples$BUGSoutput$sims.array[, 2, item.name.beta])) # chain 2
eta.list <- list(mcmc(samples$BUGSoutput$sims.array[, 1, item.name.beta]), # chain 1
mcmc(samples$BUGSoutput$sims.array[, 2, item.name.beta])) # chain 2
# compute gelmans R
convergenceTable[i, "Item"] <- i
convergenceTable[i, "Scale.reduction.factor"] <- round(gelman.diag(beta.list)$psrf[1],5)
convergenceTable[i, "Upper.CI"] <- round(gelman.diag(beta.list)$psrf[2],5)
}
return(convergenceTable)
}
# Calculate median plus HDI
estimateMedian <- function(samples) {
medianPlusHDI <- data.frame(lower = NA, median = NA, upper = NA)
for(i in 1:ncol(samples)){
d <- quantile(samples[,i], c(0.025, 0.5, 0.975))
medianPlusHDI[i,] <- d
}
return(medianPlusHDI)
}
## Customize beta and eta priors
# # adjust priors for beta and eta
# beta.percentages <- c(50, 60, 10, 90)
# eta.percentages <- c(50, 60, 20, 80)
# # beta and eta priors
# beta.a <- beta.percentages/10
# beta.b <- 10- beta.a
# eta.a <- eta.percentages/10
# eta.b <- 10- eta.a
# n.iter <- 120000
# n.thin <- 10
# n.burnin <- 5000
# frequentist Plot functions
NRplot <- function(N.R){
# N.R: named vector with N.R values
sets <- c("0", "a", "b", "c", "d", "ab", "ac", "ad", "bc",
"bd", "cd", "abc", "abd", "acd", "bcd", "abcd")
N.R <- N.R[order(match(names(N.R), sets))] # match order of sets and P.K
N.R <- as.table(N.R)
barplot(N.R,
horiz = TRUE,
names.arg = names(N.R),
axes = FALSE,
cex.names = 1.3,
las = 2,
col = "grey36",
xlim = c(0, max(N.R) + 5))
# include axis
axis(side = 3,
cex.axis = 1.3)
}
betaPlot <- function(beta.vector){
# beta.vector: vector with beta values (beta.a, beta.b, beta.c, beta.d)
beta <- as.table(rev(beta.vector)) # Reverse Eta
barplot(beta,
horiz = TRUE,
names.arg = expression(beta[d], beta[c], beta[b], beta[a]),
axes = FALSE,
cex.names = 2,
las = 2,
col = "grey36")
# include axis
axis(side = 3,
cex.axis = 1.5,
at = c(0, 0.1, 0.2),
labels = c("0", ".1", ".2"))
}
etaPlot <- function(eta.vector){
# eta.vector: vector with eta values (eta.a, eta.b, eta.c, eta.d)
eta <- as.table(rev(eta.vector)) # Reverse Eta
barplot(eta,
horiz = TRUE,
names.arg = expression(eta[d], eta[c], eta[b], eta[a]),
axes = FALSE,
cex.names = 2,
las = 2,
col = "grey36")
axis(side = 3, cex.axis = 1.5,
at = c(0, 0.1, 0.2),
labels = c("0", ".1", ".2"))
}
ksPlot <- function(P.K){
# P.K: vector with P.K values
hasNoNames <- is.null(names(P.K))
if(hasNoNames){
names(P.K) <- c("0000", "0110", "0101", "1110", "1011", "1101", "1111")
}
sets <- c("bc", "abc", "acd", "abcd", "bd", "abd")
P.K <- P.K[order(match(names(P.K), sets))] # match order of sets and P.K
x0 <- c(4, 1, 7, 1, 4, 7, 4) # coordinates
y0 <- c(1, 4, 4, 7, 7, 7, 10)-0.5
# Plot Knowledge Structure
# Empty Plot
plot(1, type = "n", xlab = "", ylab = "",
xlim = c(0, 10),
ylim = c(0, 12),
axes = FALSE)
# Add Lines
polygon(x = c(5, 2, 2, 5, 5, 5, 8, 8, 5),
y = c(1, 4, 7, 10, 7, 1, 4, 7, 10),
border = "grey36",
lwd = 2)
# Add Points
symbols(x = c(2, 2, 5, 5, 5, 8, 8),
y = c(4, 7, 1, 7, 10, 4, 7),
circles = rep(0.6, 7),
add = TRUE,
inches = FALSE,
bg = "lightgrey",
fg = "grey36",
lwd = 2)
# Add Knowledge States
text(x = c(2, 2, 5, 5, 8, 8),
y = c(4, 7, 7, 10, 4, 7),
labels = sets,
col = "darkred",
cex = 1.5)
text(x = 5,
y = 1,
labels = expression(symbol("\306")), # Empty Set
col = "darkred",
cex = 1.5)
# Barplots for P.K values
for(i in 1:7){
rect(x0[i] + 0.1, y0[i], x0[i] + 0.3, y0[i] + (P.K[i] * 4),
col = "red3",
border = "grey36")
}
# Add Axis of Barplots
segments(x0 = x0, y0 = y0, x1 = x0, y1 = y0 + 0.8) # Axis
segments(x0 = x0, y0 = y0, x1 = x0 - 0.1, y1 = y0) # Ticks
segments(x0 = x0, y0 = y0 + 0.8, x1 = x0 - 0.1, y1 = y0 + 0.8)
segments(x0 = x0, y0 = y0 + 0.4, x1 = x0 - 0.1, y1 = y0 + 0.4)
text(x0, y0, pos = 2, labels = rep("0", 8)) # Labels
text(x0, y0 + 0.8, pos = 2, labels = rep(".2", 8))
}
blimPlot <- function(beta.vector, eta.vector, P.K){
# define layout
matrix <- matrix(c(1, 1, 3, 3, 3, 3, 3, 3, 3,
1, 1, 3, 3, 3, 3, 3, 3, 3,
1, 1, 3, 3, 3, 3, 3, 3, 3,
2, 2, 3, 3, 3, 3, 3, 3, 3,
2, 2, 3, 3, 3, 3, 3, 3, 3,
2, 2, 3, 3, 3, 3, 3, 3, 3), 6, 9, byrow = TRUE)
layout(matrix)
betaPlot(beta.vector)
etaPlot(eta.vector)
ksPlot(P.K)
}
# Bayesian Plot functions
bayesTrace <- function(samples){
# plot the traceplots for beta and eta for all items
par(mfrow=c(6, 2),
cex.axis = 1.3, las = 1, bty="n", mgp = c(3.5, 1, 0),
cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5,
font.lab = 2, cex.axis = 1.3, bty = "n")
for(i in 1:4){
item.name.beta <- paste0("beta[", i, "]")
item.name.eta <- paste0("eta[", i, "]")
# samples for beta and eta parameter
beta.vector <- samples$BUGSoutput$sims.array[, , item.name.beta]
eta.vector <- samples$BUGSoutput$sims.array[, , item.name.eta]
plot(beta.vector[, 1], type = 'l', col = "lightcoral", bty ="n", ylim = c(0,1),
ylab = eval(bquote(expression(beta[.(letters[i])]))))
lines(beta.vector[, 2], col = "lightblue3")
plot(eta.vector[, 1], type = 'l', col = "lightcoral", bty = "n", ylim = c(0,1),
ylab = eval(bquote(expression(eta[.(letters[i])]))))
lines(eta.vector[, 2], col = "lightblue3")
}
}
bayesHist <- function(samples, item.number){
par(mfrow = c(4, 2))
for(i in item.number){
plotPriorPosterior(samples$beta$prior.samples, samples$beta$post.samples,
info = samples$beta$info, param = "beta", item.number = i)
plotPriorPosterior(samples$eta$prior.samples, samples$eta$post.samples,
info = samples$eta$info, param = "eta", item.number = i)
}
}
plotPriorPosterior <- function(prior.samples, post.samples, info, param = c("beta", "eta"), item.number = 1:4){
par(cex.axis = 1.3, las = 1, bty="n", mgp = c(3.5, 1, 0),
cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5,
font.lab = 2, cex.axis = 1.3, bty = "n")
for(i in item.number){
# parameter
y.high <- round(max(max(density(prior.samples[,i])$y), max(density(post.samples[,i])$y)))
y.high.10 <- y.high + ceiling(y.high/100*10) # plus 10%
y.high.20 <- y.high + ceiling(y.high/100*20) # plus 20%
x.legend.pos <- ifelse(median(density(post.samples[,i])$x) > 0.5, 0, 1)
x.lab <- ifelse(param == "beta", eval(bquote(expression(beta[.(letters[i])]))), eval(bquote(expression(eta[.(letters[i])]))))
# plot
plot(density(post.samples[,i], adjust = 2),
xlim = c(0, 1), ylim = c(0, y.high.20),
xlab = x.lab, main = "", lwd = 2)
par(new = TRUE)
plot(density(prior.samples[,i], adjust = 2),
axes = FALSE, xlim = c(0, 1),
xlab = "", ylab = "", main = "",
ylim = c(0, y.high.20), lwd = 2, lty = 3)
# add legend
legend(x.legend.pos, y.high.20, legend = c("Posterior", "Prior"),
lty = c(1, 3), bty = "n", lwd = c(2, 2), cex = 1.2,
xjust = x.legend.pos, yjust = 1, x.intersp = 0.6, seg.len = 1.2)
# add parameter information
mtext(paste("median =", round(info$median,2)[i], sep = " "), side = 3, line = 2.5, adj = 1)
mtext(paste("95% CI: [", round(info$lower,2)[i], ", ", round(info$upper,2)[i], "]", sep = ""), side = 3, line = 1, adj = 1)
arrows(info$lower[i], y.high.10, info$upper[i], y.high.10, code = 3, length = 0.05, angle = 90)
}
}
plotBayesSummary <- function(info, param = c("beta", "eta", "PK")){
# general settings
K <- as.binmat(c("0000", "0110", "0101", "1110", "1011", "1101", "1111"))
# create axis labels
x <- NULL
if(param == "PK"){
x <- as.pattern(K, as.letters=TRUE)
} else {
for(i in 1:4){
x[i] <- ifelse(param == "beta", eval(bquote(expression(beta[.(letters[i])]))), eval(bquote(expression(eta[.(letters[i])]))))
}
x <- unlist(x)
}
# add label and text
info$label <- 1:nrow(info)
info$text <- paste(sprintf("%.2f", round(info$median, 2)), " [",
sprintf("%.2f", round(info$lower, 2)), ", ",
sprintf("%.2f", round(info$upper, 2)), "]", sep = "")
# Create summary plot
p <- ggplot(info, aes(x = info$label, y = median, ymin = lower, ymax = upper)) +
scale_y_continuous("Posterior Parameter Estimates", breaks = seq(0, 1, by = 0.2), limits=c(0,1)) +
coord_flip() +
scale_x_reverse(name = "", breaks = 1:nrow(info), labels = x, minor_breaks = waiver()) +
geom_pointrange() +
geom_point(shape = 21, color = "black", fill = "grey", size = 4.5) +
theme_classic(base_size = 20) +
# extend grid
theme(plot.margin = grid::unit(c(2, 10.5, 1, 1), "lines")) +
# annotations on right side
geom_text(aes(y = Inf, x = info$label, label = info$text), size = 5, hjust = -0.1) +
annotate("text", x = -Inf, y=Inf, label = "median 95% CI", hjust = -0.1, vjust = -1, size = 5)
# Override clipping
gt <- ggplot_gtable(ggplot_build(p))
gt$layout$clip[gt$layout$name == "panel"] <- "off"
grid::grid.draw(gt)
}