Extracting the data

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

df = data.frame(
  year = 2009:2023,
  rate = c(1.36, 2.04, 3.44, 2.79, 3.04, 4.73, 2.56, 3.51, 
           2.99, 2.81, 2.82, 6.59, 4.57, 5.02, 3.90)
)

Plotting and linear regression

Now we can plot the points and also to show the linear regression line:

library(ggplot2)

df |> ggplot(aes(x=year, y=rate)) +
  geom_point() +
  geom_smooth(method=lm) +
  theme_classic()

## `geom_smooth()` using formula = 'y ~ x'

However, if we want to show the exact coefficients:

rate.lm = lm(rate ~ year, data = df)
summary(rate.lm)

## 
## Call:
## lm(formula = rate ~ year, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.2199 -0.7689  0.0320  0.2869  2.3629 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -374.09000  126.62713  -2.954   0.0112 *
## year           0.18729    0.06281   2.982   0.0106 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.051 on 13 degrees of freedom
## Multiple R-squared:  0.4061, Adjusted R-squared:  0.3605 
## F-statistic: 8.891 on 1 and 13 DF,  p-value: 0.0106

Confidence intervals (for the mean)

new.data <- data.frame(year = 2025)
predict(rate.lm, newdata = new.data, interval = 'confidence')

##        fit      lwr      upr
## 1 5.163571 3.808887 6.518256

Prediction intervals

predict(rate.lm, newdata = new.data, interval = 'prediction')

##        fit      lwr      upr
## 1 5.163571 2.519554 7.807589

Confidence interval for the coefficients

confint(rate.lm)

##                     2.5 %       97.5 %
## (Intercept) -647.65129323 -100.5287068
## year           0.05159094    0.3229805

Calculating confidence intervals

mean.pred.intervals <- function(x, y, pred.x) {
  n <- length(y) # Find sample size
  lm.model <- lm(y ~ x) # Fit linear model
  y.fitted <- lm.model$fitted.values # Extract the fitted values of y
  
  # Coefficients of the linear model, beta0 and beta1
  b0 <- lm.model$coefficients[1]
  b1 <- lm.model$coefficients[2]
  
  pred.y <- b1 * pred.x + b0 # Predict y at the given value of x (argument pred.x)
  
  # Find SSE and MSE
  sse <- sum((y - y.fitted)^2)
  mse <- sse / (n - 2)
  
  t.val <- qt(0.975, n - 2) # Critical value of t
  
  mean.se.fit <- (1 / n + (pred.x - mean(x))^2 / (sum((x - mean(x))^2))) # Standard error of the mean estimate
  pred.se.fit <- (1 + (1 / n) + (pred.x - mean(x))^2 / (sum((x - mean(x))^2))) # Standard error of the prediction
  
  # Mean Estimate Upper and Lower Confidence limits at 95% Confidence
  mean.conf.upper <- pred.y + t.val * sqrt(mse * mean.se.fit)
  mean.conf.lower <- pred.y - t.val * sqrt(mse * mean.se.fit)
  
  # Prediction Upper and Lower Confidence limits at 95% Confidence
  pred.conf.upper <- pred.y + t.val * sqrt(mse * pred.se.fit)
  pred.conf.lower <- pred.y - t.val * sqrt(mse * pred.se.fit)
  
  # Beta 1 Upper and Lower Confidence limits at 95% Confidence
  b1.conf.upper <- b1 + t.val * sqrt(mse) / sqrt(sum((x - mean(x))^2))
  b1.conf.lower <- b1 - t.val * sqrt(mse) / sqrt(sum((x - mean(x))^2))
  
  # Build data.frame of upper and lower limits calculated above, as well as the predicted y and beta 1 values
  upper <- data.frame(rbind(round(mean.conf.upper, 3), round(pred.conf.upper, 3), round(b1.conf.upper, 3)))
  lower <- data.frame(rbind(round(mean.conf.lower, 3), round(pred.conf.lower, 3), round(b1.conf.lower, 3)))
  fit <- data.frame(rbind(round(pred.y, 3), round(pred.y, 3), round(b1, 3)))
  
  # Collect all into data.frame and rename columns
  results <- data.frame(cbind(lower, upper, fit), row.names = c('Mean', 'Prediction', 'Coefficient'))
  colnames(results) <- c('Lower', 'Upper', 'Fit')
  
  return(results)
}

mean.pred.intervals(df$year, df$rate, new.data)

##             Lower Upper   Fit
## Mean        3.809 6.518 5.164
## Prediction  2.520 7.808 5.164
## Coefficient 0.052 0.323 0.187

R-squared

How do we find it from the previous output?