Calculating Confidence Intervals for two datasets

Question

Quite the number of questions I've made today.

I'd like to calculate the Confidence Interval (99% level, not 95) for the mean value of variable age of two dataframes, infert_control and infert_patient where:

infert_control = subset(infert$age, infert$case == 0)
infert_patient = subset(infert$age, infert$case == 1) 

infert is a built-in R dataset, for those not familiar with it, here it is: case 0 denominates the control-group patients, case 1 the actual ones.

> infert
    education age parity induced case spontaneous stratum pooled.stratum
1      0-5yrs  26      6       1    1           2       1              3
2      0-5yrs  42      1       1    1           0       2              1
3      0-5yrs  39      6       2    1           0       3              4
4      0-5yrs  34      4       2    1           0       4              2
5     6-11yrs  35      3       1    1           1       5             32
6     6-11yrs  36      4       2    1           1       6             36
7     6-11yrs  23      1       0    1           0       7              6
8     6-11yrs  32      2       0    1           0       8             22
9     6-11yrs  21      1       0    1           1       9              5
10    6-11yrs  28      2       0    1           0      10             19
11    6-11yrs  29      2       1    1           0      11             20
...
239   12+ yrs  38      6       0    0           2      74             63
240   12+ yrs  26      2       1    0           1      75             49
241   12+ yrs  31      1       1    0           0      76             45
242   12+ yrs  31      2       0    0           1      77             53
243   12+ yrs  25      1       0    0           1      78             41
244   12+ yrs  31      1       0    0           1      79             45
245   12+ yrs  34      1       0    0           0      80             47
246   12+ yrs  35      2       2    0           0      81             54
247   12+ yrs  29      1       0    0           1      82             43
248   12+ yrs  23      1       0    0           1      83             40

What would be the correct way to solve this?

I've already calculated the Mean value of column age for both infert_control and infert_patient, plus the standard deviation of each subset.

Answer 1

You could easily calculate the confidence interval manually:

infert_control <- subset(infert$age, infert$case == 0)

# calculate needed values
m <- mean(infert_control)
s <- sd(infert_control)
n <- length(infert_control)

# calculate error for normal distribution (choose you distribution here, e.g. qt for t-distribution)
a <- 0.995 # 99% CI => 0.5% on both sides
error <- qnorm(a)*s/sqrt(n)

# calculate CI
ci_lower <- m-error
ci_upper <- m+error

See also http://en.wikipedia.org/wiki/Confidence_interval (sorry for a wikipedia link, but it has a good explanation and shows you the formula)

Answer 2

You could use bootstrap for this:

library(boot)
set.seed(42)
boot_mean <- boot(infert_control, function(x, i) mean(x[i]), R=1e4)
quantile(boot_mean$t, probs=c(0.005, 0.995))
#      0.5%    99.5% 
#  30.47273 32.58182 

Or if you don't want to use a library:

set.seed(42)
R <- 1e4
boot_mean <- colMeans(
                matrix(
                   sample(infert_control, R * length(infert_control), TRUE), 
                   ncol=R))
quantile(boot_mean, probs=c(0.005, 0.995))
#    0.5%    99.5% 
#30.42424 32.55152 

Answer 3

So many answers...

The mean value of a random sample has a t-distribution, not normal, although t -> N as df -> Inf.

cl <- function(data,p) {
  n  <- length(data)
  cl <- qt(p/2,n-1,lower.tail=F)*sd(data)/sqrt(n)
  m  <- mean(data)
  return(c(lower=m-cl,upper=m+cl))
}
cl.control <- cl(infert_control,0.01)
cl.control
#    lower    upper 
# 30.42493 32.55689 

cl.patient <- cl(infert_patient,0.01)
cl.patient
#    lower    upper 
# 30.00221 33.05803

aggregate(age~case,data=infert,cl,p=0.01)  # much better way...
#   case age.lower age.upper
# 1    0  30.42493  32.55689
# 2    1  30.00221  33.05803

Also, the quantile functions (e.q. qt(...) and qnorm(...)) return the lower tail by default, so your limits would be reversed unless you set lower.tail=F

Answer 4

... or as small function:

cifun <- function(data, ALPHA){
  c(mean(data) - qnorm(1-ALPHA/2) * sd(data)/sqrt(length(data)),
    mean(data) + qnorm(1-ALPHA/2) * sd(data)/sqrt(length(data)))
}

cifun(infert_control, 0.01)