\documentclass[11pt]{article} \usepackage{indentfirst} \usepackage{natbib} \usepackage{url} \begin{document} \title{Fuzzy Error Bars} \author{Charles J. Geyer} \maketitle \section{Introduction} This document tries out fuzzy error bars'' using the R package \verb@ump@. The idea was mentioned in \citet{GeyerMeeden05}, but was not adequately explored due to lack of space. \section{Example 1: Keeping It Simple} For our first example, we make up some data. <>= set.seed(42) fred <- c("red", "orange", "yellow", "green", "blue", "indigo", "violet", "cyan", "magenta", "white", "black") n <- seq(15, by = 5, length = length(fred)) x <- rbinom(length(n), n, 0.5) @ and calculate the fuzzy confidence intervals <>= library(ump) p <- seq(0, 1, 0.0001) alpha <- 0.05 phi <- matrix(NaN, length(p), length(fred)) for (i in 1:ncol(phi)) phi[ , i] <- umpu.binom(x[i], n[i], p, alpha) dimnames(phi) <- list(as.character(p), fred) @ While we are at it, we calculate the core <>= make.core <- function(phi) { nam <- dimnames(phi)[[2]] p <- as.numeric(dimnames(phi)[[1]]) stopifnot(all(0 <= p & p <= 1)) foo <- matrix(NaN, ncol(phi), 2) for (i in 1:ncol(phi)) { phii <- phi[, i] foo[i, 1] <- min(p[1 - phii == 1]) foo[i, 2] <- max(p[1 - phii == 1]) } dimnames(foo) <- list(nam, c("low", "high")) return(foo) } make.core(phi) @ and support <>= make.support <- function(phi) { nam <- dimnames(phi)[[2]] p <- as.numeric(dimnames(phi)[[1]]) stopifnot(all(0 <= p & p <= 1)) foo <- matrix(NaN, ncol(phi), 2) for (i in 1:ncol(phi)) { phii <- phi[, i] foo[i, 1] <- min(p[1 - phii > 0]) foo[i, 2] <- max(p[1 - phii > 0]) } dimnames(foo) <- list(nam, c("low", "high")) return(foo) } make.support(phi) @ Figure~\ref{fig:one} is produced by the following code <>= par(mar = c(5, 6, 1, 1) + 0.1) delta <- 0.3 foo <- make.support(phi) plot(range(foo), c(0, length(fred) - delta), type = "n", axes = FALSE, ylab = "", xlab = "") axis(side = 1) box() for (i in 1:ncol(phi)) { phii <- phi[ , i] ppp <- p[1 - phii > 0] xxx <- c(ppp, rev(range(ppp))) yyy <- (1 - phii)[1 - phii > 0] yyy <- c(yyy, 0, 0) yyy <- i - 1 + yyy * (1 - delta) polygon(xxx, yyy) lines(rep(x[i] / n[i], 2), c(i - 1, max(yyy))) } axis(side = 2, at = seq(along = fred) - 1 + (1 - delta) / 2, lab = fred, las = 1, padj = 1) title(xlab = expression(p)) @ and appears on p.~\pageref{fig:one}. \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Multiple Fuzzy Error Bars. Fuzzy confidence intervals associated with the UMPU test for the binomial distribution. Sample sizes range between \Sexpr{min(n)} and \Sexpr{max(n)}. Height of the error bar denotes the amount of fuzziness (full height, no fuzziness). Line in the middle of each bar denotes the MLE.} \label{fig:one} \end{figure} \section{Example 2: Pointy Haired Bars (PHB)} The associate editor handling \citet{GeyerMeeden05} quite rightly pointed out that it is not \emph{always} the case that our fuzzy error bars have full height, so the implied vertical axis for a fuzzy error bar goes from zero at the bottom to one at the top. However these special cases are self-labeling by their special form. We call them \emph{pointy haired bars} (PHB), a Dilbert allusion. To see how this works, let's do an example of the kind of confidence interval explanation seen in many introductory statistics textbooks where a single plot shows the result of a small simulation study, and (sure enough) about 95\% of the 95\% confidence intervals for simulated data cover the simulation truth'' parameter value and about 5\% don't cover. Intro stats books usually do this with normal data (or at least quasicontinuous data). We'll use binomial. Also (to be fair to the intro stats books) we obey the more than five expected in each cell'' rule for contingency tables with $n = 100$ and $p_0 = 0.05$ for the simulation truth. Thus we are in large $n$'' territory according to the time honored (but dumb) rule. Some newer textbooks, such as \citet{WildSeber00} do strongly criticize the time honored rule of thumb. Let's do it. <>= n <- 100 p0 <- 0.05 nsim <- 50 x <- rbinom(nsim, n, p0) print(x) @ O.~K. We got a zero, so we'll see the problem which is the ostensible point of this section. Figure~\ref{fig:two} is produced by the code very similar to Figure~\ref{fig:one}, so we will not show the actual code. <>= p <- seq(0, 1, 0.0001) alpha <- 0.10 phi <- matrix(NaN, length(p), length(x)) for (i in 1:ncol(phi)) phi[ , i] <- umpu.binom(x[i], n, p, alpha) dimnames(phi) <- list(as.character(p), NULL) @ <>= par(mar = c(5, 1, 1, 1) + 0.1) delta <- 0.3 foo <- make.support(phi) plot(range(foo), c(0, length(x) - delta), type = "n", axes = FALSE, ylab = "", xlab = "") axis(side = 1) box() for (i in 1:ncol(phi)) { phii <- phi[ , i] ppp <- p[1 - phii > 0] xxx <- c(ppp, rev(range(ppp))) yyy <- (1 - phii)[1 - phii > 0] yyy <- c(yyy, 0, 0) yyy <- i - 1 + yyy * (1 - delta) polygon(xxx, yyy) lines(rep(x[i] / n, 2), c(i - 1, max(yyy))) } title(xlab = expression(p)) abline(v = p0, lty = 2) @ The figure appears on p.~\pageref{fig:two}. \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Multiple Fuzzy Error Bars. Fuzzy confidence intervals associated with the UMPU test for the binomial distribution. Height of the error bar denotes the amount of fuzziness (full height, no fuzziness). All bars with flat top have full height. The bar without flat top (number \Sexpr{seq(along = x)[x == 0]} counting from bottom) only goes up to the confidence level, which is \Sexpr{1 - alpha}. Line in the middle of each bar denotes the MLE. Vertical dashed line is the simulation truth'' parameter value \Sexpr{p0}. All sample sizes $n = \Sexpr{n}$.} \label{fig:two} \end{figure} Hmmmmm. A bit cluttered with \Sexpr{nsim} intervals, but we'll leave it. The fuzziness can still be seen. Oh. Almost forgot. How well did we do? <>= coverage <- 1 - phi[p == 0.05, ] sum(coverage == 1) sum(coverage == 0) sort(coverage[0 < coverage & coverage < 1]) mean(coverage) sd(coverage) / sqrt(nsim) @ So we fully covered \Sexpr{sum(coverage == 1)} times, completely missed \Sexpr{sum(coverage == 0)} times, and partially covered \Sexpr{sum(0 < coverage & coverage < 1)} times. The empirical expected coverage was \Sexpr{round(mean(coverage), 4)} as opposed to the nominal \Sexpr{1 - alpha} level of the intervals. Not bad for simulation sample size \Sexpr{nsim}, only a little more than one standard error (Monte Carlo standard error \Sexpr{round(sd(coverage) / sqrt(nsim), 4)}) away from the theoretically correct result. \section{Example 3: Technical Quibbles} It is tempting to describe the exceptional interval in Figure~\ref{fig:two} as abutting the end of the data space (zero), but that's not a good description. Because the interval two above it, \emph{also} satisfies that description, but is not problematic. Too see this clearly, here's Figure~\ref{fig:two} redone to just show the intervals for $x = 0$ or $x = 1$. <>= par(mar = c(5, 1, 1, 1) + 0.1) delta <- 0.3 phi <- phi[ , x == 0 | x == 1] x <- x[x == 0 | x == 1] foo <- make.support(phi) plot(range(foo), c(0, length(x) - delta), type = "n", axes = FALSE, ylab = "", xlab = "") axis(side = 1) box() for (i in 1:ncol(phi)) { phii <- phi[ , i] ppp <- p[1 - phii > 0] xxx <- c(ppp, rev(range(ppp))) yyy <- (1 - phii)[1 - phii > 0] yyy <- c(yyy, 0, 0) yyy <- i - 1 + yyy * (1 - delta) polygon(xxx, yyy) lines(rep(x[i] / n, 2), c(i - 1, max(yyy))) } title(xlab = expression(p)) abline(v = p0, lty = 2) @ This is Figure~\ref{fig:three}, which appears on p.~\pageref{fig:three}. \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Fuzzy Error Bars that Abut the End. Fuzzy confidence intervals from Figure~\ref{fig:two} that have zero for their left endpoint. The bar with flat top has full height. The bar with pointed top only goes up to the confidence level, which is \Sexpr{1 - alpha}.} \label{fig:three} \end{figure} Thus we see that the correct description of \emph{full height} fuzzy error bars is \emph{having flat top}, not \emph{abutting the end of the range}. There's another technical reason for this quibble as the following case shows. <>= n <- 5 x <- 2 alpha <- 2 / 3 phii <- umpu.binom(x, n, p, alpha) @ <>= par(mar = c(5, 1, 1, 1) + 0.1) delta <- 0.3 phi <- as.matrix(phii) dimnames(phi) <- list(as.character(p), NULL) foo <- make.support(phi) plot(range(foo), c(0, length(x) - delta), type = "n", axes = FALSE, ylab = "", xlab = "") axis(side = 1) box() for (i in 1:ncol(phi)) { phii <- phi[ , i] ppp <- p[1 - phii > 0] xxx <- c(ppp, rev(range(ppp))) yyy <- (1 - phii)[1 - phii > 0] yyy <- c(yyy, 0, 0) yyy <- i - 1 + yyy * (1 - delta) polygon(xxx, yyy) lines(rep(x[i] / n, 2), c(i - 1, max(yyy))) } title(xlab = expression(p)) @ This is Figure~\ref{fig:four}, which appears on p.~\pageref{fig:four}. \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Another Fuzzy Error Bar without Full Height. Fuzzy confidence interval associated with the UMPU test for the binomial distribution. Data $x = \Sexpr{x}$, $n = \Sexpr{n}$, confidence level \Sexpr{round(100 * (1 - alpha), 2)}\%. The value of the membership function of the fuzzy interval goes from zero at the bottom to \Sexpr{round(max(1 - phii), 3)} at the apex.} \label{fig:four} \end{figure} Admittedly a \Sexpr{round(100 * (1 - alpha), 2)}\% confidence interval has no conceivable application. We show it only to make the technical point that very low confidence level combined with small $n$ can give another case where the fuzzy error bar does not have full height. And it is also another case where \emph{full height} goes with \emph{flat top} and \emph{non-full height} or \emph{empty core} goes with \emph{pointed top}. For the curious, the reason for the pointed top is that in order to reject at such a high $\alpha$ we must reject, at least partially, at \emph{all} points, even the point $x$, even when $p = \hat{p}$, where the fuzzy confidence interval membership function always achieves its maximum (possibly one, possibly less than one). This latter claim is a theorem from unpublished notes of the author: the membership function of a fuzzy confidence interval associated with the UMPU test for a one-parameter discrete exponential family achieves its maximum when the parameter value equals the MLE. After that bit of theory, we see that abstractly the behavior in the bottom interval of Figure~\ref{fig:three} is just the same as the behavior of the interval in Figure~\ref{fig:four}. The only difference is that in Figure~\ref{fig:three}, the MLE of the mean value parameter is on the boundary of the sample space, and when that happens, the MLE distribution concentrates all probability at the point on the boundary, and thus we see the irregularity for all $\alpha$ levels. But when the phenomenon occurs in the interior of the sample sample space, we see the irregularity only for high $\alpha$ levels, hence low confidence levels $$1 - \alpha < f_{\hat{\mu}}(x)$$ where $\hat{\mu} = x$ is the MLE of the mean value parameter. See equations (3.7a) and (3.7b) and the following discussion in \citet{GeyerMeeden05} for explanation. \begin{thebibliography}{} \bibitem[{Geyer and Meeden(2005)}]{GeyerMeeden05} \textsc{Geyer, C.~J.} and \textsc{Meeden, G.~D.} (2005). \newblock Fuzzy and randomized confidence intervals and p-values. \newblock Revised and resubmitted (thrice) to \emph{Statistical Science}, \url{http://www.stat.umn.edu/geyer/fuzz}. \bibitem[{Wild and Seber(2000)}]{WildSeber00} \textsc{Wild, C.~J.} and \textsc{Seber, G.~A.~F.} (2000). \newblock \emph{Chance Encounters: A First Course in Data Analysis and Inference}. \newblock Wiley. \end{thebibliography} \end{document}