Confidence Intervals Around a Mean: biased (non-centered) confidence interval? (an exercise using R)

I've been playing around with the module "Confidence Intervals Around a Mean" (meanCI) from statsTeachR (www.statsteachr.org), authored by Eric A Cohen (unfortunately, author's contact information was not available).

It "focuses on understanding and calculating confidence intervals around a sample mean. This includes understanding the concept of a confidence interval; understanding when using either the normal or the t-distribution (or neither) is appropriate; writing R code to calculate confidence intervals; and writing code to test coverage of the true population mean by confidence intervals" (www.statsteachr.org/modules/19).

Building on this module, I produced this plot:

It shows the average number of "mistakes" (confidence intervals that do not contain the population mean) as the sample size (n) increases from 2 to 1000 (the population size N), distinguishing between "below mistakes" and "above mistakes" (above mistake: the lower limit of the CE is above the population mean; below mistake: the upper limit of the CE is below the population mean).

The figure has 4 panels: Top left: When the population follows a chi square distribution, and confidence interval using the normal distribution. Top right: When the population follows a chi square distribution, and confidence interval using the t distribution. Bottom left: When the population is normally distributed, and confidence interval using the normal distribution. Bottom right: When the population is normally distributed, and confidence interval using the t distribution.

My question is, why there seems to be some sort of bias, with below mistakes consistently higher than the above mistakes when the population follows a chi squared distribution and above mistakes consistently higher than below mistakes when the population is normally distributed?

Is it really like that?, I mean, is this a standard result?, or it is just an error in my code? (I swear I've checked my code many many times, and I think everything is ok). If there is no error in my code, doesn't it mean that we could improve confidence interval's performance trying to correct for such 'bias'?

I think I've been doing my homework and I cannot reconcile what I find in my exercise using R with the theory. Checking out my introductory statistics book, it doesn't say anything about non-centered confidence interval. In fact, it says that the mean estimator is consistent and thus, it gets closer to the true population mean as the sample size increases, and given that the confidence interval has the mean estimator in its center, it seems to me that the above and below mistakes should be relatively similar.

I also searched here in CrossValidated and found this VERY interesting post-answers (What, precisely, is a confidence interval?), but nobody there discussed anything about whether the confidence interval is indeed centered or not.

Any thought is more than welcome.

Here's my code:

There are three core functions as follows. The first two (CIs_normal_dist, CIs_t_dist) calculate "by hand" the confidence interval using the normal approximation and the t distribution respectively. Note that both functions assume the population variance is unknown and thus, estimate the standard error using the sample values (as discussed here Confidence intervals for mean, when variance is unknown and here http://www.r-tutor.com/elementary-statistics/interval-estimation/interval-estimate-population-mean-unknown-variance). The third function produces sets of CIs and calculates how many mistakes (population mean not in the CI) there are in a set, and repeat that many times to calculate the average number of mistakes.

# BASED ON:
#   statsTeachR module: Confidence Intervals around a Mean
#   Author: Eric A. Cohen

# Following statsTeachR module, this function takes many (num_CIs) samples
# from the data and calculates the sample mean and the confindence interval
# around it, using the normal approximation.
CIs_normal_dist <- function(data, num_CIs=100, alpha=0.1, sample_size=40,
replace_sampling=FALSE){
CIs <- matrix(data=NA, nrow=num_CIs, ncol=3)
for (i in 1:num_CIs) {
sample_i <- sample(x=data, size=sample_size, replace=replace_sampling)
CIs[i, 1] <- mean(sample_i)
#   Calculate the sample standard deviations
sample_sd <- sd(sample_i)
sample_mean_sd <- sample_sd/sqrt(sample_size)
#   Use the normal distribution to calculate the confidence interval:
delta <- qnorm(p=1-alpha/2, mean=0, sd=1) * sample_mean_sd
CIs[i, 2] <- CIs[i, 1] - delta
CIs[i, 3] <- CIs[i, 1] + delta
#just to check that the CI generated by our function
#is the same as the z.test function in the BSDA package
#print(c(CIs[i, 2], CIs[i, 3]))
#if(sample_size>3){
#print(z.test(x=sample_i, conf.level=1-alpha, sigma.x=sample_sd)$conf.int) #} } CIs } # Following statsTeachR module, this function takes many (num_CIs) samples # from the data and calculates the sample mean and the confindence interval # around it, using t distribution CIs_t_dist <- function(data, num_CIs=100, alpha=0.1, sample_size=40, replace_sampling=FALSE){ CIs <- matrix(data=NA, nrow=num_CIs, ncol=3) for (i in 1:num_CIs) { sample_i <- sample(x=data, size=sample_size, replace=replace_sampling) CIs[i, 1] <- mean(sample_i) # Calculate the sample standard deviations sample_sd <- sd(sample_i) sample_mean_sd <- sample_sd/sqrt(sample_size) # Use the t distribution to calculate the confidence interval: delta <- qt(p=1-alpha/2, df=sample_size-1) * sample_mean_sd CIs[i, 2] <- CIs[i, 1] - delta CIs[i, 3] <- CIs[i, 1] + delta #just to check that the CI generated by our function #is the same as the built-in t.test function #print(c(CIs[i, 2], CIs[i, 3])) #print(t.test(x=sample_i, conf.level=1-alpha)$conf.int)
}
CIs
}

# Use the previous function (CIs_normal_dist) to produce sets of CIs
# Calculate how many mistakes (population mean not in the CI) there are
# in a set of num_CIs CIs (distinguish mistakes above and below mistakes,
# above: the lower limit of the CE is above the population mean
# below: the upper limit of the CE is below the population mean)
# Do that num_CI_sets times, for each sample size between 1 and max_sample_size
# And take the average number of mistakes for each sample size
# The function returns a data frame with 4 columns
# mean_mistakes[,1] sample size
# mean_mistakes[,2] average number of mistakes above
# mean_mistakes[,3] average number of mistakes below
# mean_mistakes[,4] average number of total mistakes
# (averages of mistakes in num_CI_sets CIs)
mistakes_mean <- function(data, population_mean,
max_sample_size = 2, num_CI_sets = 100,
num_CIs = 100, alpha = 0.1, replace_sampling = FALSE,
FUN = CIs_t_dist){
mean_mistakes <- matrix(data=NA, nrow=max_sample_size-1, ncol=4)
for(sample_size in seq(from=2, to=max_sample_size, by=1)){
mistakes <- matrix(data=NA, nrow=num_CI_sets, ncol=4)
for(i in 1:num_CI_sets){
CIs <- FUN(data = data, num_CIs = num_CIs, alpha = alpha,
sample_size = sample_size, replace_sampling=FALSE)
mistakes[i,1] <- sample_size
mistakes[i,2] <- sum(CIs[,2]>population_mean) #above mistake
mistakes[i,3] <- sum(CIs[,3]<population_mean) #below mistake
}
mistakes[,4] <- mistakes[,2]+mistakes[,3]
mean_mistakes[sample_size-1,1] <- sample_size
mean_mistakes[sample_size-1,2] <- mean(mistakes[,2])
mean_mistakes[sample_size-1,3] <- mean(mistakes[,3])
mean_mistakes[sample_size-1,4] <- mean(mistakes[,4])
}
df <- as.data.frame(mean_mistakes)
names(df) <- c("Sample size","Above mistakes", "Below mistakes", "Total mistakes")
df
}


Now, I just use those functions for the exercise and summarize the results in a plot:

rm(list = ls(all = TRUE))
source("CE_bias_func.R") #I saved the functions above in this file
library(BSDA)
library(Hmisc)

# set the parameters for the exercise
pop_mean <- 2
alpha <- 0.05
N <- 1000

# Let's generate observations (our finite universe) chisq distributed
universe <- rchisq(n=N, df=pop_mean)
hist(universe)

mm1 <- mistakes_mean(data = universe, population_mean = pop_mean,
max_sample_size = N,
num_CI_sets = 100, num_CIs = 100,
alpha = 0.05, replace_sampling = FALSE,
FUN = CIs_t_dist)

mm2 <- mistakes_mean(data = universe, population_mean = pop_mean,
max_sample_size = N,
num_CI_sets = 100, num_CIs = 100,
alpha = 0.05, replace_sampling = FALSE,
FUN = CIs_normal_dist)

# generate observations (our finite universe), now, normally  distributed
pop_mean <- -15
universe <- rnorm(n=N, mean=pop_mean, sd=10)

mm3 <- mistakes_mean(data = universe, population_mean = pop_mean,
max_sample_size = N,
num_CI_sets = 100, num_CIs = 100,
alpha = 0.05, replace_sampling = FALSE,
FUN = CIs_t_dist)

mm4 <- mistakes_mean(data = universe, population_mean = pop_mean,
max_sample_size = N,
num_CI_sets = 100, num_CIs = 100,
alpha = 0.05, replace_sampling = FALSE,
FUN = CIs_normal_dist)

# Describe each simulation
mm1$"Description" <- "Population: chisq dist, Confidence Interval: t dist" mm2$"Description" <- "Population: chisq dist, Confidence Interval: normal dist"
mm3$"Description" <- "Population: normal dist, Confidence Interval: t dist" mm4$"Description" <- "Population: normal dist, Confidence Interval: normal dist"

#
mm <- rbind(mm1, mm2, mm3, mm4)
mm$"Total mistakes" <- NULL library(ggplot2) # to plot library(reshape2) # to reshape the data as ggplot2 like dfm <- melt(mm, id.vars=c("Sample size", "Description"), measure.vars=c("Below mistakes", "Above mistakes")) names(dfm) <- c("Sample_size", "Description", "variable", "value") gp <- ggplot(data=dfm, aes(x=Sample_size, y=value, color=variable)) gp <- gp + geom_point() gp <- gp + labs(x="Sample size", y="Avg. # of mistakes in 100 trials of 100 CIs", title="Average numer of mistakes in 100 trials of 100 confidence intervals each") gp <- gp + theme(legend.title=element_blank()) #turn off legen title gp <- gp + theme(legend.key=element_rect(fill=NA)) #get rid off legend's boxes gp <- gp + theme(legend.position="top") gp <- gp + facet_wrap(~Description, nrow=2, ncol=2) gp  UPDATE: After very useful comments from @swmo and @heropup I re-run the exercise, now taking samples from the whole distribution instead of samples from a finite population, so as to have truly independent samples. And this is the result: pretty much what you would expect in the first place. i) no considerable differences in above and below mistakes and ii) total number of CIs that do not include the true parameter (above+below mistakes) close to 5%. Here's my code (for now, I only did it for a normal dist and using a t.test, but very easily could be adapted to cover the 4 cases in my original post): # A function that takes 100 samples from the distribution, # given the sample size (n) received as argument # then calculate the CI for each of them # and calculate the number of mistakes (above and below) for each of them # It returns a data.frame with one row and three variables for sample size, # above mistakes and below mistakes trial <- function(n){ # Take 100 samples, of sample size n, from the whole distribution normal dist # and not from a finite population samples <- lapply(rep(n, 100), FUN=rnorm, mean=2, sd=1) # For each sample, calculate the confidence interval CIs <- lapply(samples, FUN = function(x) t.test(x=x, conf.level=0.95)$conf.int)

# For each confidence interval, determine above and below mistakes
above_mistakes <- sapply(CIs, FUN = function(x) if (x[1] > 2) 1 else 0)
below_mistakes <- sapply(CIs, FUN = function(x) if (x[2] < 2) 1 else 0)

data.frame(sample_size = n,
above_mistakes = sum(above_mistakes),
below_mistakes = sum(below_mistakes))
}

# Now create a function that replicates the trial 100 times
# for a given sample size (n) and aggregates the results (rbind) in a data.frame
trialx100 <- function(n) do.call("rbind", replicate(trial(n), n=100, simplify=FALSE))

# Finally, for sample sizes between 10 and 1000, use the trailx100 function.
sample_sizes <- seq(from = 10, to = 1000, by = 10)
mistakes_all <- do.call("rbind", sapply(sample_sizes, FUN = trialx100, simplify=FALSE))

# And take the mean of the mistakes
mistakes_agg <- aggregate(mistakes_all, by = list(mistakes_all$sample_size), FUN = mean) library(ggplot2) # to plot library(reshape2) # to reshape the data as ggplot2 like # Melt the data as ggplot fancy dfm <- melt(mistakes_agg, id.vars=c("sample_size"), measure.vars=c("above_mistakes", "below_mistakes")) gp <- ggplot(data=dfm, aes(x=sample_size, y=value, color=variable)) gp <- gp + geom_line(fill="black") gp <- gp + labs(x="Sample size", y="Avg. # of mistakes in 100 trials of 100 CIs", title="Average numer of mistakes in 100 trials of 100 confidence intervals each") gp <- gp + theme(legend.title=element_blank()) #turn off legen title gp <- gp + theme(legend.position="bottom") gp  • The results you show in your plots don't make sense to me, especially the lower two plots. If you calculate a$100(1-\alpha)\%$CI for normally distributed data using either a normal or t distributed sampling distribution, all distributions involved are symmetric and therefore the CI will be equally likely to overestimate the true mean as it is to underestimate it. Moreover, this property should be independent of the sample size for a fixed$\alpha$level. Sep 15, 2014 at 22:36 • Thanks for the comment. Yeah, that's exactly my question. Indeed I expected the "CI [to] be equally likely to overestimate the true mean as it is to underestimate it", but the results does not show that. Of course, my first thought was "I am doing something wrong" and then I checked my code literally a thousand times (well, not literally 1k, but surely many many times) and I could not find anything wrong. My last test was to use R built-in function to calculate the confidence interval and compare it with the CI calculated "by hand" and they are always exactly the same. Sep 16, 2014 at 5:00 • Your checking might not take into account methodological errors. When I performed a simulation of my own, using Mathematica, my empirical results demonstrated that, in fact, underestimation and overestimation occurred with roughly equal frequency for symmetric distributions, and was constant across sample size. I have not taken a close look at your R code. Sep 16, 2014 at 5:15 1 Answer First of all, I agree with the comments left by heropup. I'll add some details. The reason why your simulation breaks down may be a little subtle. At least I spend some time reading your code to find the source of the problem. Please notice, that you only simulate once for each of the cases. Then your CIs functions resample this initial data set. This clearly gives a lot of dependence between all of the samples. For instance, if you draw a sample of 1000 of the original data set, there is only one way to do this. If you draw a sample of 999 an overwhelming majority of the data set will still be the same between resamples. You'll need to do independent resampling. Otherwise, the 100 samples are essentially the same when you let$n$get large. Turning to your question, a confidence interval as the ones you do above are based on a distributional assumption, for instance your observations are normally distributed. If that's the case the confidence interval will be 'centered' in the sense you talk about when you construct the confidence interval symmetrically. This is evident from the symmetry of the distribution and of the procedure of constructing a confidence interval. In the above, you also calculate a confidence interval when the distributional assumption you make is not correct. Then an confidence interval need not be centered, even if the distribution assumed is symmetric. This can be seen simulated observations from a chi squared distribution and calculating confidence intervals based on a normal distribution. However, using a central limit theorem we can argue that the mean of the chi squared observations will be approximately normally distributed for large enough sample sizes. Finally, I just want to note that a confidence interval (or more generally, a confidence set) is basically and loosely speaking just some subset of a parameter set such that when you calculate this set a lot of times (hypothetically) it will contain the true parameter value for example 95% of the times. There's no claim of this set being 'centered' or symmetric around a parameter estimate. It can be chosen to have all sorts of strange forms. This is just not very intuitive and most of the time not very helpful. • Many thanks.I wanted to simulate what you say in your last paragraph in the sense that if you calculate the CI "a lot of times (hypothetically) it will contain the true parameter value for example 95% of the times". That is why the CI functions calculate the CI many times (num_CIs=100), each of them using a different sample of the same population. Since it is a new sample every time, I assumed them independent of each other; after your comment I see the dependency, but how can I do it differently? it seems to me that this dependency arises as a result of sampling the same population, isn't it? Sep 16, 2014 at 10:31 • Yes, if you sample the same finite population as you do, there'll be that kind of dependency. But if you just recall e.g.${\tt rchisq }\$ each time to make independent samples, there'll be no problem. You want to estimate the probability of getting a confidence interval that contains the true parameter. This should be done using independent samples as you are basically using the law of large numbers to estimate a probability (the mean of a indicator function).
– swmo
Sep 16, 2014 at 10:46
• As you say, I expected 95 out the 100 CIs to contain the true parameter. And they actually do roughly in those proportions (when the distributional assumption holds or when it does not hold but the central limit theorem kick-in). And this was actually the original motivation for my exercise, but then I found this discrepancy between the below and above mistakes, and I cannot seem to get a handle on it. Sep 16, 2014 at 10:46
• Thanks again, ..., one question. Recalling rchisq would be like generating a new universe (finite population). isn't it? I thought that to simulate the concept of a CI, the best way would be to take many -seemingly independent- samples from the same population. But yeah, maybe the way to think on this would be to take the whole distribution to be the population. I will try that. Sep 16, 2014 at 10:51
• Many many thanks, ..., I cannot vote because I am new here and my reputation is still not high enough, but I really appreciate your answer and comments. I'll work on it and I'll post my new results later. Sep 16, 2014 at 10:53