Example using R and ggplot of Measures of Risk and Error (MORE) plot for simulation output based on Nelson (2008)

MOREplotR.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# \"MORE Plot for hospital reception simulation\"\n",
    "- Louis Luangkesorn <[email protected]>\n",
    "- 8/11/2019"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This introduces the Measure of Risk and Error (MORE) plot [Nelson08] of simulation output which simultaneously shows the analyst a sense of the risk, expressed as a histogram along with the mean, the 0.05 quantile, and the 0.95 quantile), and a sense of the error of the simulation output (expressed as 95% C.I. for the mean, the 0.05 quantile, and the 0.95 quantile).\n",
    "\n",
    "This example uses a hospital receptionist simulation orginally found in [Nelson13]. The implementation of the simulation uses the Simmer library in R [simmer19]\n",
    "\n",
    "\n",
    "Nelson08: Nelson, B., (2008), THE MORE PLOT: DISPLAYING MEASURES OF RISK & ERROR FROM SIMULATION OUTPUT, 2008 Winter Simulation Conference, https://www.informs-sim.org/wsc08papers/048.pdf\n",
    "\n",
    "Nelson13: Nelson, B., (2013), Foundations and Methods of Stochastic Simulation, Springer, http://users.iems.northwestern.edu/~nelsonb/IEMS435/\n",
    "\n",
    "simmer19: Ucar I, Smeets B, Azcorra A (2019). “simmer: Discrete-Event Simulation for R.” Journal of Statistical Software, 90(2), 1–30. doi: 10.18637/jss.v090.i02. https://www.jstatsoft.org/article/view/v090i02/0 https://r-simmer.org/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "Attaching package: 'dplyr'\n",
      "\n",
      "The following objects are masked from 'package:stats':\n",
      "\n",
      "    filter, lag\n",
      "\n",
      "The following objects are masked from 'package:base':\n",
      "\n",
      "    intersect, setdiff, setequal, union\n",
      "\n",
      "\n",
      "Attaching package: 'simmer'\n",
      "\n",
      "The following object is masked from 'package:dplyr':\n",
      "\n",
      "    select\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(ggplot2)\n",
    "library(dplyr)\n",
    "library(magrittr)\n",
    "library(simmer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hospital simulation code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Global constants\n",
    "G <- list(RANDOM_SEED = 42,\n",
    "          INTERVAL_patientS = 10.0, # Generate new patients roughly every x seconds\n",
    "          maxTime = 8 * 60,          # minutes\n",
    "          timeReceptionist = 0.8,   # mean, minutes\n",
    "          phases = 3,\n",
    "          ARRint = 1.0             # mean, minutes\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visitor trajectory\n",
    "visitor <- trajectory(\"Arriving visitor\") %>%\n",
    "  seize(\"receptionist\", amount=1) %>%\n",
    "  timeout(function () rgamma(n=1, \n",
    "                             scale=G[['timeReceptionist']]/G[['phases']], \n",
    "                             shape = G[['phases']])\n",
    "  ) %>%\n",
    "  release(\"receptionist\", amount=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Main Simulation code\n",
    "reps = 10\n",
    "set.seed(G[['RANDOM_SEED']])\n",
    "envs <- lapply(1:reps, function(i) {\n",
    "                simmer(\"HospitalDES\") %>%\n",
    "                add_resource(\"receptionist\", 1) %>%\n",
    "                add_generator(\"visitor\", visitor, \n",
    "                              function() rexp(n=1,rate = (1/G[['ARRint']]))) %>%\n",
    "                run(until=G[['maxTime']])\n",
    "})\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Collect statistics with monitors\n",
    "receptionmonitor <- lapply(1:reps, function(i){\n",
    "  result <- envs[[i]] %>% \n",
    "  get_mon_resources()\n",
    "})\n",
    "patientmonitor <- lapply(1:reps, function(i){\n",
    "        result <- envs[[i]] %>% \n",
    "        get_mon_arrivals %>%\n",
    "        mutate(flow_time = end_time - start_time)\n",
    "      }\n",
    ")\n",
    "# Note that the result of lapply on the simmer output is a list of dataframes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Measures of Risk and Error\n",
    "\n",
    "##  Risk vs. Error\n",
    "-  *Measures of risk*\n",
    "\t-  Directly support decision making.\n",
    "\t-  What is the range of outcomes that can occur and their probabilities if the system is designed in a specific way.\n",
    "\t\t-  If $F_Y(5)$ (probability of completing a project in 5 days) is too small, we may decide not to undertake the project as designed.\n",
    "-  *Measures of error*\n",
    "\t-  Directly support experiment design\n",
    "\t-  Did we expend enough simulation effort (e.g. replications) to be confident in our estimates of system performance\n",
    "\t\n",
    "##  Measure of Risk and Error (MORE) Plot\n",
    "\n",
    "-  Measures of error.   The C.I. on the average or the quantiles are measure of error.\n",
    "    -  As a simulation analyst we can reduce error (how?)\n",
    "    -  How to choose which performance measure's C.I. we should use in specifying error.\n",
    "    -  What does the C.I. on these statistics (average, quantile) tell us about our system and experimental design?\n",
    "-  Measures of risk\n",
    "    -  How can we as engineers reduce risk?\n",
    "    -  How do we measure risk?  What do these types of statistics actually tell us about our system and experimental design?\n",
    "    \n",
    "## Creating a MORE plot\n",
    "\n",
    "-  Start with a histogram.\n",
    "-  Add the sample average and 0.05 and 0.95 quantiles.\n",
    "-  This range is considered *likely*, beyond that is considered *unlikely*.\n",
    "\t-  Measure of *risk*\n",
    "-  Include 95% C.I. for the average and both quantiles.\n",
    "\t-  Measure of *error* at the center and the extremes of the range of outcomes.\n",
    "-  Goal:  to represent both risk and error in an intuitively understandable way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading required package: foreach\n"
     ]
    }
   ],
   "source": [
    "# Creating MORE plot\n",
    "\n",
    "require(foreach)\n",
    "require(ggplot2)\n",
    "\n",
    "moreplot <- function(simoutput, dataelement){\n",
    "  # Measure of Risk and Error (MORE) plot from Nelson, et. al.\n",
    "  # Generate histogram using all of the simulation data from all replications\n",
    "  # Add vertical lines with mean and 5% and 95% quantile and error bars for confidence intervals from the replications\n",
    "  # INPUTS: \n",
    "  # - simoutput: list of dataframes of simmer monitor output from simulation model replications\n",
    "  # - dataelement:  quoted name of column with performance data to be plotted\n",
    "  # OUTPUT: \n",
    "  # creates a MORE plot with histogram of performance measure, mean, and 0.05 and 0.95 quantiles with 95% confidence intervals\n",
    "  # EXAMPLE OF USE: \n",
    "  # moreplot(patientmonitor, flow_time)\n",
    "\n",
    "  reps <- length(simoutput)\n",
    "  allpatients <- bind_rows(simoutput)\n",
    "  dataindex <- which(colnames(allpatients)==dataelement)\n",
    "  # get mean and averages for mean and the 5% and 95% quantiles \n",
    "  getquantiles <- function(datavector){\n",
    "    mean <- mean(datavector)\n",
    "    quantiles <- quantile(datavector, probs=c(0.05, 0.95))\n",
    "    return(c(mean, quantiles[[1]], quantiles[[2]]))\n",
    "  }\n",
    "  \n",
    "  # get mean and quantiles for each replication, \n",
    "  # then calculate the mean and C.I. for the mean and both quantiles\n",
    "  \n",
    "  samplestats <- foreach(i=1:length(simoutput)) %do% getquantiles(simoutput[[i]][[dataelement]])  \n",
    "  samplestatsdf <- as.data.frame(t(as.data.frame(lapply(samplestats, as.data.frame))))\n",
    "  names(samplestatsdf) <- c('mean', 'q05', 'q95')\n",
    "  \n",
    "  # create a data frame for the ci of each statistic for use in the ggplot to follow\n",
    "  meanstats <- data.frame(mean = mean(samplestatsdf$mean),\n",
    "                          lowerci = qt(0.05,reps) * sd(samplestatsdf$mean)/length(samplestatsdf$mean) + \n",
    "                                      mean(samplestatsdf$mean),\n",
    "                          upperci = qt(0.95, reps) * sd(samplestatsdf$mean)/length(samplestatsdf$mean) + \n",
    "                                      mean(samplestatsdf$mean))\n",
    "  q05stats <- data.frame(mean = mean(samplestatsdf$q05),\n",
    "                         lowerci = qt(0.05, reps) * sd(samplestatsdf$q05)/length(samplestatsdf$q05) + \n",
    "                                    mean(samplestatsdf$q05),\n",
    "                         upperci = qt(0.95, reps) * sd(samplestatsdf$q05)/length(samplestatsdf$q05) + \n",
    "                                    mean(samplestatsdf$q05))  \n",
    "  q95stats <- data.frame(mean = mean(samplestatsdf$q95),\n",
    "                         lowerci = qt(0.05, reps) * sd(samplestatsdf$q95)/length(samplestatsdf$q95) + \n",
    "                                     mean(samplestatsdf$q95),\n",
    "                         upperci = qt(0.95, reps) * sd(samplestatsdf$q95)/length(samplestatsdf$q95) + \n",
    "                                     mean(samplestatsdf$q95))\n",
    "  # Generate base histogram\n",
    "  flowplot <- ggplot(allpatients) +\n",
    "    geom_histogram(aes(x=flow_time, y=..density..)) +\n",
    "    ggtitle(\"Measure of risk and error plot\") +\n",
    "    xlab(paste(dataelement, \"(q05, mean, q95)\"))\n",
    "  ggp_build <- ggplot_build(flowplot)\n",
    "  maxdensity <- max(ggp_build[[\"data\"]][[1]][[\"density\"]])\n",
    "  # Add vertical lines for mean, q05, q95, along with the C.I. error bars\n",
    "  flowplot <- flowplot +\n",
    "    geom_segment(data = meanstats, aes(x = mean, y = 0, xend = mean, yend = (ceiling(maxdensity*100)/100))) +\n",
    "    geom_segment(data = q05stats, aes(x = mean, y = 0, xend = mean, yend = (ceiling(maxdensity*100)/100))) +\n",
    "    geom_segment(data = q95stats, aes(x = mean, y = 0, xend = mean, yend = (ceiling(maxdensity*100)/100))) +\n",
    "    geom_errorbarh(data = meanstats, aes(xmax = upperci, xmin = lowerci, \n",
    "                                         y = -0.01, height = 0.01)) +\n",
    "    geom_errorbarh(data = q05stats, aes(xmax = upperci, xmin = lowerci, \n",
    "                                        y = -0.01, height = 0.01)) +\n",
    "    geom_errorbarh(data = q95stats, aes(xmax = upperci, xmin = lowerci, \n",
    "                                        y = -0.01, height = 0.01)) +\n",
    "    annotate(\"text\", label = \"C.I. mean\", x = meanstats$mean, y = -0.02) +\n",
    "    annotate(\"text\", label = \"C.I. Q0.05\", x = q05stats$mean, y = -0.02) +\n",
    "    annotate(\"text\", label = \"C.I. Q0.95\", x = q95stats$mean, y = -0.02)\n",
    "  flowplot\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.\n",
      "`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAYFBMVEUAAAAzMzNHR0dNTU1Z\nWVlgYGBoaGhycnJ8fHyBgYGMjIyOjo6ampqkpKSnp6eurq6ysrK3t7e9vb3AwMDHx8fIyMjP\nz8/Q0NDW1tbZ2dnd3d3h4eHp6enr6+vw8PD////sZuebAAAACXBIWXMAABJ0AAASdAHeZh94\nAAAgAElEQVR4nO2di3bayLZF1eGGOI4P7XbS7Th2zP//5ZW0eRQSotiqBaZ25hrjJLCBqeWS\nphGPnG7WhJDiNB9dgJAIQSRCBEEkQgRBJEIEQSRCBEEkQgRBJEIEQSRCBEEkQgRBJEIEKRNp\ncf/0trn49nS/KK9zPO8PTbM8HDXN6euZOO/uIb0/qNCkopQdUE3TrDYXV43u4BzmvmXfDbd8\n+nomFxTpcstAbjilIi23T0OL5eWOoKZ5zd/FSZxdJkdCpD8ypSI9Ni/9pZf20gVFUtyl4O4e\nEiL9kSkV6W1zbrdqXjdH0Mtdex72bLev2qepOzPttX2ds3h4We+OtO6vpnlfNveDB60391+0\n93/t79gkp42bh9hgAH1bLB93gP2229tWi2Zp+Oe7ZvF4cLTvt71jb1rtO+xnG+DTYveYI23n\nriepNqUirTfndovF5pD6Ycd9r9fCLncH3Euzu3wg0n1/1/RBluf0/oci9Q/pB0PoInkllWy7\n9WR3edVfekiO9WTbO7a1SjrsZpsOj8ljjrQtWlRSY4pFWvXndi/bY/u1fxZ4veuOqceme354\n7I/uZfNj3R20y4FId++DB1nap7fH9/V7e9i/rQ/PljYP6UcD6F3iUbrt9lnrZf1+311uue1D\nnhd7YrrtHXvTKumwmW07LIzyahse3LNoSUmdKRbppf+93OnUH0Grpj/e3rvToKVd7ufJ4XUg\nUn/ulT7Istr89n/Y/9bfPfplBzmErprF2+764bafe3rTcZ/6LS12j0y3vWO/jDpsZtuN2YYW\nK9vC4J5nLh6JlGKRupO67kRq9ySxTX+H1+fHu/7iffsS4sfb7kHrrUj9xcGD+ond963//OhQ\npP1fB9Cn/skmyW7byUM2frWP3G9pv+2UPehwoMf2me9+uaVO3ZP8KSkXqXsyemkedr/9Eyee\nFruLb/3F5dP6mEjNSKTdpeEzz8HBfgBtDj+0TbY93mB7FrjnTYl00OFQpM2Gls3gxtE9yZ+S\ncpGe2/OZ1e71fnoYtc8Rd6sfb5vR80N3tD4eF+kYeH9hSqQD6OJHs3+NdbDtIyIdUXbEzot0\nh0hkk3KR1s2iP7PbvnG2/+x0mbyc6fP60N73yHG9GH3get6p3QH0rX2JskwAo5dSyQbf9w9P\nt332qd3mtRyndmQbgUgPzfPmzbH+Sv9Vs9duYofU81CDpn+d8pwc1+mDLCubHH2zIf0rgfZ6\nvBze7Xko0r1p87R/+Ljw9q+DDocibd6g2L/ZMHVP8qdEIFL3IcqP7ZXX/tOV10V3lrXs3iLb\nvNVs71Svut/ld839+2a8fUsieZClPSdb2RvKr+tpkRJof321f9sv3XbykKf+jesfyauxdNuH\nIh10GIi0eO7p73bD4J5va/LHRSBS+0TQ7D7a2X422T2RPG1fxr/sPjvt3qC2i6v09UTyoPXB\npEm+PZBscvNXArUm+/O0dNvpQ+yj2fTrTMm2B892hx/IHvzY94c3pPdc8oHsnxiBSNsPQre/\nyLuv4/Sf1nTvnC0eXp7tS0D9t3n6X9Yvy+23dHaHXPKgTZKv50yf2u2hdv1p/5SUbPvgIT9G\nXxHab3t42njwFaHDH7v9PZDekNzzddlc7N+TkJsNvzxnhKccMgyHxIwgEhmGQ2JGEIkMwyEx\nI4hEhuGQIEQQRCJEEEQiRBBEIkQQRCJEEEQiRBBEIkSQEpF+++K9/8di4cLNcxEJLlwBF5Hg\nwhVwEQkuXAEXkeDCFXARCS5cAReR4MIVcBEJLlwBF5HgwhVwEQkuXAEXkeDCFXARCS5cAReR\n4MIVcBEJLlwBF5HgwhVwEQkuXAEXkeDCFXARCS5cAReR4MIVcBEJLlwBF5HgwhVwEQkuXAEX\nkeDCFXARCS5cAReR4MIVcBEJLlwBF5HgwhVwEQkuXAEXkeDCFXARCS5cAReR4MIVcBEJLlwB\nF5HgwhVwEQkuXAH3FkRqZv5AGayMft0dXVtfuDZEpJnc8sToC9eGiDSTW54YfeHaEJFmcssT\noy9cGyLSTG55YvSFa0NEmsktT4y+cG2ISDO55YnRF64NEWkmtzwx+sK1ISLN5JYnRl+4NkSk\nmdzyxOgL14aINJNbnhh94doQkWZyyxOjL1wbItJMbnli9IVrQ0SayS1PjL5wbYhIM7nlidEX\nrg0RaSa3PDH6wrUhIs3klidGX7g2RKSZ3PLE6AvXhog0k1ueGH3h2hCRZnLLE6MvXBsi0kxu\neWL0hWtDRJrJLU+MvnBtiEgzueWJ0ReuDW9UpP87Hv8PWN2BWVtfuDZEpGwQCW6Wi0j5IBLc\nLBeR8kEkuFkuIuWDSHCzXETKB5HgZrmIlA8iwc1yESkfRIKb5X60SJ0dzYQ1iKTklgfuCS4i\nuepKE6MvXBsikqeuNDH6wrUhInnqShOjL1wbIpKnrjQx+sK1ISJ56koToy9cGyKSp640MfrC\ntSEieepKE6MvXBsikqeuNDH6wrUhInnqShOjL1wbIpKnrjQx+sK1ISJ56koToy9cGyKSp640\nMfrCtSEieepKE6MvXBsikqeuNDH6wrUhInnqShOjL1wbSkSaH6dIM7bwMT/X/NTWlwzCM5KP\nW54YfeHaEJE8daWJ0ReuDRHJU1eaGH3h2hCRPHWlidEXrg0RyVNXmhh94doQkTx1pYnRF64N\nEclTV5oYfeHaEJE8daWJ0ReuDRHJU1eaGH3h2hCRPHWlidEXrg0RyVNXmhh94doQkTx1pYnR\nF64NEclTV5oYfeHaEJE8daWJ0ReuDRHJU1eaGH3h2hCRPHWlidEXrg0RyVNXmhh94doQkTx1\npYnRF64NEclTV5oYfeHaEJE8daWJ0ReuDRHJU1eaGH3h2hCRPHWlidEXrg0RyVNXmhh94doQ\nkTx1pYnRF64NEclTV5oYfeHaEJE8daWJ0ReuDRHJU1eaGH3h2hCRPHWlidEXrg0RyVNXmhh9\n4doQkTx1pYnRF64NEclTV5oYfeHaEJE8daWJ0ReuDRHJU1eaGH3h2hCRPHWlidEXrg0RyVNX\nmhh94doQkTx1pYnRF64NEclTV5oYfeHaEJE8daWJ0ReuDRHJU1eaGH3h2hCRPHWlidEXrg0R\nyVNXmhh94doQkTx1pYnRF64NEclTV5oYfeHaEJE8daWJ0ReuDRHJU1eaGH3h2hCRPHWlidEX\nrg0RyVNXmhh94doQkTx1pYnRF64NEclTV5oYfeHaEJE8daWJ0ReuDRHJU1eaGH3h2hCRPHWl\nidEXrg0RyVNXmhh94doQkTx1pYnRF64NEclTV5oYfeHaEJE8daWJ0ReuDRHJU1eaGH3h2hCR\nPHWlidEXrg0RyVNXmhh94doQkTx1pYnRF64NEclTV5oYfeHaEJE8daWJ0ReuDRHJU1eaGH3h\n2hCRPHWlidEXrg0RyVNXmhh94doQkTx1pYnRF64NEclTV5oYfeHaEJE8daWJ0ReuDRHJU1ea\nGH3h2hCRPHWlidEXrg0RyVNXmhh94doQkTx1pYnRF64NEclTV5oYfeHaEJE8daWJ0ReuDRHJ\nU1eaGH3h2hCRPHWlidEXrg0RyVNXmhh94doQkTx1pYnRF64NEclTV5oYfeHaEJE8daWJ0Reu\nDRHJU1eaGH3h2hCRPHWlidEXrg0RyVNXmhh94doQkTx1pYnRF64NEclTV5oYfeHaEJE8daWJ\n0ReuDRHJU1eaGH3h2hCRPHWlidEXrg0RyVNXmhh94doQkTx1pYnRF64NEclTV5oYfeHaEJE8\ndaWJ0ReuDRHJU1eaGH3h2hCRPHWlidEXrg0RyVNXmhh94doQkTx1pYnRF64NEclTV5oYfeHa\nEJE8daWJ0ReuDRHJU1eaGH3h2hCRPHWlidEXrg0RyVNXmhh94doQkTx1pYnRF64NEclTV5oY\nfeHaEJE8daWJ0ReuDRHJU1eaGH3h2hCRPHWlidEXrg0RyVNXmhh94doQkTx1pYnRF64NEclT\nV5oYfeHaEJE8daWJ0ReuDRHJU1eaGH3h2hCRPHWlidEXrg0dIi3aHLuMSEpubX3h2vB8kRa7\nPw4vI5KUW1tfuDZEJE9daWL0hWvDeSKtDy/Pb4JIw9TWF64NS0X6q0vuodNxijRjCyVvonxE\nautLBnGLxJsNF+HW1heuDRHJU1eaGH3h2nC2SOlJ3vwmiDRMbX3h2nCuSKlH1xPppF8xDsza\n+sK14UyRDjxCJCG3tr5wbXi+SLtvMyzs4uID3v5GpNnc8sA9wfWINJ35TRBpmNr6wrUhInnq\nShOjL1wbIpKnrjQx+sK1ISJ56koToy9cGyKSp640MfrCtSEieepKE6MvXBsikqeuNDH6wrUh\nInnqShOjL1wbIpKnrjQx+sK1ISJ56koToy9cGyKSp640MfrCtSEieepKE6MvXBsikqeuNDH6\nwrUhInnqShOjL1wbIpKnrjQx+sK1ISJ56koToy9cGyKSp640MfrCtSEieepKE6MvXBsikqeu\nNDH6wrUhInnqShOjL1wbIpKnrjQx+sK1ISJ56koToy9cGyKSp640MfrCtSEieepKE6MvXBsi\nkqeuNDH6wrUhInnqShOjL1wbIpKnrjQx+sK1ISJ56koToy9cGyKSp640MfrCtSEieepKE6Mv\nXBsikqeuNDH6wrUhInnqShOjL1wbIpKnrjQx+sK1ISJ56koToy9cGyKSp640MfrCtSEieepK\nE6MvXBsikqeuNDH6wrUhInnqShOjL1wbIpKnrjQx+sK1ISJ56koToy9cGyKSp640MfrCtSEi\neepKE6MvXBsikqeuNDH6wrUhInnqShOjL1wbIpKnrjQx+sK1ISJ56koToy9cGyKSp640MfrC\ntSEieepKE6MvXBsikqeuNDH6wrUhInnqShOjL1wbIpKnrjQx+sK1ISJ56koToy9cGyKSp640\nMfrCtSEieepKE6MvXBsikqeuNDH6wrUhInnqShOjL1wbIpKnrjQx+sK1ISJ56koToy9cGyKS\np640MfrCtSEieepKE6MvXBsikqeuNDH6wrUhInnqShOjL1wbIpKnrjQx+sK1ISJ56koToy9c\nGyKSp640MfrCtSEieepKE6MvXBsikqeuNDH6wrUhInnqShOjL1wbIpKnrjQx+sK1ISJ56koT\noy9cGyKSp640MfrCtSEieepKE6MvXBuGEOmUYLUdmLX1hWtDRPLUlQaRAnERyVVXGkQKxEUk\nV11pECkQF5FcdaVBpEBcRHLVlQaRAnERyVVXGkQKxEUkV11pECkQF5FcdaVBpEBcRHLVlQaR\nAnERyVVXGkQKxEUkV11pECkQF5FcdaVBpEBcRHLVlQaRAnERyVVXGkQKxEUkV11pECkQF5Fc\ndaVBpEBcRHLVlQaRAnERyVVXGkQKxEUkV11pECkQF5FcdaVBpEBcRHLVlQaRAnERyVVXGkQK\nxEUkV11pECkQF5FcdaVBpEBcRHLVlQaRAnERyVVXGkQKxEUkV11pECkQF5FcdaVBpEBcRHLV\nlQaRAnERyVVXGkQKxEUkV11pECkQF5FcdaVBpEBcRHLVlQaRAnERyVVXGkQKxEUkV11pECkQ\nF5FcdaVBpEBcRHLVlQaRAnE1Is3PZUX6uJ9rfmrrSwbhGemMutLwjBSIi0iuutIgUiAuIrnq\nSoNIgbiI5KorDSIF4iKSq640iBSIi0iuutIgUiAuIrnqSoNIgbiI5KorDSIF4iKSq640iBSI\ni0iuutIgUiAuIrnqSoNIgbiI5KorDSIF4iKSq640iBSIi0iuutIgUiAuIrnqSoNIgbiI5Kor\nDSIF4iKSq640iBSIi0iuutIgUiAuIrnqSoNIgbiI5KorDSIF4iKSq640iBSIi0iuutIgUiAu\nIrnqSoNIgbiI5KorDSIF4iKSq640iBSIi0iuutIgUiAuIrnqSoNIgbiI5KorDSIF4iKSq640\niBSIi0iuutIgUiDupEjLxxdEGtWVBpECcSdFappm8fCMSL8RCe4Z3EmR3n/cty41dz/eEKmY\n4OHW1heuDadE6vK8WrQuLfPPS/ObINIwtfWFa8NTIq3fVk3/tIRIlwgiBeKeEun1vn86erlr\n7hHpAkGkQNxpkZ7vdmd1Te6t8flNEGmY2vrCteGUSMumuX/d3rRApAsEkQJxJ0VqVq/rszO/\nCSINU1tfuDacEun9fI0QScmtrS9cG04+I22uL3KndYik5dbWF64Nj4q0aJIgUjHBw62tL1wb\nHhXpKfHoCZGKCR5ubX3h2vCoSOsz3vJGpMtwa+sL14ZTIrkyvwkiDVNbX7g2PCpS+3TEa6Rj\ndaVBpEBcRHLVlQaRAnGPi+TN/CaINExtfeHaEJE8daVBpEDcaZGeFuv1S7N4RCREgpvlTor0\n1L44eus+mD3HpPlNEGmY2vrCteGUSMvmpf3f02v2m9+IpOXW1heuDadEap+QnpvlmR/Mzm+C\nSMPU1heuDadEWjRvD81r9yoJkYoJHm5tfeHacEqkx+7/j6t7QlohUjHBw62tL1wbTom0XjWL\n5/aJ6RyPEEnIra0vXBtOiuTJ/CaINExtfeHaEJE8daVBpEDcaZFWC75rN6orDSIF4k6KtOJL\nq0fqSoNIgbiTIi3O+qexiCTn1tYXrg2nROJfyB6rKw0iBeJOinTfOP4PueY3QaRhausL14ZT\nIr0t7s75D7ogkppbW1+4NpwSiX8he6yuNIgUiItIrrrSIFIg7qRIrsxvgkjD1NYXrg0RyVNX\nGkQKxD0h0tN9e1p3d9Z/k2J+E0Qapra+cG04JdL7sn991DQviFRM8HBr6wvXhlMiPTSr7kPZ\nH/n/gOwNizRBVyycIIgUiDspUvdu3fZ/iHSRIFIgLiLNXDhBECkQd1KkzandqnlAJOf9y7iI\nVCV3UqT3zT9HWpzzRaH5TRBpGESqkjsp0nr9uGya5eqsr67Ob4JIwyBSldwTIjkyvwkiDYNI\nVXIRaebCCYJIgbjHRUr/W8y8a4dIcLNcRJq5cIIgUiDucZG63Hf/sO/t7v4MjxBpVhApEHdS\npO0/NW/OMWl+E0QaBpGq5E6KtDmle+fUDpHg5rmTIt01dmrHMxIiwc1zJ0V645sNJxdOEEQK\nxJ0Uaf2+6r7Z8Mg3GxAJbp47LZIn85sg0jCIVCX3jxRpKp6FEwSRAnERKYln4QRBpEBcREri\nWThBECkQF5GSeBZOEEQKxEWkJJ6FEwSRAnERKYln4QRBpEBcREriWThBECkQF5GSeBZOEEQK\nxEWkJJ6FEwSRAnERKYln4QRBpEBcREriWThBECkQF5GSeBZOEEQKxEWkJJ6FEwSRAnERKYln\n4QRBpEBcREriWThBECkQF5GSeBZOEEQKxEWkJJ6FEwSRAnERKYln4QRBpEBcREriWThBECkQ\nF5GSeBZOEEQKxEWkJJ6FEwSRAnERKYln4QRBpEBcREriWThBECkQF5GSeBZOEEQKxEWkJJ6F\nEwSRAnE9Ii3aJNcQqTCIFIjrEGlxoM8CkUqDSIG4c0Va8IxUHEQKxJ39jIRIxUGkQNxikf7q\nknvodJyH+ow46PN/DGFK3vQhNxCekXhGgjubi0hJPAsnCCIF4iJSEs/CCYJIgbiIlMSzcIIg\nUiAuIiXxLJwgiBSI6xBp982GBSJJgkiBuB6RpjO/ifNQnxFEUgXuCS4iJfEsnCCIFIiLSEk8\nCycIIgXiIlISz8IJgkiBuIiUxLNwgiBSIC4iJfEsnCCIFIiLSEk8CycIIgXiIlISz8IJgkiB\nuIiUxLNwgiBSIC4iJfEsnCCIFIiLSEk8CycIIgXiIlISz8IJgkiBuIiUxLNwgiBSIC4iJfEs\nnCCIFIiLSEk8CycIIgXiIlISz8IJgkiBuIiUxLNwgiBSIC4iJfEsnCCIFIiLSEk8CycIIgXi\nIlISz8IJgkiBuIiUxLNwgiBSIC4iJfEsnCCIFIiLSEk8CycIIgXiIlISz8IJgkiBuIiUxLNw\ngiBSIC4i5YNIcLNcRMoHkeBmuYiUDyLBzXIRKR9EgpvlIlI+iAQ3y0WkfBAJbpaLSPkgEtws\nF5HyQSS4WS4i5YNIcLNcRMoHkeBmuYiUDyLBzXIRKR9EgpvlIlI+iAQ3y0WkfBAJbpaLSPkg\nEtwsF5HyQSS4WS4i5YNIcLNcRMoHkeBmuYiUDyLBzXIRKR9EgpvlIlI+iAQ3y0WkfBAJbpaL\nSPkgEtwsF5HyQSS4WS4i5YNIcLNcRMoHkeBmuYiUDyLBzXIRKR9EgpvlIlI+iAQ3y0WkfBAJ\nbpaLSPkgEtwsF5HyQSS4WS4i5YNIcLNcRMoHkeBmuYiUDyLBzXIRKR9EgpvlIlI+iAQ3y0Wk\nfBAJbpaLSPkgEtwsF5HyQSS4WS4i5YNIcLNcRMoHkeBmuYiUDyLBzXIRKR9EgpvlIlI+iAQ3\ny0WkfBAJbpaLSLNTsCtGy5AEkarkItLsFOyK0TIkQaQquYg0OwW7YrQMSRCpSi4izU7Brhgt\nQxJEqpKLSLNTsCtGy5AEkarkItLsFOyK0TIkQaQquYg0OwW7YrQMSRCpSi4izU7BrhgtQxJE\nqpKLSLNTsCtGy5AEkarkItLsFOyK0TIkQaQquYg0OwW7YrQMSRCpSi4izU7BrhgtQxJEqpKL\nSLNTsCtGy5AEkarkItLsFOyK0TIkQaQquYg0OwW7YrQMSRCpSi4izU7BrhgtQxJEqpKLSLNT\nsCtGy5AEkarkItLsFOyK0TIkQaQquYg0OwW7YrQMSRCpSi4izU7BrhgtQxJEqpKLSLNTsCtG\ny5AEkarkItLsFOyK0TIkQaQquYg0OwW7YrQMSRCpSi4izU7BrhgtQxJEqpKLSLNTsCtGy5AE\nkarkItLsFOyK0TIkQaQquYg0OwW7YrQMSRCpSi4izU7BrhgtQxJEqpKLSLNTsCtGy5AEkark\nItLsFOyK0TIkQaQquYg0OwW7YrQMSRCpSi4izU7BrhgtQxJEqpKLSLNTsCtGy5AEkarkItLs\nFOyK0TIkQaQquYg0OwW7YrQMSRCpSi4izU7BrhgtQxJEqpKLSLNTsCtGy5AEkarkItLsFOyK\n0TIkQaQquYg0OwW7YrQMSRCpSq5GpPm56KH+fxelX2hFPmY/EFl4RnKm4HfaaBmS8IxUJReR\nZqdgV4yWIQkiVclFpNkp2BWjZUiCSFVyEWl2CnbFaBmSIFKVXESanYJdMVqGJIhUJReRZqdg\nV4yWIQkiVclFpNkp2BWjZUiCSFVyEUmesh2CSHVyEUmesh2CSHVyEUmesh2CSHVyEUmesh2C\nSHVyEUmesh2CSHVyrybSiSMPkdIgUpVcRJKnbIcgUp1cRJKnbIcgUp1cRJKnbIcgUp1cRJKn\nbIcgUp1cRJKnbIcgUp1cRJKnbIcgUp1cRJKnbIcgUp1cRJKnbIcgUp1cRJKnbIcgUp1cRJKn\nbIcgUp1cRJKnbIcgUp1cRJKnbIcgUp1cRJKnbIcgUp1cRJKnbIcgUp1cRJKnbIcgUp1cRJKn\nbIcgUp1cRJKnbIcgUp1cRJKnbIcgUp1cRJKnbIcgUp1cRJKnbIcgUp1cRJKnbIcgUp1cRJKn\nbIcgUp1cRJKnbIcgUp1cRJKnbIcgUp1cRJKnbIcgUp1cRJKnbIcgUp1cRLpWztwhiFQnF5Gu\nlTN3CCLVyUWka+XMHYJIdXIR6Vo5c4cgUp1cRLpWztwhiFQnF5GulTN3CCLVyUWka+XMHYJI\ndXIR6Vo5c4cgUp1cRLpWztwhiFQnF5GulTN3CCLVyUWka+XMHYJIdXIR6Vo5c4cgUp1cRLpW\nztwhiFQnF5GulTN3CCLVyUWka+XMHYJIdXIR6Vo5c4cgUp1cRLpWztwhiFQnF5GulTN3CCLV\nyUWkD85ooRCpSi4ifXBGC4VIVXIR6YMzWihEqpKLSB+c0UIhUpVcRPrgjBYKkarkItIHZ7RQ\niFQlF5E+OKOFQqQquYj0wRktFCJVyUWkD85ooRCpSi4ifXBGC4VIVXIR6YMzWihEqpKLSB+c\n0UIhUpVcRPrgjBYKkarkItIHZ7RQiFQlF5E+OKOFQqQquYj0wRktFCJVyUWkD85ooRCpSi4i\n3Vo2fdU7WhC4J7iIdGtBpCq5iHRrQaQquYh0a0GkKrmIdGtBpCq5iHRrQaQquYh0a0GkKrmI\ndGtBpCq5iHRrQaQquYh0a0GkKrmIdGtBpCq5iHRrQaQquYh0a0GkKrmIdGtBpCq5iHRrOd13\n9o4WBO4JLiLdWhCpSi4i3VoQqUouIt1aEKlKLiLdWhCpSi4i3VoQqUouIt1aEKlKLiLdWhCp\nSi4i3Vrm9c3uaEHgnuAi0q0FkarkItKt5Y8T6YyfaF4QSRpEUgWRTnAR6daCSOOfaF4QSRpE\nUgWRTnAR6daCSOOfaF4QSRpEUuW6Ip3xg35MX0SqI9K+lz+AfLv6qj8RIkmDSBc9gHy7+qo/\nESJJg0gXPYB8u/pyKejrDSJVkXpEUhYtzhl9L7UONkSkGwsizcoZfS+1DjZEpBvLDYqkbHSp\nHDu2NT/9eVxEurUg0qwcO7Y1P/15XES6tVyj74yddOs594AXpFikRZtjlxFJmdr63nbEDpkz\nR4fni7TY/XF4GZGkqa3vbcd5OCLSR9H1qa3vbQeRZKntwKyt721nxuGYpVxGpL+65B5KyB+T\nCz8jZZUuT2XfeoYbiYtIcOEKuIgEF66Ai0hw4Qq4iAQXroDrEGn3bYZFchmR4ML9fcXv2mWb\nlCfEDoFbJxeR4MIVcBEJLlwBF5HgwhVwEQkuXAH3FkRqmpk/0Wns7+boxXKuJI2be/bPEOLA\nrI2LSH6uJIgUi4tIfq4kiBSLi0h+riSIFIuLSH6uJIgUi4tIfq4kiBSLi0h+riSIFIuLSH6u\nJIgUi4tIfq4kiBSLi0h+riSIFIuLSH6uJIgUi4tIfq4kiBSLG1qkJEquJI0/H9oX7mluaJGO\nXiznSsIzUiwuIvm5kiBSLC4i+bmSIFIsLiL5uZIgUiwuIvm5kiBSLC4i+bmSIBSsZpkAAAkm\nSURBVFIsLiL5uZIgUiwuIvm5kiBSLC4i+bmSIFIsLiL5uZIgUizuLYikSogdArdOLiLBhSvg\nIhJcuAIuIsGFK+AiEly4Ai4iwYUr4CISXLgC7oVF+vtL03z5Z3d184nR+vD2v+3yv+3lr//u\nbtpfPfXvQ9MNrHUfSR1EtkOOr8bg9uOr8d+X5tPDL3tU7l/Lhjgwa+NeVKSfn2yvf9kOBiL9\n/Gy3f/rZXvlll//b3La/+t/0oXO4gRsXaWo1drefWI3/drecWA1xX7ge7kVF+tR8a3f990/N\n5rfsUKTPzZd/+9+9n9orfzcPv39/a/63uW1/9Xt36ZwN3LhIU6uxzanV+NY9qJ+dWA1xX7ge\n7iVF+qf52v/9vT80uhyK9H37y/lLd5h8tWefr5u77q/+3ezPhk5uYCRSe+1re5f2V/3X7qzo\n17em+dafHv37tf31/tDf4+dXu3Qimh0yuRqbnFwNu2t3j+nVEPeF6+JeUqSvzeYcf3uCMhBp\nd3t/DH1qknukV7813782n3eM6Q0cEakVpvmnO2X61iPbfO631+ehu8enzaUT0eyQydUY3n5s\nNZozVkPcF66Le0mRxudZhyLtb+8uNVOHzlc77H9mN3BEpG/tE0GryT/d/H+dLw/dr/vP3W/1\n/2yjX361v+Q//T4VzQ6ZXI3RtSOr0T9LPZxejV1CHJi1cSsQqWm+d2c037IbOCLSz+6PXzb/\nbMT+bOnn9/99MfbPo0WHXEGKRPq3f5fi9GrsEuLArI17IyJ9PiHSFOwckdI/kv+fxS/bS8Ot\nHM3VRRqvxu/v3eu806uxS4gDszbuJUXanfX/3v498RrpX3u50t/j8JX4/qTr2KEz2MDZIn1r\nPv/9/eeVRZpcjeHtx1ajz8/Jd84PE+LArI17SZG271P9+2l7InIo0j/9gfHl688v3dnX58P3\nqfZXN6+7P2c3kBHp8+Hv819XFmlyNXa3n1gNe2Rn2PRqiPvCdXEvKdL+k5Ptm0zHPkfqPofs\n3tF96M77v+0+ZNlf3XyG8vfvcQ43kBHpoTsO+8O1aX/5//pyZZEmV2ObU6vxuX1l9Kt/S/zE\naoj7wvVwLyrS9rP6/t3l3XHbNLtvNmw+6+9PajYf2f/a3HV/9efubeuzN7C5+VCkX5+2XxZ4\n2J3kXVOkTNmTq/H37pEnVkPcF66He1GR2t+/3z5tv112RKT2AGkPri/fv/VvQnXfLvvy7/au\nydWfX9tj6JdnA5sbD0X6/fPblthfuLZImbK/T65Ge9NnexY6tRrivnDP515YpDOb/Dr5AeNs\nrCjX5pauRpR1qIp7GyJpEmKHwK2Ti0hw4Qq4iAQXroCLSHDhCriIBBeugItIcOEKuIgEF66A\ni0hw4Qq4iAQXroCLSHDhCriIBBeugItIcOEKuIgEF66Ai0hw4Qq4iAQXroCLSHDhCriIBBeu\ngItIcOEKuIgEF66Ai0hw4Qq4iAQXroCLSHDhCriIBBeugItIcOEKuIgEF66AqxHJmb+utylF\nKqtL3wsn0xeRplJZXfpeOIg0M5XVpe+Fg0gzU1ld+l44tyMSIXGDSIQIgkiECIJIhAiCSIQI\ngkiECHI1kRZtrrUtQRZV9bWm9VTe962h8LbmybrXEmmx+6OOVFS128P9n7s/bj2b47GKrvt1\nPb2+iHQ8FVVtu9Yl0mKNSMVtKkk9TfvUJVJVVTdBpHmp5fR9kzpFqmmNEWleqq1bR+fa+vY1\nEWlu6qlb24G5OHLppoNIJamnLiJdNmesLyIdTbV16+hcWd/kTUZE8mX3VF5HKjsw91WrqLtI\n/vpwkSr62L1PXXXr/WbDBxc5J4vt24s38c0GQkIHkQgRBJEIEQSRCBEEkQgRBJEIEQSRCBEE\nkQgRBJEIEQSR5Hm7a5plc87CPnUflJ9zz+eH9NpLu4H7l3X/2C5zSk7n9a5ZrN5T+sOzdgsh\ng0jyLM4+vPs7nXHPl4Nvprwb/7U95i8gkjEXbyl98SLdRMggkjzn6nHmndosntJrT82qfZZo\nHtsnqu6SOA/N02YTe/pTDV+J++Agkjr2a7x35O2haR7e2iOyOzN7abozpPvm+eCeZlP7v/vm\nfv22bO67k6r37nHvu/ut+uO4PeW67+98b09G993x/uPI5idQL/ftE82qv8fbvV1KsqOb3M3d\nAX14bzIKIqmzF+m9P8lbvK+bzoRVY0fx4T23IrVHefOjfWXVO9c/brm923v/wB626u68sEO9\ne/HSPN83y9fB5o+jnu1EbdXdY7G5tM+e3hyjr5r3NTkZRJJnd2q36n6v37VH7EP3HLLodHpt\nHo7dszvof3TH9o/u2mN3kK+a7fncY/8k1sHe7w4P9XuT4+1w68dRy+755dUef/fePt0cnK7t\n6XfdnVcD+nN3IklOBZHk2emx7A7Ct/b5oDsOX9pf+C/t8ft87J794dp0v/btcf2N95t73ffH\ncg97PRSp6WhPB26eQL09P97tNjZ4fbanv3Ty3A3ob7syZCKIJE+ix/Zq+/u+PTtqnx3um+l7\nJtcO3o3bY4bY9PYRc4i6214aPn6wjefu1dWQrn6PPV5YIHmOiPTQvC/u1/eL9eFv9iKRFgeU\nMXOAemiWT89vZ4jU5a07KT16C5kICyTP+NSuO7drX6L8aM/ufhy95+HRv2zGvC2sv3n7rt3m\nbYflUeYA1Y/ep0Ta0/u8dC+tUjoi5cICyTN+s6F/ZnjrP0l9P3rPw6O/f3/vx+5JwV4jPXZv\nEdzZzQ+bj3s2n/k8HWUOUE37Cm34ZsU+e/qyfWX03r/DntB5jZQNIsmzO4p3b3/3J1br7iC9\nG9xzcUwke1yzfVvb3rXr38juD/XNNw5a7Nv+ffLDE8EjqNXuJC8RaW/Tjv60fZc8pfOuXTaI\nJM/+EN18ILvefkvgcfAB6tNxkfrH3e2+lmOfI63f7ps7O/nqvmtnN7ezZvu9uNHWB6j+wqRI\ne/rTsn01tT6k8zlSNoh0+1ntd9Lka5W3Wd/iGTxqkn6BryJFCyJVkP137SYP9eWs75UOHjVF\n57t2+SDStdPscvZDXibe6lZnis63v/NBpGtnhkj7f4/0ISLx75HOCCIRIggiESIIIhEiCCIR\nIggiESIIIhEiCCIRIggiESLI/wP5D8U4yBvIfgAAAABJRU5ErkJggg==",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "flowmore <- moreplot(patientmonitor, \"flow_time\")\n",
    "flowmore"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this MORE plot, the solid lines are the mean and 0.05 and 0.95 quantiles. This provides a measure of the risk, what is the range of values that will generally be experienced by entities in this system (visitors to a hospital going to the information desk).  The error bars below the x axis provide measures of the error, which are the confidence intervals of the estimate of the mean, the 0.05 quantile, and the 0.95 quantile.  Note that the measure of error is not the same for each of the three.  Depending on what performance measure is most useful for the system, the analyst would then decide how many replications would be needed to reduce that particular performance measure's C.I. to an acceptable level."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}