{"id":3069,"date":"2023-05-15T05:49:08","date_gmt":"2023-05-15T05:49:08","guid":{"rendered":"http:\/\/optimumsportsperformance.com\/blog\/?p=3069"},"modified":"2023-05-15T12:02:51","modified_gmt":"2023-05-15T12:02:51","slug":"different-ways-of-calculating-intervals-of-uncertainty","status":"publish","type":"post","link":"https:\/\/optimumsportsperformance.com\/blog\/different-ways-of-calculating-intervals-of-uncertainty\/","title":{"rendered":"Different ways of calculating intervals of uncertainty"},"content":{"rendered":"

I’ve talked a lot in this blog about making predictions (see HERE<\/a><\/strong><\/span>, HERE<\/a><\/span><\/strong>, and HERE<\/a><\/span><\/strong>) as well as the difference between confidence intervals and prediction intervals and why you’d use one over the other (see HERE<\/a><\/span><\/strong>). Tonight I was having a discussion with a colleague about some models he was working on and he was building some confidence intervals around his predictions. That got me to thinking about the various ways we can code confidence intervals, quantile intervals, and prediction intervals in R. So, I decided to put together this quick tutorial to provide a few different ways of constructing these values (after all, unless we can calculate the uncertainty in our predictions, point estimate predictions are largely useless on their own).<\/span><\/p>\n

The full code is available on my GITHUB page<\/a><\/span><\/strong>.<\/p>\n

Load packages, get data, and fit regression model<\/strong><\/span><\/p>\n

The only package we will need is {tidyverse<\/strong>}, the data will be the mtcars<\/strong> dataset and the model will be a linear regression which attempts to predict mpg from wr and carb.<\/p>\n

\r\n## Load packages\r\nlibrary(tidyverse)\r\n\r\ntheme_set(theme_classic())\r\n\r\n## Get data\r\nd <- mtcars d %>%\r\n  head()\r\n\r\n## fit model\r\nfit_lm <- lm(mpg ~ wt + carb, data = d)\r\nsummary(fit_lm)\r\n<\/pre>\n
\"\"<\/a><\/p>\n
Get some data to make predictions on<\/strong><\/span><\/p>\n
We will just grab a random sample of 5 rows from the original data set and use that to make some predictions on.<\/p>\n
\r\n## Get a few rows to make predictions on\r\nset.seed(1234)\r\nd_sample <- d %>%\r\n  sample_n(size = 5) %>%\r\n  select(mpg, wt, carb)\r\n\r\nd_sample\r\n<\/pre>\n
\"\"<\/a><\/p>\n
Confidence Intervals with the predict() function<\/strong><\/span><\/p>\n
Using preidct() <\/strong>we calculate the predicted value with 95% Confidence Intervals.<\/p>\n
\r\n## 95% Confidence Intervals\r\nd_sample %>%\r\n  bind_cols(\r\n    predict(fit_lm, newdata = d_sample, interval = "confidence", level = 0.95)\r\n  )\r\n<\/pre>\n
\"\"<\/a>
\nCalculate confidence intervals by hand<\/strong><\/span><\/p>\n
Instead of using the R function, we can calculate the confidence intervals by hand (and obtain the same result).<\/p>\n
\r\n## Calculate the 95% confidence interval by hand\r\nlevel <- 0.95\r\nalpha <- 1 - (1 - level) \/ 2\r\nt_crit <- qt(p = alpha, df = fit_lm$df.residual) \r\n\r\nd_sample %>%\r\n  mutate(pred = predict(fit_lm, newdata = .),\r\n         se_pred = predict(fit_lm, newdata = ., se = TRUE)$se.fit,\r\n         cl95 = t_crit * se_pred,\r\n         lwr = pred - cl95,\r\n         upr = pred + cl95)\r\n<\/pre>\n
\"\"<\/a>
\nCalculate confidence intervals with the qnorm() function<\/strong><\/span><\/p>\n
Above, we calculated a 95% t-critical value for the degrees of freedom of our model. Alternatively, we could calculate 95% confidence intervals using the standard z-critical value for 95%, 1.96, which we obtain with the qnorm()<\/strong> function.<\/p>\n
\r\nd_sample %>%\r\n  mutate(pred = predict(fit_lm, newdata = .),\r\n         se_pred = predict(fit_lm, newdata = ., se = TRUE)$se.fit,\r\n         lwr = pred + qnorm(p = 0.025, mean = 0, sd = 1) * se_pred,\r\n         upr = pred + qnorm(p = 0.975, mean = 0, sd = 1) * se_pred)\r\n<\/pre>\n
\"\"<\/a><\/p>\n
Calculate quantile intervals via simulation<\/strong><\/span><\/p>\n
Finally, we can calculate quantile intervals by simulating predictions using the predicted value and standard error for each of the observations. We simulate 1000 times from a normal distribution and then use the quantile()<\/strong> function to get our quantile intervals.<\/p>\n
If all we care about is a predicted value and the lower and upper intervals, we can use the rowwise()<\/strong> function to indicate that we are going to do a simulation for each row and then store the end result (our lower and upper quantile intervals) in a new column.<\/p>\n
\r\n## 95% Quantile Intervals via Simulation\r\nd_sample %>%\r\n  mutate(pred = predict(fit_lm, newdata = .),\r\n         se_pred = predict(fit_lm, newdata = ., se = TRUE)$se.fit) %>%\r\n  rowwise() %>%\r\n  mutate(lwr = quantile(rnorm(n = 1000, mean = pred, sd = se_pred), probs = 0.025),\r\n         upr = quantile(rnorm(n = 1000, mean = pred, sd = se_pred), probs = 0.975))\r\n<\/pre>\n
\"\"<\/a><\/p>\n
While that is useful, there might be times where we want to extract the full simulated distribution. We can create a simulated distribution (1000 simulations) for each of the 5 observations using a for()<\/strong> loop.<\/p>\n
\r\n## 95% quantile intervals via Simulation with full distribution\r\nN <- 1000\r\npred_sim <- list()\r\n\r\nset.seed(8945)\r\nfor(i in 1:nrow(d_sample)){\r\n  \r\n  pred <- predict(fit_lm, newdata = d_sample[i, ])\r\n  se_pred <- predict(fit_lm, newdata = d_sample[i, ], se = TRUE)$se.fit\r\n  \r\n  pred_sim[[i]] <- rnorm(n = N, mean = pred, sd = se_pred)\r\n  \r\n}\r\n\r\nsim_df <- tibble( sample_row = rep(1:5, each = N), pred_sim = unlist(pred_sim) ) \r\n\r\nsim_df %>%\r\n  head()\r\n<\/pre>\n
\"\"<\/a>
\nNext we summarize the simulation for each observation.<\/p>\n
\r\n# get predictions and quantile intervals\r\nsim_df %>%\r\n  group_by(sample_row) %>%\r\n  summarize(pred = mean(pred_sim),\r\n         lwr = quantile(pred_sim, probs = 0.025),\r\n         upr = quantile(pred_sim, probs = 0.975)) %>%\r\n  mutate(sample_row = rownames(d_sample))\r\n<\/pre>\n
\"\"<\/a>
\nWe can then plot the entire posterior distribution for each observation.<\/p>\n
\r\n# plot the predicted distributions\r\nsim_df %>%\r\n  mutate(actual_value = rep(d_sample$mpg, each = N),\r\n         sample_row = case_when(sample_row == 1 ~ "Hornet 4 Drive",\r\n                                sample_row == 2 ~ "Toyota Corolla",\r\n                                sample_row == 3 ~ "Honda Civic",\r\n                                sample_row == 4 ~ "Ferrari Dino",\r\n                                sample_row == 5 ~ "Pontiac Firebird")) %>%\r\n  ggplot(aes(x = pred_sim)) +\r\n  geom_histogram(color = "white",\r\n                 fill = "light grey") +\r\n  geom_vline(aes(xintercept = actual_value),\r\n             color = "red",\r\n             size = 1.2,\r\n             linetype = "dashed") +\r\n  facet_wrap(~sample_row, scale = "free_x") +\r\n  labs(x = "Predicted Simulation",\r\n       y = "count",\r\n       title = "Predicted Simulation with actual observation (red line)",\r\n       subtitle = "Note that the x-axis are specific to that simulation and not the same")\r\n<\/pre>\n
\"\"<\/a><\/p>\n
Prediction Intervals with the predict() function<\/strong><\/span><\/p>\n
Next we turn attention to prediction intervals, which will be wider than the confidence intervals because they are incorporating additional uncertainty.<\/p>\n
The predict()<\/strong> function makes calculating prediction intervals very convenient.<\/p>\n
\r\n## 95% Prediction Intervals\r\nd_sample %>%\r\n  bind_cols(\r\n    predict(fit_lm, newdata = d_sample, interval = "predict", level = 0.95)\r\n  )\r\n<\/pre>\n
\"\"<\/a><\/p>\n
Prediction Intervals from a simulated distribution<\/strong><\/span><\/p>\n
Similar to how we simulated a distribution for calculating quantile intervals, above, we will perform the same procedure here. The difference is that we need to get the residual standard error (RSE) from our model as we need to add this additional piece of uncertainty (on top of the predicted standard error) to each of the simulated predictions.<\/p>\n
\r\n## 95% prediction intervals from a simulated distribution \r\n# store the model residual standard error\r\nsigma <- summary(fit_lm)$sigma\r\n\r\n# run simulation\r\nN <- 1000\r\npred_sim2 <- list()\r\n\r\nset.seed(85)\r\nfor(i in 1:nrow(d_sample)){\r\n  \r\n  pred <- predict(fit_lm, newdata = d_sample[i, ])\r\n  se_pred <- predict(fit_lm, newdata = d_sample[i, ], se = TRUE)$se.fit\r\n  \r\n  pred_sim2[[i]] <- rnorm(n = N, mean = pred, sd = se_pred) + rnorm(n = N, mean = 0, sd = sigma)\r\n  \r\n}\r\n\r\n# put results in a data frame\r\nsim_df2 <- tibble( sample_row = rep(1:5, each = N), pred_sim2 = unlist(pred_sim2) ) \r\n\r\nsim_df2 %>%\r\n  head()\r\n<\/pre>\n
\"\"<\/a><\/p>\n
We summarize our predictions and their intervals.<\/p>\n
\r\n# get predictions and intervals\r\nsim_df2 %>%\r\n  group_by(sample_row) %>%\r\n  summarize(pred = mean(pred_sim2),\r\n            lwr = quantile(pred_sim2, probs = 0.025),\r\n            upr = quantile(pred_sim2, probs = 0.975)) %>%\r\n  mutate(sample_row = rownames(d_sample))\r\n<\/pre>\n
\"\"<\/a><\/p>\n
Finally, we plot the simulated distributions for each of the observations.<\/p>\n
\"\"<\/a><\/p>\n
Wrapping Up<\/strong><\/span><\/p>\n
Uncertainty is important to be aware of and convey whenever you share your predictions. The point estimate prediction is one a single value of many plausible values given the data generating process. This article provided a few different approaches for calculating uncertainty intervals. The full code is available on my GITHUB page<\/a><\/span><\/strong>.<\/p>\n","protected":false},"excerpt":{"rendered":"
I’ve talked a lot in this blog about making predictions (see HERE, HERE, and HERE) as well as the difference between confidence intervals and prediction intervals and why you’d use one over the other (see HERE). Tonight I was having a discussion with a colleague about some models he was working on and he was […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[47,43],"tags":[],"class_list":["post-3069","post","type-post","status-publish","format-standard","hentry","category-model-building-in-r","category-sports-analytics"],"_links":{"self":[{"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/posts\/3069","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/comments?post=3069"}],"version-history":[{"count":1,"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/posts\/3069\/revisions"}],"predecessor-version":[{"id":3083,"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/posts\/3069\/revisions\/3083"}],"wp:attachment":[{"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/media?parent=3069"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/categories?post=3069"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/tags?post=3069"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}