Skip to main content
deleted 52 characters in body
Source Link
kjetil b halvorsen
  • 82.8k
  • 32
  • 201
  • 663
set.seed(1)
n_population <- 1e6
xx_population <- runif(n_population)
param <- 0.5
yy_population <- 2+param*xx_population+rnorm(n_population,0,0.5)

n_analyses <- 100
n_sample <- 30

CIs <- matrix(NA,nrow=n_analyses,ncol=3)
for ( ii in 1:n_analyses ) {
    index <- sample(1:n_population,n_sample)
    model <- lm(yy_population[index]~xx_population[index])
    CIs[ii,] <- c(confint(model)[2,1],coef(model)[2],confint(model)[2,2])
}

opar <- par(mai=c(.5,.1,.1,.1))
ii <- 1
plot(range(CIs),c(ii,ii),type="n",xlab="",ylab="",yaxt="n")
lines(CIs[ii,c(1,3)],rep(ii,2),col=2-(CIs[ii,1]<param&param<CIs[ii,3]))
points(CIs[ii,2],ii,pch=19,col=2-(CIs[ii,1]<0.5&0.5<CIs[ii,3]))
abline(v=param,lty=2,lwd=2)

plot(range(CIs),c(1,n_analyses),type="n",xlab="",ylab="",yaxt="n")
sapply(1:n_analyses,function(ii)lines(CIs[ii,c(1,3)],rep(ii,2),col=2-(CIs[ii,1]<param&param<CIs[ii,3])))
points(CIs[,2],1:n_analyses,pch=19,col=2-(CIs[,1]<0.5&0.5<CIs[,3]))
abline(v=param,lty=2,lwd=2)
set.seed(1)
n_population <- 1e6
xx_population <- runif(n_population)
param <- 0.5
yy_population <- 2+param*xx_population + rnorm(n_population, 0, 0.5)

n_analyses <- 100
n_sample <- 30

CIs <- matrix(NA, nrow=n_analyses, ncol=3)
for ( ii in 1:n_analyses ) {
    index <- sample(1:n_population,n_sample)
    model <- lm(yy_population[index]~xx_population[index])
    CIs[ii,] <- c(confint(model)[2,1],coef(model)[2],confint(model)[2,2])
}

opar <- par(mai=c(.5,.1,.1,.1))
ii <- 1
plot(range(CIs), c(ii,ii), type="n", xlab="", ylab="", yaxt="n")
lines(CIs[ii,c(1,3)], rep(ii,2), 
      col=2-(CIs[ii,1] < param&param < CIs[ii,3]))
points(CIs[ii,2], ii, pch=19, col=2-(CIs[ii,1]<0.5&0.5<CIs[ii,3]))
abline(v=param,lty=2,lwd=2)

plot(range(CIs),c(1,n_analyses),type="n",xlab="",ylab="",yaxt="n")
sapply(1:n_analyses,function(ii)lines(CIs[ii,c(1,3)],rep(ii,2),
        col=2-(CIs[ii,1]<param&param<CIs[ii,3])))
points(CIs[,2],1:n_analyses,pch=19,col=2-(CIs[,1]<0.5&0.5<CIs[,3]))
abline(v=param,lty=2,lwd=2)
set.seed(1)
n_population <- 1e6
xx_population <- runif(n_population)
param <- 0.5
yy_population <- 2+param*xx_population+rnorm(n_population,0,0.5)

n_analyses <- 100
n_sample <- 30

CIs <- matrix(NA,nrow=n_analyses,ncol=3)
for ( ii in 1:n_analyses ) {
    index <- sample(1:n_population,n_sample)
    model <- lm(yy_population[index]~xx_population[index])
    CIs[ii,] <- c(confint(model)[2,1],coef(model)[2],confint(model)[2,2])
}

opar <- par(mai=c(.5,.1,.1,.1))
ii <- 1
plot(range(CIs),c(ii,ii),type="n",xlab="",ylab="",yaxt="n")
lines(CIs[ii,c(1,3)],rep(ii,2),col=2-(CIs[ii,1]<param&param<CIs[ii,3]))
points(CIs[ii,2],ii,pch=19,col=2-(CIs[ii,1]<0.5&0.5<CIs[ii,3]))
abline(v=param,lty=2,lwd=2)

plot(range(CIs),c(1,n_analyses),type="n",xlab="",ylab="",yaxt="n")
sapply(1:n_analyses,function(ii)lines(CIs[ii,c(1,3)],rep(ii,2),col=2-(CIs[ii,1]<param&param<CIs[ii,3])))
points(CIs[,2],1:n_analyses,pch=19,col=2-(CIs[,1]<0.5&0.5<CIs[,3]))
abline(v=param,lty=2,lwd=2)
set.seed(1)
n_population <- 1e6
xx_population <- runif(n_population)
param <- 0.5
yy_population <- 2+param*xx_population + rnorm(n_population, 0, 0.5)

n_analyses <- 100
n_sample <- 30

CIs <- matrix(NA, nrow=n_analyses, ncol=3)
for ( ii in 1:n_analyses ) {
    index <- sample(1:n_population,n_sample)
    model <- lm(yy_population[index]~xx_population[index])
    CIs[ii,] <- c(confint(model)[2,1],coef(model)[2],confint(model)[2,2])
}

opar <- par(mai=c(.5,.1,.1,.1))
ii <- 1
plot(range(CIs), c(ii,ii), type="n", xlab="", ylab="", yaxt="n")
lines(CIs[ii,c(1,3)], rep(ii,2), 
      col=2-(CIs[ii,1] < param&param < CIs[ii,3]))
points(CIs[ii,2], ii, pch=19, col=2-(CIs[ii,1]<0.5&0.5<CIs[ii,3]))
abline(v=param,lty=2,lwd=2)

plot(range(CIs),c(1,n_analyses),type="n",xlab="",ylab="",yaxt="n")
sapply(1:n_analyses,function(ii)lines(CIs[ii,c(1,3)],rep(ii,2),
        col=2-(CIs[ii,1]<param&param<CIs[ii,3])))
points(CIs[,2],1:n_analyses,pch=19,col=2-(CIs[,1]<0.5&0.5<CIs[,3]))
abline(v=param,lty=2,lwd=2)
Source Link
Stephan Kolassa
  • 130.7k
  • 22
  • 264
  • 497

I'll be honest there: I don't think the actual distinction is all that important. Yes, saying that "the probability of the estimated parameter being included in the confidence interval is 95%" is incorrect, for the precise reason you give. However, I do not think it is a major problem. (I would be interested in any other point of view. Has this incorrect manner of writing ever led to "real" problems?)

If you run a single experiment and get a single CI, then yes, it either contains or does not contain the true value of the parameter:

single CI

As you write, there is no probability involved any more. The correct interpretation of a CI only comes in if we (explicitly or implicitly) run precisely the same experiment many times and collect all CIs:

100 CIs

And here, we see that (approximately) 95% of the CIs do contain the correct parameter. (The CI from the single experiment pictured above is the one at the bottom in this second plot.)

Yes, it would be better if everyone used the correct nomenclature, or at least had the correct interpretation involving many re-runs of the experiment in the back of their head while they were writing sloppily. But people don't.

And I honestly don't think this is a truly big deal.

R code:

set.seed(1)
n_population <- 1e6
xx_population <- runif(n_population)
param <- 0.5
yy_population <- 2+param*xx_population+rnorm(n_population,0,0.5)

n_analyses <- 100
n_sample <- 30

CIs <- matrix(NA,nrow=n_analyses,ncol=3)
for ( ii in 1:n_analyses ) {
    index <- sample(1:n_population,n_sample)
    model <- lm(yy_population[index]~xx_population[index])
    CIs[ii,] <- c(confint(model)[2,1],coef(model)[2],confint(model)[2,2])
}

opar <- par(mai=c(.5,.1,.1,.1))
ii <- 1
plot(range(CIs),c(ii,ii),type="n",xlab="",ylab="",yaxt="n")
lines(CIs[ii,c(1,3)],rep(ii,2),col=2-(CIs[ii,1]<param&param<CIs[ii,3]))
points(CIs[ii,2],ii,pch=19,col=2-(CIs[ii,1]<0.5&0.5<CIs[ii,3]))
abline(v=param,lty=2,lwd=2)

plot(range(CIs),c(1,n_analyses),type="n",xlab="",ylab="",yaxt="n")
sapply(1:n_analyses,function(ii)lines(CIs[ii,c(1,3)],rep(ii,2),col=2-(CIs[ii,1]<param&param<CIs[ii,3])))
points(CIs[,2],1:n_analyses,pch=19,col=2-(CIs[,1]<0.5&0.5<CIs[,3]))
abline(v=param,lty=2,lwd=2)