9

How can I calculate a 95% interval to estimate the actual proportion of SUV's in the city in R? I would like to calculate the interval on this data:

vehicleType <- c("suv", "suv", "minivan", "car", "suv", "suv", "car", "car", "car", "car", "minivan", "car", "truck", "car", "car", "car", "car", "car", "car", "car", "minivan", "car", "suv", "minivan", "car", "minivan", "suv", "suv", "suv", "car", "suv", "car", "car", "suv", "truck", "truck", "minivan", "suv", "car", "truck", "suv", "suv", "car", "car", "car", "car", "suv", "car", "car", "car", "suv", "car", "car", "car", "truck", "car", "car", "suv", "suv", "minivan", "suv", "car", "car", "car", "car", "car", "minivan", "suv", "car", "car", "suv", "minivan", "car", "car", "car", "minivan", "minivan", "minivan", "car", "truck", "car", "car", "car", "suv", "suv", "suv", "car", "suv", "suv", "car", "suv", "car", "minivan", "car", "car", "car", "car", "car", "car", "car")

Thanks in advance.

  • Besides speed, I'll bet that the non-coverage probability in at least one of the two tails is more accurate than the bootstrap. – Frank Harrell Apr 25 '17 at 03:12

3 Answers3

16

First, remember that an interval for a proportion is given by:

p_hat +/- z * sqrt(p_hat * (1-p_hat)/n)

With that being said, we can use R to solve the formula like so:

# Set CI alpha level (1-alpha/2)*100%
alpha = 0.05

# Load Data
vehicleType = c("suv", "suv", "minivan", "car", "suv", "suv", "car", "car", "car", "car", "minivan", "car", "truck", "car", "car", "car", "car", "car", "car", "car", "minivan", "car", "suv", "minivan", "car", "minivan", "suv", "suv", "suv", "car", "suv", "car", "car", "suv", "truck", "truck", "minivan", "suv", "car", "truck", "suv", "suv", "car", "car", "car", "car", "suv", "car", "car", "car", "suv", "car", "car", "car", "truck", "car", "car", "suv", "suv", "minivan", "suv", "car", "car", "car", "car", "car", "minivan", "suv", "car", "car", "suv", "minivan", "car", "car", "car", "minivan", "minivan", "minivan", "car", "truck", "car", "car", "car", "suv", "suv", "suv", "car", "suv", "suv", "car", "suv", "car", "minivan", "car", "car", "car", "car", "car", "car", "car")

# Convert from string to factor
vehicleType = factor(vehicleType)

# Find the number of obs
n = length(vehicleType)

# Find number of obs per type
vtbreakdown = table(vehicleType)

# Get the proportion
p_hat = vtbreakdown['suv']/n

# Calculate the critical z-score
z = qnorm(1-alpha/2)

# Compute the CI
p_hat + c(-1,1)*z*sqrt(p_hat*(1-p_hat)/n)

So, we have:

0.1740293 0.3459707

For the p_hat of:

0.26
coatless
  • 296
  • 3
  • 6
  • Should the "Compute the CI" code be `p_hat + c(-1.96,1.96)*sqrt(p_hat*(1-p_hat)/n)` for a 95% confidence interval (where `c(-1.96,1.96)` = -/+ the value of `z`)? – eipi10 Apr 17 '16 at 04:46
  • No need to do that @eipi10. Just forgot to include the `z`. @jason-todd note the correction needed. – coatless Apr 17 '16 at 04:48
  • Nice answer. It may be worth noting that the score interval tends to perform better in terms of coverage probabilities. – Josh Jan 22 '18 at 15:27
16

@Coatless's method will get the job done in most cases (including the OP's case). However, for completeness, I thought I'd add a couple of other options.

Bootstrap Method

The function below draws n resamples from the data vector. For each resample, it calculates the proportion of "successes" and then calculates the overall mean and 95% confidence interval

bp = function(x, lev, n = 1e3, alpha=0.05) {
  res = replicate(n, sum(sample(x, length(x), replace=TRUE) == lev)/length(x))
  return(list(mean=mean(res),
              `95% CI`=quantile(res, c(0.5*alpha,1-0.5*alpha))))
}

bp(vehicleType, "suv")

$mean
[1] 0.259628

$`95% CI`
  2.5% 97.5% 
  0.18  0.35

binom Package

The binom package will run the test in @Coatless's answer, which assumes the errors are normally distributed. This can result in incorrect values when the proportion of "successes" is near zero or one and/or if the sample is relatively small. binom.confint from the binom package has other options that avoid this pitfall.

In the output below, the asymptotic test is the same as the one coded by @Coatless. You can get the results for just one of the methods by using, for example, the methods="exact" argument. Also, binom.test() uses the exact (Pearson-Klopper) test by default.

library(binom)

binom.confint(sum(vehicleType=="suv"), length(vehicleType))

          method  x   n      mean     lower     upper
1  agresti-coull 26 100 0.2600000 0.1836007 0.3541561
2     asymptotic 26 100 0.2600000 0.1740293 0.3459707
3          bayes 26 100 0.2623762 0.1788095 0.3485750
4        cloglog 26 100 0.2600000 0.1787357 0.3485852
5          exact 26 100 0.2600000 0.1773944 0.3573121
6          logit 26 100 0.2600000 0.1835016 0.3545416
7         probit 26 100 0.2600000 0.1818365 0.3526030
8        profile 26 100 0.2600000 0.1808127 0.3513344
9            lrt 26 100 0.2600000 0.1808329 0.3513338
10     prop.test 26 100 0.2600000 0.1797427 0.3590222
11        wilson 26 100 0.2600000 0.1840470 0.3537099
eipi10
  • 637
  • 5
  • 16
  • 6
    I recommend Wilson's method which you can also get from the R `Hmisc` package `binconf` function. – Frank Harrell Sep 21 '16 at 22:20
  • 2
    @FrankHarrell, is there a reason to prefer Wilson's method over the bootstrap (or vice versa)? If the answer is too long for a comment, please feel free to link to something you or others have written about this. – eipi10 Apr 24 '17 at 21:06
0

Because we are using a continuous normal distribution as an approximation of a discrete binomial distribution there should be a correction term added (0.5/N) to the above calculations:

enter image description here

See here for more details

Kneef
  • 11
  • 3