Example code and comments for creating a confidence interval at a particular X0 for the mean of Y in a simple linear regression problem.

ConfidenceBandExample.R

# To illustrate the use of confidence bands, let's just make up 
# a data example. First, to ensure that everyone's simulated data
# looks the same in case you want to work together, let's tell R
# to use the same starting point for random number generation. 
# (don't worry about how this works)

set.seed(12345)

# Ok, now let's simulate some X values
# Don't worry about the function 'rpois' if you're not interested. 
# It just simulates random poisson variables (we aren't studying the 
# poisson distribution)
n = 200
X = rpois(n=n, 50)

# Now let's pick a true value for the slope and intercept

TrueBeta0 = 10
TrueBeta1 = 0.5

# Cool, now all that is left is to simulate Y with this TrueBeta and 
# some true value of sigma 

TrueSigma = 4
Y = rnorm(n=length(X), mean = TrueBeta0 + TrueBeta1*X, sd = TrueSigma)

# For good measure, let's plot X and Y along with their regression line

plot(X, Y, main = "Fancy Scatterplot")
myFancyRegression = lm(Y~X)
abline(reg = myFancyRegression)

# For convenience, let's grab the estimated Beta0 and Beta1, 
# along with the estimated sigma

EstBeta0 = myFancyRegression$coefficients[1]
EstBeta1 = myFancyRegression$coefficients[2]
sigma = summary(myFancyRegression)$sigma

# Alright, now let's pick an X value (called X0) at which to make a 
# confidence interval for the mean of Y|X

X0 = 55

# Let's draw a vertical line at X0 to remind ourselves of where we're looking:

abline(v=X0, lty = 2, col = "red")

# Now we can apply the formulas from the notes, but we need to pick 
# a true alpha 

alpha = 0.01
Yhat_at_X0 = EstBeta0 + EstBeta1*X0

# Try to come up with the following formula from the notes. You can do it!
SE_of_mean_at_X0 = (sigma* sqrt(1/n + ((X0-mean(X))^2)/((n-1)*var(X))))

# Finally, we can get the confidence interval 


T_quantiles = qt(c(alpha/2, 1-alpha/2), df = n-2)
CI_for_mean_at_X0 = Yhat_at_X0 + T_quantiles*SE_of_mean_at_X0

# Now let's proudly print our result:

NicelyFormattedOutput = paste("My CI for the mean of Y at X =",
                              X0,
                              "was: [",
                              round(CI_for_mean_at_X0[1],3),
                              ",",
                              round(CI_for_mean_at_X0[2],3),
                              "]")
print(NicelyFormattedOutput, quote = FALSE)



# Alternatively, you could just be boring and do 

print(CI_for_mean_at_X0)


# Last of all, it would be nice to see our CI on the graph. Let's plot It

abline(h=CI_for_mean_at_X0, lty = 3, col = "blue")

# Now update this to do the prediction band from the notes. It's easy, I promise!