[R-gregmisc-commits] r2058 - in pkg/gmodels: R man

[noreply at r-forge.r-project.org]

Author: warnes
Date: 2015-07-19 04:30:05 +0200 (Sun, 19 Jul 2015)
New Revision: 2058

Modified:
   pkg/gmodels/R/ci.R
   pkg/gmodels/man/ci.Rd
Log:
ci.binom() was using an incorrect method for calculating binomial confidence interval.   The revised code calculates the Clopper-Pearson 'exect' interval, which is *conservative* due to the discrete nature of the binomial distribution.

Modified: pkg/gmodels/R/ci.R
===================================================================
--- pkg/gmodels/R/ci.R	2015-07-14 21:19:03 UTC (rev 2057)
+++ pkg/gmodels/R/ci.R	2015-07-19 02:30:05 UTC (rev 2058)
@@ -15,20 +15,24 @@
                  "CI upper"=ci.high,
                  "Std. Error"=stderr
                  )
-    
+
     retval
   }
 
 ci.binom <- function(x, confidence=0.95,alpha=1-confidence,...)
   {
-    if( !(all(x) %in% c(0,1)) ) stop("Binomial values must be either 0 or 1.")
+      if( !(all(x %in% c(0,1))) ) stop("Binomial values must be either 0 or 1.")
+      if( all(x==0) || all(x==1) )
+          warning("All observed values are ", as.numeric(x[1]), ", so estimated Std. Error is 0.")
 
     est  <-  mean(x, na.rm=TRUE)
     n <- nobs(x)
+    x <- sum(x)
     stderr <- sqrt(est*(1-est)/n)
-    ci.low  <- qbinom(p=alpha/2, prob=est, size=n)/n
-    ci.high <- qbinom(p=1-alpha/2, prob=est, size=n)/n
 
+    ci.low  <- qbeta(   alpha/2, x  , n + 1 - x)
+    ci.high <- qbeta(1- alpha/2, x+1, n-x      )
+
     retval  <- cbind(Estimate=est,
                      "CI lower"=ci.low,
                      "CI upper"=ci.high,
@@ -48,7 +52,7 @@
                      "CI upper"=ci.high,
                      "Std. Error"= coef(x)[,2],
                      "p-value" = coef(x)[,4])
-    
+
     retval
   }
 
@@ -68,18 +72,18 @@
     retval
   }
 
-ci.mer <- function (x, 
-                    confidence = 0.95, 
-                    alpha = 1 - confidence, 
+ci.mer <- function (x,
+                    confidence = 0.95,
+                    alpha = 1 - confidence,
                     n.sim = 1e4,
-                    ...) 
+                    ...)
 {
     x.effects <- x at fixef
     n <- length(x.effects)
 
-    retval <- gmodels:::est.mer(obj = x, 
+    retval <- gmodels:::est.mer(obj = x,
                                 cm = diag(n),
-                                beta0 = rep(0, n), 
+                                beta0 = rep(0, n),
                                 conf.int = confidence,
                                 show.beta0 = FALSE,
                                 n.sim = n.sim)
@@ -90,7 +94,7 @@
                      ]
     colnames(retval)[c(2:3, 5)] <- c("CI lower", "CI upper", "p-value")
     rownames(retval) <- names(x.effects)
-    
+
     retval
 }
 
@@ -106,6 +110,6 @@
                      "p-value" = x$"Pr(>|t|)"
                      )
     rownames(retval) <- rownames(x)
-    
+
     retval
   }

Modified: pkg/gmodels/man/ci.Rd
===================================================================
--- pkg/gmodels/man/ci.Rd	2015-07-14 21:19:03 UTC (rev 2057)
+++ pkg/gmodels/man/ci.Rd	2015-07-19 02:30:05 UTC (rev 2058)
@@ -20,8 +20,8 @@
   \method{ci}{binom}(x, confidence=0.95, alpha=1-confidence, ...)
   \method{ci}{lm}(x, confidence=0.95, alpha=1-confidence, ...)
   \method{ci}{lme}(x, confidence=0.95, alpha=1-confidence, ...)
-  \method{ci}{mer}(x, confidence=0.95, alpha=1-confidence, n.sim=10000, ...) 
-  \method{ci}{estimable}(x, confidence=0.95, alpha=1-confidence, ...) 
+  \method{ci}{mer}(x, confidence=0.95, alpha=1-confidence, n.sim=10000, ...)
+  \method{ci}{estimable}(x, confidence=0.95, alpha=1-confidence, ...)
 }
 \arguments{
   \item{x}{ object from which to compute confidence intervals. }
@@ -32,9 +32,12 @@
   \item{\dots}{Arguments for methods}
   \item{n.sim}{Number of samples to take in \code{mcmcsamp}.}
 }
-%\details{
-%  ~~ If necessary, more details than the __description__  above ~~
-%}
+\details{
+  \code{ci.binom} computes binomial confidence intervals using the
+  Clopper-Pearson 'exact' method based on the binomial quantile
+  function.  Due to the discrete nature of the binomial distribution,
+  this interval is conservative.
+}
 \value{
   vector or matrix with one row per model parameter and elements/columns
   \code{Estimate}, \code{CI lower}, \code{CI upper}, \code{Std. Error},
@@ -50,7 +53,7 @@
 
 \examples{
 
-# mean and confidence interval 
+# mean and confidence interval
 ci( rnorm(10) )
 
 # binomial proportion and exact confidence interval
@@ -62,7 +65,7 @@
 # confidence intervals for regression parameteres
 data(state)
 reg  <-  lm(Area ~ Population, data=as.data.frame(state.x77))
-ci(reg) 
+ci(reg)
 
 %\dontrun{
 # mer example