Capture percentage
Confidence interval and capture percentage
March 6, 2016
####################################################################################################The 95% CI of the sample mean is an interval for which the true mean will fall 95% of the time, on repeated measurements (see: http://neurostatscyrilpernet.blogspot.it/). However, this says nothing on what is the chance that a 95% CI of a single sample mean will capture future 95% CIs. According to Cumming & Maillardet (Psychological Methods, 2006; 11, 217-227), “On average, a 95% CI will include just 83.4% of future replication means”. Someone calls this the Capture Percentage (CP). See for example Daniel Lakens (Eindhoven University of Technology, The Netherlands): http://daniellakens.blogspot.it/2016/03/the-difference-between-confidence.html
I used Lakens’ simulation to see whether the Capture Percentage is influenced by sample size. I made minimal changes to Lakens’ codes. The final plot is mine (hence every error, if there is one).
What I’ve found? Essentially, the 95% CI, as calculated according the simulation, is insensitive to sample size: I mean, you will always found that in future measurement, 95% of CIs include the “true mean”. The Capture Percentage is lower, around 84%, and does not vary with sample size when the sample size is maintained constant in the simulation, albeit it does vary a bit more than the 95% CI.
library(ggplot2)
### nine repetated simulations with different sample size, from 10 to 810
par(mfcol = c (3,3))
### creation of an empty dataframe
data <- data.frame()
for (k in 1:9) {
N <- 100
for (j in 1:N) {
### original code by Daniel Lakens at https://gist.githubusercontent.com/Lakens/4ac8d88e5b7828a9f521/raw/49e79a98b1a30580068cdba910be38af02d12804/CI_vs_CP.R
## n=20 #set sample size
## nSims<-100000 #set number of simulations
### modified by me just for sample size and number of simulation
set.seed(j) ### set seed for reproducibility
n = 10*k*k #set sample size
nSims<-1000 #set number of simulations
x<-rnorm(n = n, mean = 100, sd = 15) #create sample from normal distribution
#95%CI
CIU<-mean(x)+qt(0.975, df = n-1)*sd(x)*sqrt(1/n)
CIL<-mean(x)-qt(0.975, df = n-1)*sd(x)*sqrt(1/n)
#plot data
#png(file="CI_mean.png",width=2000,height=2000, res = 300)
ggplot(as.data.frame(x), aes(x)) +
geom_rect(aes(xmin=CIL, xmax=CIU, ymin=0, ymax=Inf), fill="#E69F00") +
geom_histogram(colour="black", fill="grey", aes(y=..density..), binwidth=2) +
xlab("IQ") + ylab("number of people") + ggtitle("Data") + theme_bw(base_size=20) +
theme(panel.grid.major.x = element_blank(), axis.text.y = element_blank(), panel.grid.minor.x = element_blank()) +
geom_vline(xintercept=100, colour="black", linetype="dashed", size=1) +
coord_cartesian(xlim=c(50,150)) + scale_x_continuous(breaks=c(50,60,70,80,90,100,110,120,130,140,150)) +
annotate("text", x = mean(x), y = 0.02, label = paste("Mean = ",round(mean(x)),"\n","SD = ",round(sd(x)),sep=""), size=6.5)
#dev.off()
#Simulate Confidence Intervals
CIU_sim<-numeric(nSims)
CIL_sim<-numeric(nSims)
mean_sim<-numeric(nSims)
for(i in 1:nSims){ #for each simulated experiment
x<-rnorm(n = n, mean = 100, sd = 15) #create sample from normal distribution
CIU_sim[i]<-mean(x)+qt(0.975, df = n-1)*sd(x)*sqrt(1/n)
CIL_sim[i]<-mean(x)-qt(0.975, df = n-1)*sd(x)*sqrt(1/n)
mean_sim[i]<-mean(x) #store means of each sample
}
#Save only those simulations where the true value was inside the 95% CI
CIU_sim<-CIU_sim[CIU_sim<100]
CIL_sim<-CIL_sim[CIL_sim>100]
# cat((100*(1-(length(CIU_sim)/nSims+length(CIL_sim)/nSims))),"% of the 95% confidence intervals contained the true mean")
#Calculate how many times the observed mean fell within the 95% CI of the original study
mean_sim<-mean_sim[mean_sim>CIL&mean_sim<CIU]
# cat("The capture percentage for the plotted study, or the % of values within the observed confidence interval from",CIL,"to",CIU,"is:",100*length(mean_sim)/nSims,"%")
conf <- (100*(1-(length(CIU_sim)/nSims+length(CIL_sim)/nSims)))
capt <- 100*length(mean_sim)/nSims
### collect the data in a dataframe
data <- rbind(data, c(conf, capt))
names(data) <- c("95% CI", "Capture %")
}
### check the result
### head(data)
### Plot the results
set <- colMeans(data) ### calculate the mean on the long run
### the "\n" helps to print the title in two lines
plot(jitter(data[,1]),jitter(data[,2]), pch=19, col = "lightblue",
xlab = "95% CI", ylab = "Capture %",
main =(paste0("Include true mean: ", "95% CI = ", round(set[1],1), "%", "\n", "Capture percentage = ", round(set[2],1), "%")), cex.main= 1.25)
text(min(data[,1])+1, min(data[,2])+3, paste0("Sample size: n = ", n))
}
Let’ s see what is the Capture percentage of a 99.5% CI.
library(ggplot2)
### nine repetated simulations with different sample size, from 10 to 810
par(mfcol = c (3,3))
### creation of an empty dataframe
data <- data.frame()
for (k in 1:9) {
N <- 100
for (j in 1:N) {
### original code by Daniel Lakens at https://gist.githubusercontent.com/Lakens/4ac8d88e5b7828a9f521/raw/49e79a98b1a30580068cdba910be38af02d12804/CI_vs_CP.R
## n=20 #set sample size
## nSims<-100000 #set number of simulations
### modified by me just for sample size and number of simulation and the CI (99.5% rather than 95%)
set.seed(j) ### set seed for reproducibility
n = 10*k*k #set sample size
nSims<-1000 #set number of simulations
x<-rnorm(n = n, mean = 100, sd = 15) #create sample from normal distribution
#95%CI
CIU<-mean(x)+qt(0.9975, df = n-1)*sd(x)*sqrt(1/n)
CIL<-mean(x)-qt(0.9975, df = n-1)*sd(x)*sqrt(1/n)
#plot data
#png(file="CI_mean.png",width=2000,height=2000, res = 300)
ggplot(as.data.frame(x), aes(x)) +
geom_rect(aes(xmin=CIL, xmax=CIU, ymin=0, ymax=Inf), fill="#E69F00") +
geom_histogram(colour="black", fill="grey", aes(y=..density..), binwidth=2) +
xlab("IQ") + ylab("number of people") + ggtitle("Data") + theme_bw(base_size=20) +
theme(panel.grid.major.x = element_blank(), axis.text.y = element_blank(), panel.grid.minor.x = element_blank()) +
geom_vline(xintercept=100, colour="black", linetype="dashed", size=1) +
coord_cartesian(xlim=c(50,150)) + scale_x_continuous(breaks=c(50,60,70,80,90,100,110,120,130,140,150)) +
annotate("text", x = mean(x), y = 0.02, label = paste("Mean = ",round(mean(x)),"\n","SD = ",round(sd(x)),sep=""), size=6.5)
#dev.off()
#Simulate Confidence Intervals
CIU_sim<-numeric(nSims)
CIL_sim<-numeric(nSims)
mean_sim<-numeric(nSims)
for(i in 1:nSims){ #for each simulated experiment
x<-rnorm(n = n, mean = 100, sd = 15) #create sample from normal distribution
CIU_sim[i]<-mean(x)+qt(0.9975, df = n-1)*sd(x)*sqrt(1/n)
CIL_sim[i]<-mean(x)-qt(0.9975, df = n-1)*sd(x)*sqrt(1/n)
mean_sim[i]<-mean(x) #store means of each sample
}
#Save only those simulations where the true value was inside the 95% CI
CIU_sim<-CIU_sim[CIU_sim<100]
CIL_sim<-CIL_sim[CIL_sim>100]
# cat((100*(1-(length(CIU_sim)/nSims+length(CIL_sim)/nSims))),"% of the 99.5% confidence intervals contained the true mean")
#Calculate how many times the observed mean fell within the 99.5% CI of the original study
mean_sim<-mean_sim[mean_sim>CIL&mean_sim<CIU]
# cat("The capture percentage for the plotted study, or the % of values within the observed confidence interval from",CIL,"to",CIU,"is:",100*length(mean_sim)/nSims,"%")
conf <- (100*(1-(length(CIU_sim)/nSims+length(CIL_sim)/nSims)))
capt <- 100*length(mean_sim)/nSims
### collect the data in a dataframe
data <- rbind(data, c(conf, capt))
names(data) <- c("99.5% CI", "Capture %")
}
### check the result
### head(data)
### Plot the results
set <- colMeans(data) ### calculate the mean on the long run
### the "\n" helps to print the title in two lines
plot(jitter(data[,1]),jitter(data[,2]), pch=19, col = "pink",
xlab = "99.5% CI", ylab = "Capture %",
main =(paste0("Include true mean: ", "99.5% CI = ", round(set[1],1), "%", "\n", "Capture percentage = ", round(set[2],1), "%")), cex.main = 1.25)
text(min(data[,1])+0.5, min(data[,2])+3, paste0("Sample size: n = ", n))
}
If I have did not mess up with the calculations, the Capture percentage of the 99.5% CI of a sample mean is a pretty satisfying 95%, more or less, depending on the sample size.
So, if you use a stricter threshold (alpha at p < .005, as suggested by Johnson, PNAS, 2013; 110, 19313–19317) than the used one (alpha at p < .05), you can be “confident” than your results could be replicated (replicability) when you will reply the experiment in an indipendent sample from the same population, and as well be repeated (repeatability) when an independent researcher will reply the experiment in an independent sample from a similar population, and you can hope they could be generalized (generalizability) when the experiment will be repeated in an sample from a different population (sharing some point of continuity with the original population).
Of course, with a study based on a single sample you cannot know whether your findings are from the fortunate 95% (or 99.5%) side or from the unfortunate one.
Solutions? Trust n° 1: Trust no one, as Mulder would say.
from: http://it.fanpop.com/clubs/fox-mulder/images/25366898/title/fox-mulder-wallpaper
Never trust a single sample study….
References
Disclaimer: Not really a fan of Geoff Cumming.
Cumming G, Maillardet R. Confidence intervals and replication: where will the next mean fall? Psychol Methods. 2006 Sep;11(3):217-27.
Johnson VE. Revised standards for statistical evidence. Proc Natl Acad Sci USA. 2013 Nov 26;110(48):19313-7.
See also: Cumming G. Understanding replication: Confidence intervals, p values, and what’s likely to happen next time. In: Proceedings of the seventh internazional conference international conference on teaching statistics. International Association for Statistical Education. 2006.
http://iase-web.org/documents/papers/icots7/7D3_CUMM.pdf
I hope you’ve enjoyed this!
Have a nice day!
