# How do I calculate single variable confidence intervals assuming different distributions with R?

I'm trying to convert a bunch of Stata commands to their R equivalents, and I'm struggling with how R does handle confidence intervals for inferential statistics of single variables.

In Stata, I can use ci variable to calculate normal confidence intervals, ci variable, b to calculate binomial intervals, and ci variable, p to calculate the intervals assuming the variable is distributed as Poisson.

R seems slightly more complicated, though (confint() only works with models…). It appears that I have to run the respective tests to get the confidence intervals (t.test, binom.test, and poisson.test). I'm fine doing that, but these functions seem to be less generalized than Stata's ci command.

For example, to calculate a normal 95% confidence interval I use t.test:

set.seed(1234)
x <- rnorm(100)
t.test(x)

# 95 percent confidence interval:
#  -0.35605755  0.04253406

However, it seems to be a lot more complicated to calcuate binomial or Poisson intervals. For example, if I have a column of binary data (yes/no; 1/0) like x below, Stata appears to convert it into count data automatically when running ci x, b. In R, I have to convert it to count data on my own with sum():

set.seed(1234)
x <- sample(0:1, 100, replace=TRUE)  # lots of 0s and 1s
binom.test(sum(x), length(x))

# 95 percent confidence interval:
# 0.3503202 0.5527198

Is manually feeding in the sum and the length of the variable the official R-esque way of calculating the confidence interval assuming a binomial distribution, or is there a better way?

Likewise, what's the most R-esque way to calculate confidence intervals for a variable assumed to follow a Poisson distribution--the equivalent of ci x, p:

set.seed(1234)
x <- rpois(100, 5)  # random Poisson distribution
# ... magic R voodoo ...
# Confidence interval!

So, I guess in summary, what's the best R equivalent for Stata's ci x, ci x, b, and ci x, p commands for single variables?

-
I really prefer the R framework. There are many ways to construct a confidence interval for a Gaussian mean, a confidence interval for a binomial proportion, a confidence interval for a Poisson rate... For example the binom.test() function provides the Clopper-Pearson binomial test with its associated confidence interval, as shown by ?binom.test. And there is no reason to impose a "default" confidence interval to statisticians who want to understand what they are doing. –  Stéphane Laurent Mar 3 '12 at 14:45
add comment

## 1 Answer

One way is to do a regression of your variable on just an intercept, i.e.

x <- rpois(100, 5)
glm1 <- glm(x~1, family=poisson)
confint(glm1) # or confint.default(), if you want less-clever intervals

This gives a confidence interval for the log of the mean of the X's. If you want a confidence interval around something else, either use the link argument in your family= statement, or transform the interval provided by the code above, e.g.

exp(confint(glm1))

As Stéphane notes, there's also the poisson.test() function. Its syntax is a little removed from typical regression functions, and for this example works out to be;

poisson.test(sum(x), 100)

In large samples there will be little difference between the methods given.

-
We could prefer poisson.test() which provides an "exact" confidence interval (similarly to binom.test() for a binomial proportion) –  Stéphane Laurent Mar 3 '12 at 14:48
So how would you use poisson.test() with the x given in this answer (rpois(100, 5)? This seems overly complicated and wrong: poisson.test(as.integer(mean(x))) –  Andrew Mar 3 '12 at 18:35
add comment