What are we going to do?


This app refers to the presentation 'Probabilistic Knowledge Structures' by Heller and Wickelmaier. It will demonstrate parameter estimation in probabilistic knowledge structures:

  • If you want to be reminded of the three estimation procedures, proceed to the tab 'Estimation methods' above (For a more comprehensive overview, see the presentation by Heller and Wickelmaier)
  • If you are already confident in your knowledge about parameter estimation in probabilistic knowledge structures, you can directly proceed to 'Application' on the left and see how changes to the observed data and the assumed knowledge structure will affect parameter estimation
  • In the end, you can test your knowledge with a small quiz

The three estimation methods


With a given knowledge structure and observed response patterns, there are three methods to estimate the parameters of a basic local independence model (BLIM):

  • Maximum Likelihood (ML) Estimation
  • Minimum Discrepancy (MD) Estimation
  • Minimum Discrepancy ML Estimation (MDML)

Method Principle Pros Cons
ML estimates parameters that maximize the probability of the observed data
  • driven by the likelihood of the data
  • (approximately) unbiased estimates
  • iterative (EM algorithm)
  • might inflate error rates for good fit
MD assumes that any response pattern is generated by the knowledge state closest to it
  • computationally efficient (explicit estimators)
  • avoids inflating the error rates
  • ignores the likelihood of the data
  • estimates not unbiased
MDML ML estimation under certain MD restrictions
  • minimizes the expected number of response errors
  • maximizes the likelihood under this constraint
  • reference for quantifying the amount of fit obtained by inflating error rates
  • estimates not unbiased

For this example, we will be working with a knowledge domain of four items: Q = {a,b,c,d}


1) Set the observed response frequencies of ...



2) Set your knowledge structure 𝒦

You have selected the following knowledge structure:

3) Calculate your parameter estimates

The error rate estimates for each item are:


The estimates of the state probabilities are:

4) Calculated expected frequencies with each parameter set


5) Simulate response patterns with each parameter set

To which of the three parameter estimation methods do these statements apply?

Your results

show with app
# app.R
# Visualisation of Parameter Estimation
#
# Enter response frequencies and visualize parameter estimates
#
# Claudia Glemser, last edited: 26/Jan/17

library(shiny)
library(shinydashboard)
library(pks)
library(xtable)

############### UI ###############
ui <- dashboardPage(
  dashboardHeader(
    title = "Parameter estimation in probabilistic knowledge structures",
    titleWidth = "100%"
  ),

  dashboardSidebar(
    sidebarMenu(
      menuItem("Introduction", tabName = "theory",
	       icon = icon("info-circle")),
      menuItem("Application", tabName = "turn", badgeLabel = "Try it!",
	       badgeColor = "blue", icon = icon("pencil-square-o")),
      menuItem("Quiz", tabName = "quiz", badgeLabel = "Got it?",
	       badgeColor = "blue", icon = icon("check-square-o"))  # question-circle-o
    )
  ),
  
  dashboardBody(
    includeScript("../../../Matomo-tquant.js"),
    tabItems(
      tabItem(tabName = "theory",
        tabsetPanel(
          tabPanel("Intro",
	    HTML(
              "<h3 align = 'center'><b> What are we going to do? </h3></b><br>
	       <div> This app refers to the presentation 'Probabilistic
	       Knowledge Structures' by Heller and Wickelmaier. It will
	       demonstrate parameter estimation in probabilistic knowledge
	       structures:</div><br><ul>
	         <li> If you want to be reminded of the three estimation
	              procedures, proceed to the tab 'Estimation methods'
	              above (For a more comprehensive overview, see the
		      presentation by Heller and Wickelmaier)</li>
	         <li> If you are already confident in your knowledge about
		      parameter estimation in probabilistic knowledge
		      structures, you can directly proceed to 'Application' on
		      the left and see how changes to the observed data and
		      the assumed knowledge structure will affect parameter
		      estimation </li>
		 <li> In the end, you can test your knowledge with a small quiz
		 </li>
	         </ul>"
	     )
          ),


	  tabPanel("Estimation methods",
	    HTML(  
	      "<h3 align ='center'><b> The three estimation methods</b></h3>
	       <br>
	       With a given knowledge structure and observed response
	       patterns, there are three methods to estimate the parameters
	       of a basic local independence model (BLIM):
	       <br><br><ul>
	         <li> Maximum Likelihood (ML) Estimation</li>
	         <li> Minimum Discrepancy (MD) Estimation</li>
		 <li> Minimum Discrepancy ML Estimation (MDML)</li>
	       </ul><br></h5>
	       
	       <table style = 'width:100%' id = 't01'>
	       <tr>
	         <th>Method</th>
		 <th>Principle</th>
		 <th>Pros</th>
		 <th>Cons</th>
	       </tr>
	       <tr bgcolor = '#ddd'>
		 <td>ML</td>
		 <td>estimates parameters that maximize the probability of the
		      observed data</td>
		  <td><ul>
		        <li> driven by the likelihood of the data</li>
		        <li> (approximately) unbiased estimates </li>
		      </ul>
		  </td>
		  <td><ul>
		        <li> iterative (EM algorithm) </li>
		        <li> might inflate error rates for good fit </li>
		      </ul>
		  </td>
		</tr>
		<tr bgcolor = '#fff'>
		  <td>MD</td>
		  <td>assumes that any response pattern is generated by the
		  knowledge state closest to it</td>
		  <td><ul>
		    <li> computationally efficient (explicit estimators) </li>
		    <li> avoids inflating the error rates </li>
		  </ul>
		  <td><ul>
		    <li> ignores the likelihood of the data </li>
		    <li> estimates not unbiased </li>
		  </ul></td>
		</tr>
		<tr bgcolor = '#ddd'>
		  <td>MDML</td>
		  <td>ML estimation under certain MD restrictions</td>
		  <td><ul>
		    <li>minimizes the expected number of response errors</li>
		    <li>maximizes the likelihood under this constraint</li>
		    <li>reference for quantifying the amount of fit obtained
		        by inflating error rates</li>
		  </ul></td>
		  <td><ul> 
		    <li> estimates not unbiased </li>
		  </ul></td>
		</tr>
	      </table>"),
	      
            tags$head(tags$style(HTML(
	      "table, th, td {
                border: 1px solid black;
                border-collapse: collapse;
              }

              th {
	        text-align: center;
	      }
	      
	      th, td {
	        padding: 10px;
	      }
	      
	      table#t01 th {
	        background-color: #367fa9;
	        color: white;
              }")))	      
          )
        )
      ),

      tabItem(tabName = "turn",
        tabsetPanel(
          tabPanel("Response frequencies",
            br(), HTML(
	      "<h4 align='center'> For this example, we will be working with a
	      knowledge domain of four items: Q = {a,b,c,d}</h4>"),
	    br(),
	    h4(strong("1) Set the observed response frequencies of ...")),
	    br(), 
            fluidRow(
              column(1,
                numericInput("NR1", HTML("&#8709"), value = 15,
                  min = 0, max = 100, step = 1),
                numericInput("NR2", "{a}", value = 4,
                  min = 0, max = 100, step = 1)),

	      column(1,
                numericInput("NR3", "{b}", value = 6, 
                  min = 0, max = 100, step = 1),
                numericInput("NR4", "{c}", value = 4,
                  min = 0, max = 100, step = 1)),

              column(1,
	        numericInput("NR5", "{d}", value = 5,
                  min = 0, max = 100, step = 1),
                numericInput("NR6", "{a,b}", value = 2,
                  min = 0, max = 100, step = 1)),

	      column(1,
                numericInput("NR7", "{a,c}", value = 3,
                  min = 0, max = 100, step = 1),
                numericInput("NR8", "{a,d}", value = 4,
                  min = 0, max = 100, step = 1)),

              column(1,
                numericInput("NR9", "{b,c}", value = 18,
                  min = 0, max = 100, step = 1),
                numericInput("NR10", "{b,d}", value = 22,
                  min = 0, max = 100, step = 1)),
	
              column(1,
                numericInput("NR11", "{c,d}", value = 2,
                  min = 0, max = 100, step = 1),
                numericInput("NR12", "{a,b,c}", value = 39,
                  min = 0, max = 100, step = 1)),
		    
	      column(1,
                numericInput("NR13", "{a,b,d}", value = 37,
                  min = 0, max = 100, step = 1),
                numericInput("NR14", "{a,c,d}", value = 12,
                  min = 0, max = 100, step = 1)),

	      column(1,
                numericInput("NR15", "{b,c,d}", value = 7,
                  min = 0, max = 100, step = 1),
                numericInput("NR16", "{a,b,c,d}", value = 20,
                  min = 0, max = 100, step = 1)),

	      column(4,
		br(),
		actionButton("go.plot", "Submit frequencies N.R")
	      )
            ),
            
            fluidRow(
	      column(12,
	        plotOutput("plot.NR")
              )
            )	      
          ),

	  tabPanel("Knowledge structure",   
	    fluidRow(
              column(4,
	        HTML("<b><h4>2) Set your knowledge structure</b> &#119974
		     </h4>"),
                checkboxGroupInput("Ks", NULL,
                c("0" = "0000", "{a}" = "1000", "{b}" = "0100",
                  "{c}" = "0010", "{d}" = "0001",
                  "{a,b}" = "1100", "{a,c}" = "1010",
                  "{a,d}" = "1001", "{b,c}" = "0110",
                  "{b,d}" = "0101", "{c,d}" = "0011",
                  "{a,b,c}" = "1110", "{a,b,d}" = "1101",
                  "{a,c,d}" = "1011", "{b,c,d}" = "0111",
                  "{a,b,c,d}" = "1111"),
	        selected = c("0000", "0110", "0101", "1110", "1101",
	          "1011", "1111")
                ),
	      actionButton("go.K", "Submit knowledge structure")),
	
	      column(8,
		wellPanel(
	          h4("You have selected the following knowledge structure:",
		     br(), textOutput("text.str"), br(), uiOutput("text.corr")
	            )
	         )  
	      )    
	    )
          ),

          tabPanel("Parameter estimates",
            
            h4(strong("3) Calculate your parameter estimates")),
	    actionButton("go.para", "Calculate estimates"),

	    fluidRow(
	      column(8,
                tableOutput("tab"),

	        tags$head(
                  tags$style(HTML("
                    .shiny-output-error-validation {
                    color: red;
		    text-align: center;
		    text-transform: uppercase;
   		    }
                 ")))
	      ),

	      column(4,
		uiOutput("text.iter"),
		uiOutput("text.iden")
              )
	    ),

	    HTML("<h4><b>The error rate estimates for each item are:</b></h4>"
		 ),

	    fluidRow(
	      column(6,
                plotOutput("plot.eta")
              ),

	      column(6,
	        plotOutput("plot.beta")
              )
	    ),

	    br(),
	    HTML("<h4><b>The estimates of the state probabilities are:</b></h4>"
		 ),

	    plotOutput("plot.pi")
	  ),
	
	  tabPanel("Expected and simulated frequencies",
            h4(strong("4) Calculated expected frequencies with each parameter
		      set")),
	    actionButton("go.exp", "Calculate expected response frequencies"),
	    plotOutput("plot.exp"),

	    br(),

	    h4(strong("5) Simulate response patterns with each parameter set")),
	    actionButton("go.sim", "Simulate possible response pattern"),
	    plotOutput("plot.sim")
	  )
	)
      ),
      
      tabItem(tabName = "quiz",
        fluidRow(
          box(status = "primary", title = "To which of the three parameter
	    estimation methods do these statements apply?",
            solidHeader = TRUE, width = 8,
	    
	    fluidRow(
	      column(6,
                radioButtons("q1", "____ estimation tends to inflate error rates",
	          c("Maximum Likelihood (ML)" = "ml", "Minimum Discrepancy (MD)" = "md",
		    "MDML" = "mdml", "Select one" = "null"),
	          selected = "null"
	        ),
		uiOutput("r1")
              ),

	      column(6,
                radioButtons("q2", "____ estimation is a reference point for
			     quantifying the amount of fit obtained by inflating
			     error rates, when compared to other methods",
	          c("Maximum Likelihood (ML)" = "ml", "Minimum Discrepancy (MD)" = "md",
                    "MDML" = "mdml", "Select one" = "null"),
	          selected = "null"
	        ),
		uiOutput("r2")
	      )
	    ),

	    fluidRow(
	      column(6,
                radioButtons("q3", "____ estimation produces explicit estimators,
			           making it computationally efficient",
	          c("Maximum Likelihood (ML)" = "ml", "Minimum Discrepancy (MD)" = "md",
		    "MDML" = "mdml", "Select one" = "null"),
	          selected = "null"
	        ),
		uiOutput("r3")
              ),

	      column(6,
                radioButtons("q4", "____ estimation has the best fit to the data",
	          c("Maximum Likelihood (ML)" = "ml", "Minimum Discrepancy (MD)" = "md",
                    "MDML" = "mdml", "Select one" = "null"),
	          selected = "null"
	        ),
		uiOutput("r4")
	      )
	    )
	  ),

	  box(status = "primary", title = "Your results",
	    verticalLayout(
              infoBoxOutput("sBox", width = 12), 	  
	      infoBoxOutput("aBox", width = 12), 
	      infoBoxOutput("cBox", width = 12)
	    ),
	  width = 4)
	)  
      )
    )
  )
)


############### SERVER ###############

server <- function(input, output){

#### NR input & plot ####

  sets1 <- c("0000", "1000", "0100", "0010", "0001", "1100",  "1010", "1001",
            "0110", "0101", "0011", "1110", "1101", "1011", "0111", "1111")
  sets2 <- as.pattern(as.binmat(sets1), as.letters = TRUE)
 

  # creates table with N.R values and the corresponding sets
  N.R <- eventReactive(input$go.plot, {
    nr <- c(input$NR1, input$NR2, input$NR3, input$NR4, input$NR5, input$NR6,
	    input$NR7, input$NR8, input$NR9, input$NR10, input$NR11,
	    input$NR12, input$NR13, input$NR14, input$NR15, input$NR16)

    setNames(nr, sets1)
  })

  # function for the barplot of the N.R frequencies
  NRplot <- function(NR){
    par(mar = c(0, 4, 3, 1)) 	  
    barplot(NR,   
            horiz     = TRUE, 
            names.arg = sets2,
            axes      = FALSE,
            cex.names = 1.3,
            las       = 2, 
            col       = "lightblue",
            xlim      = c(0, max(NR) + 2), 
	    border    = NA,
	    width = 1.4,
	    space = .2
	    )
    # include axis
    axis(side     = 3, 
         cex.axis = 1.3)
  }

  output$plot.NR <- renderPlot({
    NRplot(NR = N.R())
  })

#### display chosen knowledge structure ####
  
  # displays knowledge structure as text
  text.str <- eventReactive(input$go.K, {
    paste(as.pattern(as.binmat(input$Ks), as.letters = TRUE), collapse = "; ")
  })

  output$text.str <- eventReactive(input$go.K, {
    paste0("{", text.str(), "}")
  })

  output$text.corr <- eventReactive(input$go.K, {
    if(input$Ks[1] != "0000" | input$Ks[length(input$Ks)] != "1111"){
      HTML("<h4 style='color:red' align='left'> Caution! &#8709 and Q have
	   to be part of your knowledge structure!")
    } else {
      HTML("<h4 style='color:green' align='left'> Well done! Your knowledge
	   structure includes &#8709 and Q!")
    }
  })


  # convert input for knowledge structures into usable variable K
  K <- eventReactive(input$go.K, {
    matrix(as.numeric(unlist(strsplit(input$Ks, split = ""))), ncol = 4,
      byrow = TRUE, dimnames = list(input$Ks, c("a", "b", "c", "d")))
  })

#### parameter estimation with ML, MD & MDML ####
  blim.ML   <- eventReactive(input$go.para, {
      withProgress(message = "Calculating ML", value = 1, {
        blim(K(), N.R(), method = "ML")}
      )
  })

  blim.MD   <- eventReactive(input$go.para, {
    withProgress(message = "Calculating MD", value = 1, {
      blim(K(), N.R(), method = "MD")}
    )
  })

  blim.MDML <- eventReactive(input$go.para, {
    withProgress(message = "Calculating MDML", value = 1, {
      blim(K(), N.R(), method = "MDML")}
    )
  })

  
#### display estimate results and feedback messages ###
  get.data <- reactive({
      dat <- data.frame(
	   G2 = c(blim.ML()$goodness.of.fit[1], blim.MD()$goodness.of.fit[1],
		  blim.MDML()$goodness.of.fit[1]),
	   iterations = round(c(blim.ML()$iter, blim.MD()$iter,
				blim.MDML()$iter), digits = 0),
	   LogLikelihood = c(blim.ML()$loglik, blim.MD()$loglik,
		             blim.MDML()$loglik),
	   mean_careless_error = c(blim.ML()$nerror[1], blim.MD()$nerror[1],
			      blim.MDML()$nerror[1]),  # mean error
	   mean_lucky_guess = c(blim.ML()$nerror[2], blim.MD()$nerror[2],
			   blim.MDML()$nerror[2]),  # mean error
	   row.names = c("ML", "MD", "MDML")
	   )
      dat$iterations <- sprintf('%1.0f', dat$iterations) # is left-justified???
      dat
  })

  # display estimates in a table
  output$tab <- renderTable({
    validate(
      need(input$go.K > 0, "Please submit your knowledge structure"),
      need(input$go.plot > 0, "Please submit your response frequencies")
    )
    get.data()
  }, include.rownames = TRUE)

  output$text.iter <- eventReactive(input$go.para, {
    if(blim.ML()$iter > 9999 | blim.MDML()$iter > 9999){
      HTML("<h4 style='color:red' align='center'> Caution! Your estimations
	   probably did not converge! (10,000 iterations!)")
    } else {
      HTML("<h4 style='color:green' align = 'center'> Well done! Your
	   estimations converged! (<10,000 iterations!)")
    }
  })

  output$text.iden <- eventReactive(input$go.para, {
    if(dim(jacobian(blim.ML()))[2] == qr(jacobian(blim.ML()))$rank){
      HTML("<h4 style = 'color:green' align = 'center'> Well done! Your
	   parameters are locally identifiable!")
    } else {
      HTML("<h4 style = 'color:red' align = 'center'> Caution! Your parameters
	   are not locally identifiable! Consider respecifying your knowledge
	   structure.")
    }
  })

#### beta and eta plot functions ####

  betaPlot <- function(beta.ML, beta.MD, beta.MDML){

    beta <- matrix(round(c(rev(beta.ML), rev(beta.MD), rev(beta.MDML)),
			digits = 2),
		  byrow = TRUE, ncol = 4, nrow = 3)
   
    par(mai = c(.1, .6, .6, .3))

    barplot(beta, beside = TRUE,
            horiz     = TRUE, 
            names.arg = expression(beta[d], beta[c], beta[b], beta[a]),
            axes      = FALSE,
            cex.names = 2,
            las       = 2, 
            col       = c("olivedrab2", "olivedrab3", "olivedrab"),
            xlim      = c(0, 0.3),
	    legend    = c("ML", "MD", "MDML"),
	    border    = "white"
	    )
    # include axis
    axis(side     = 3, 
         cex.axis = 1.5, 
         at       = seq(0, 0.3, 0.1), 
         labels = as.character(seq(0, 0.3, 0.1)))
  }

  etaPlot <- function(eta.ML, eta.MD, eta.MDML){
    # eta.vector: vector with eta values (eta.a, eta.b, eta.c, eta.d)
    
    eta <- matrix(round(c(rev(eta.ML), rev(eta.MD), rev(eta.MDML)),
			digits = 2),
		  byrow = TRUE, ncol = 4, nrow = 3)

    par(mai = c(.1, .6, .6, .3))

    barplot(eta, beside = TRUE,
            horiz     = TRUE, 
            names.arg = expression(eta[d], eta[c], eta[b], eta[a]),
            axes      = FALSE,
            cex.names = 2,
            las       = 2,
            col       = c("olivedrab2", "olivedrab3", "olivedrab"),
            xlim      = c(0, 0.3),
	    legend    = c("ML", "MD", "MDML"),
	    border    = "white")
    axis(side = 3, cex.axis = 1.5, 
         at = seq(0, 0.3, 0.1), 
         labels = as.character(seq(0, 0.3, 0.1))
	 )
  }

  piPlot <- function(pi.ML, pi.MD, pi.MDML, pi.K){
    pis <- matrix(round(c(pi.ML, pi.MD, pi.MDML),
			digits = 2),
		  byrow = TRUE, ncol = length(pi.ML), nrow = 3)
    
    pi.ceil <- ceiling(max(pis*10))/10
  
    par(mai = c(.1, .6, .6, .3), mgp = c(3, .8, 0))

    barplot(pis, beside = TRUE,
            horiz     = TRUE,
	    names.arg = as.character(as.pattern(pi.K, as.letters = TRUE)), 
            axes      = FALSE,
            cex.names = 1,
            las       = 2,
            col       = c("olivedrab2", "olivedrab3", "olivedrab"),
            xlim      = c(0, pi.ceil),
	    legend    = c("ML", "MD", "MDML"),
	    border    = "white")
    axis(side = 3, cex.axis = 1.5, 
         at = seq(0, pi.ceil, 0.1), 
         labels = as.character(seq(0, pi.ceil, 0.1))
	 )
  }
    

  output$plot.beta <- renderPlot({
    etaPlot(blim.ML()$eta, blim.MD()$eta, blim.MDML()$eta)
  })

  output$plot.eta <- renderPlot({
    betaPlot(blim.ML()$beta, blim.MD()$beta, blim.MDML()$beta)
  })

  output$plot.pi  <- renderPlot({
    piPlot(blim.ML()$P.K, blim.MD()$P.K, blim.MDML()$P.K, blim.ML()$K)
  })

### plot functions: expected and simulated frequencies ###

  dat <- data.frame(set = sets1,
		    y.co = seq(1, by = 1.68, length.out = 16),
		    p.ML = 0, p.MD = 0, p.MDML = 0, sim.ML = 0, sim.MD = 0,
		    sim.MDML = 0) 
  dat <- dat[order(dat$set),]
  dat$ord <- as.numeric(rownames(dat))
 

  simPlot <- function(simML, simMD, simMDML){

    j <- 0
    for(i in 1:16){
      if(c(dat$set %in% names(simML))[i]){
        j <- j + 1
        dat$sim.ML[i] <- simML[j]
      } else { dat$sim.ML[i] <- 0 }
    }

    j <- 0
    for(i in 1:16){
      if(c(dat$set %in% names(simMD))[i]){
        j <- j + 1
        dat$sim.MD[i] <- simMD[j]
      } else { dat$sim.MD[i] <- 0 }
    }

    j <- 0
    for(i in 1:16){
      if(c(dat$set %in% names(simMDML))[i]){
        j <- j + 1
        dat$sim.MDML[i] <- simMDML[j]
      } else { dat$sim.MDML[i] <- 0 }
    }

    dat2 <- dat[order(dat$ord),]

    par(mar = c(0, 4, 3, 1)) 	  
    barplot(N.R(),   
            horiz     = TRUE, 
            names.arg = sets2,  
            axes      = FALSE,
            cex.names = 1.3,
            las       = 2, 
            col       = "lightblue",
            xlim      = c(0, max(simML, simMD, simMDML) + 5), 
	    border    = NA,
	    width = 1.4,
	    space = .2
	    )

    # include axis
    axis(side     = 3, 
         cex.axis = 1.3)

    points(x = dat2$sim.ML,
	   y = seq(1, by = 1.68, length.out = 16),
	   col = "red", pch = 19)

    points(x = dat2$sim.MD,
	   y = seq(1, by = 1.68, length.out = 16),
	   col = "blue", pch = 8)

    points(x = dat2$sim.MDML,
	   y = seq(1, by = 1.68, length.out = 16),
	   col = "green", pch = 11)

    legend("bottomright", c("observed", "ML", "MD", "MDML"),
	   pch = c(22, 19, 8, 11), col = c("black", "red", "blue", "green"),
	   pt.bg = c("lightblue", NA, NA), pt.cex = c(1.5, 1, 1, 1))
  }

  sim.ML <- eventReactive(input$go.sim, {
    simulate(blim.ML())
  })

  sim.MD <- eventReactive(input$go.sim, {
    simulate(blim.MD())
  })

  sim.MDML <- eventReactive(input$go.sim, {
    simulate(blim.MDML())
  })

  output$plot.sim <- renderPlot({
    validate(
      need(input$go.para > 0, "Please estimate your parameters")
    )

    simPlot(simML   = sim.ML(),
            simMD   = sim.MD(),
            simMDML = sim.MDML())
  })


  expPlot <- function(pkML, pkMD, pkMDML){

    dat2 <- dat[order(dat$ord),]	  
    j <- 0
    for(i in 1:16){
      if(c(dat2$set %in% names(pkML))[i]){
        j <- j + 1
        dat2$p.ML[i] <- pkML[j]
      } else { dat2$p.ML[i] <- 0 }
    }

    j <- 0
    for(i in 1:16){
      if(c(dat2$set %in% names(pkMD))[i]){
        j <- j + 1
        dat2$p.MD[i] <- pkMD[j]
      } else { dat2$p.MD[i] <- 0 }
    }

    j <- 0
    for(i in 1:16){
      if(c(dat2$set %in% names(pkMDML))[i]){
        j <- j + 1
        dat2$p.MDML[i] <- pkMDML[j]
      } else { dat2$p.MDML[i] <- 0 }
    }
    
    par(mar = c(0, 4, 3, 1)) 	  
    barplot(N.R(),  
      horiz     = TRUE, 
      names.arg = sets2, 
      axes      = FALSE,
      cex.names = 1.3,
      las       = 2, 
      col       = "lightblue",
      xlim      = c(0, max(c(dat2$p.ML, dat2$p.MD, dat2$p.MDML)) + 5), 
      border    = NA,
      width = 1.4,
      space = .2
    )
    
    # include axis
    axis(side     = 3, 
         cex.axis = 1.3)

    points(x = dat2$p.ML,
	   y = seq(1, by = 1.68, length.out = 16),
	   col = "red", pch = 19)

    points(x = dat2$p.MD,
	   y = seq(1, by = 1.68, length.out = 16),
	   col = "blue", pch = 8)

    points(x = dat2$p.MDML,
	   y = seq(1, by = 1.68, length.out = 16),
	   col = "green", pch = 11)

    legend("bottomright", c("observed", "ML", "MD", "MDML"),
	   pch = c(22, 19, 8, 11), col = c("black", "red", "blue", "green"),
	   pt.bg = c("lightblue", NA, NA), pt.cex = c(1.5, 1, 1, 1))
  }

  # expected frequencies: fitted(blim)
 
  pk.ML <- eventReactive(input$go.exp, {
    fitted(blim.ML())
  })

  pk.MD <- eventReactive(input$go.exp, {
    fitted(blim.MD())
  })

  pk.MDML <- eventReactive(input$go.exp, {
    fitted(blim.MDML())
  })

  output$plot.exp <- renderPlot({
    validate(
      need(input$go.para > 0, "Please estimate your parameters")
    ) 
    expPlot(pk.ML(), pk.MD(), pk.MDML())
  })

  # Reactions to questions
  pro <- reactiveValues(data = numeric(4), work = numeric(4))
  
  output$r1 <- renderText({
    if(input$q1 == "ml"){  # correct
      pro$data[1] <- 1
      pro$work[1] <- 1
      HTML("<h5 style = 'color:green' align = 'left'><b>Correct!</b></h5>")
    
    } else if(input$q1 == "null"){  # none selected
      pro$work[1] <- 0
      HTML("<br>")
    
    } else {  # wrong
      pro$data[1] <- 0
      pro$work[1] <- 1
      HTML("<h5 style ='color:red' align='left'><b>Wrong! Try again!</b></h5>")	 
    }
  })


  output$r2 <- renderText({
    if(input$q2 == "mdml"){  # correct
      pro$data[2] <- 1
      pro$work[2] <- 1
      HTML("<h5 style = 'color:green' align = 'left'><b>Correct!</b></h5>")
    
    } else if(input$q2 == "null"){  # none selected
      pro$work[2] <- 0
      HTML("<br>")
    
    } else {  # wrong
      pro$data[2] <- 0
      pro$work[2] <- 1
      HTML("<h5 style ='color:red' align='left'><b>Wrong! Try again!</b></h5>")	 
    }
  })


  output$r3 <- renderText({
    if(input$q3 == "md"){  # correct
      pro$data[3] <- 1
      pro$work[3] <- 1
      HTML("<h5 style = 'color:green' align = 'left'><b>Correct!</b></h5>")
    
    } else if(input$q3 == "null"){  # none selected
      pro$work[3] <- 0
      HTML("<br>")
    
    } else {  # wrong
      pro$data[3] <- 0
      pro$work[3] <- 1
      HTML("<h5 style ='color:red' align='left'><b>Wrong! Try again!</b></h5>")	 
    }
  })

  output$r4 <- renderText({
    if(input$q4 == "ml"){  # correct
      pro$data[4] <- 1
      pro$work[4] <- 1
      HTML("<h5 style = 'color:green' align = 'left'><b>Correct!</b></h5>")
    
    } else if(input$q4 == "null"){  # none selected
      pro$work[4] <- 0
      HTML("<br>")
    
    } else {  # wrong
      pro$data[4] <- 0
      pro$work[4] <- 1
      HTML("<h5 style ='color:red' align='left'><b>Wrong! Try again!</b></h5>")	 
    }
  })



  # BoxOutputs:
  #  How many questions answered correctly?
  output$sBox <- renderInfoBox({
    infoBox("", "Solved", paste(round(100*sum(pro$data)/4, 0), "%"),
      icon = icon("check", lib = "glyphicon"),
      color = "green"
    )
  })

  # How many questions attempted?
  output$aBox <- renderInfoBox({
    infoBox("", "Attempted", paste(round(100*sum(pro$work)/4, 0), "%"),
      icon = icon("pencil", lib = "glyphicon"),
      color = "yellow"
    )
  })

  # Percentage of correct answers
  output$cBox <- renderInfoBox({
    if(sum(pro$data)/sum(pro$work) == "NaN"){
      infoBox("", "Percentage correct", paste("/"),
        icon = icon("thumbs-up", lib = "glyphicon"),
        color = "blue"
      )
    } else {
      infoBox("", "Percentage correct",
	paste(round(sum(pro$data)/sum(pro$work)*100), "%"),
        icon = icon("thumbs-up", lib = "glyphicon"),
        color = "blue"
      )
    }
  })
}



shinyApp(ui = ui, server = server)