{"id":2670,"date":"2022-10-17T01:43:44","date_gmt":"2022-10-17T01:43:44","guid":{"rendered":"http:\/\/optimumsportsperformance.com\/blog\/?p=2670"},"modified":"2022-11-08T03:32:58","modified_gmt":"2022-11-08T03:32:58","slug":"approximating-a-bayesian-posterior-prediction","status":"publish","type":"post","link":"https:\/\/optimumsportsperformance.com\/blog\/approximating-a-bayesian-posterior-prediction\/","title":{"rendered":"Approximating a Bayesian Posterior Prediction"},"content":{"rendered":"

This past week on the Wharton Moneyball Podcast<\/a><\/strong><\/span>, during Quarter 2 of the show the hosts got into a discussion about Bayesian methods, simulating uncertainty, and frequentist approaches to statistical analysis. The show hosts are all strong Bayesian advocates but at one point in the discussion, Eric Bradlow mentioned that “frequenstist can answer similar questions by building distributions from their model predictions.” <\/em>(paraphrasing)<\/p>\n

This comment reminded me of Chapter 7 in Gelman and Hill’s brilliant book, Data Analysis Using Regression and Multilevel\/Hierarchical Models<\/a><\/span><\/strong>. In this chapter, the authors do what they call an informal Bayesian approach<\/em> by simulating the predictions from a linear regression model. It’s an interesting (and easy) approach that can be helpful for getting a sense for how Bayesian posterior predictions work without building a full Bayesian model. I call it, “building a bridge to Bayes”.<\/em><\/p>\n

Since the entire book is coded in R, I decided to code an example in Python.<\/p>\n

The Jupyter notebook is accessible on my GITHUB page<\/a><\/span><\/strong>.<\/p>\n

(Special Material:<\/strong> In the Jupyter Notebook<\/a><\/span><\/strong> I also include additional material on how to calculate prediction intervals and confidence intervals by hand in python. I wont go over those entire sections here as I don’t want the post to get too long.<\/em>)<\/p>\n

Libraries & Data<\/strong><\/span><\/p>\n

We will start by loading up the libraries that we need and the data. I’ll use the iris data set, since it is conveniently available from the sklearn <\/strong>library. We will build the regression model using the statsmodels.api<\/strong> library.<\/p>\n

\r\n## import libraries\r\n\r\nfrom sklearn import datasets\r\nimport pandas as pd\r\nimport numpy as np\r\nimport statsmodels.api as smf\r\nfrom scipy import stats\r\nimport matplotlib.pyplot as plt\r\n\r\n## iris data set\r\niris = datasets.load_iris()\r\n\r\n## convert to pandas sata frame\r\ndata = iris['data']\r\ntarget = iris['target']\r\n\r\niris_df = pd.DataFrame(data)\r\niris_df.columns = ['sepal_length', 'sepal_width', 'petal_length', 'petal_width']\r\niris_df.head(3)\r\n<\/pre>\n
\"\"<\/a><\/p>\n
Build a linear regression model<\/strong><\/span><\/p>\n
Next, we build a simple ordinary least squares regression model to predict petal_width<\/strong> from petal_length<\/strong>.<\/p>\n
 <\/p>\n
\r\n## Set up our X and y variables\r\nX = iris_df['petal_length']\r\ny = iris_df['petal_width']\r\n\r\n# NOTE: statsmodels does not automatically add an intercept, so you need to do that manually\r\n\r\nX = smf.add_constant(X)\r\n\r\n# Build regression model\r\n# NOTE: the X and y variables are reversed in the function compared to sklearn\r\n\r\nfit_lm = smf.OLS(y, X).fit()\r\n\r\n# Get an R-like output of the model\r\n\r\nfit_lm.summary()\r\n<\/pre>\n
\"\"<\/a><\/p>\n
Simulate a distribution around the slope coefficient<\/strong><\/span><\/p>\n
The slope coefficient for the\u00a0petal_length<\/strong> variable, from the model output above, is 0.4158 with a standard error of 0.01. We will store these two values in their own variables and then use them to create a simulation of 10,000 samples and plot the distribution.<\/p>\n
\r\n## get summary stats\r\nmu_slope = 0.4158\r\nse_slope = 0.01\r\n\r\n## create simulation\r\nn_samples = 10000\r\ncoef_sim = np.random.normal(loc = mu_slope, scale = se_slope, size = n_samples)\r\n\r\n## plot simulated distribution\r\n\r\nplt.hist(coef_sim, bins = 60)\r\nplt.show()\r\n<\/pre>\n
\"\"<\/a><\/p>\n
We can also grab the summary statistics from the simulated distribution. We will snag the mean and the 90% quantile interval.<\/p>\n
\r\n## get summary stats from our simulation\r\nsummary_stats = {\r\n    'Mean': coef_sim.mean(),\r\n    'Low90': np.quantile(coef_sim, 0.05),\r\n    'High90': np.quantile(coef_sim, 0.95)\r\n}\r\n\r\nsummary_stats\r\n<\/pre>\n
\"\"<\/a><\/p>\n
Making a prediction on a new observation and building a posterior predictive distribution<\/strong><\/span><\/p>\n
Now that we’ve gotten a good sense for how to create a simulation in python, we can create a new observation of petal_length<\/strong> and make a prediction about what the petal_width<\/strong> would be based on our model. In addition, we will get the prediction intervals from the output and use them to calculate a standard error for the prediction, which we will use for the posterior predictive simulation.<\/p>\n
Technical Note:<\/strong> In the prediction output you will see that python returns a mean_ci_lower<\/strong><\/em> and mean_ci_upper<\/strong><\/em> and an obs_ci_lower <\/strong>and obs_ci_higher<\/strong>. The latter two are the prediction intervals\u00a0 while the former two are the confidence intervals. I previously discussed the difference HERE<\/a><\/span><\/strong><\/em> and this is confirmed in the Jupyter Notebook where I calculate these values by hand.<\/em><\/p>\n
\r\n# create a new petal_length observation\r\nnew_dat = np.array([[1, 1.7]])    # add a 1 upfront to indicate the model intercept\r\n\r\n# make prediction of petal_width using the model\r\nprediction = fit_lm.get_prediction(new_dat)\r\n\r\n# get summary of the prediction\r\nprediction.summary_frame()\r\n<\/pre>\n
\"\"<\/a><\/p>\n
Store the predicted value (0.343709) and then calculate the standard error from the lower and upper prediction intervals. Run a simulation and then plot the distribution of predictions.<\/p>\n
\r\nmu_pred = 0.343709\r\nse_pred = 0.754956 - 0.343709     # subtract the upper prediction interval from the mean to get the variability\r\nn_sims = 10000\r\n\r\npred_obs = np.random.normal(loc = mu_pred, scale = se_pred, size = n_sims)\r\n\r\nplt.hist(pred_obs, bins = 60)\r\nplt.show()\r\n<\/pre>\n
\"\"<\/a><\/p>\n
Just as we did for the simulation of the slope coefficient we can extract our summary statistics (mean and 90% quantile intervals).<\/p>\n
\r\n## get summary stats from our simulation\r\nsummary_stats = {\r\n    'Mean': pred_obs.mean(),\r\n    'Low90': np.quantile(pred_obs, 0.05),\r\n    'High90': np.quantile(pred_obs, 0.95)\r\n}\r\n\r\nsummary_stats\r\n<\/pre>\n
\"\"<\/a>Wrapping Up<\/strong><\/span><\/p>\n
That is a pretty easy way to get a sense for approximating a Bayesian posterior predictive distribution. Rather than simply reporting the predicted value and a confidence interval or prediction interval, it is sometimes nice to build an entire distribution. Aside from it being visually appealing, it allows us to answer other questions we might have, for example, what percentage of the data is greater or less than 0.5 (or any other threshold value you might be interested in)?<\/em><\/p>\n
As stated earlier, all of this code is accessible on my GITHUB page<\/a><\/span><\/strong> and the Jupyter notebook also has additional sections on how to calculate confidence intervals and prediction intervals by hand.<\/p>\n
If you notice any errors, let me know!<\/p>\n","protected":false},"excerpt":{"rendered":"
This past week on the Wharton Moneyball Podcast, during Quarter 2 of the show the hosts got into a discussion about Bayesian methods, simulating uncertainty, and frequentist approaches to statistical analysis. The show hosts are all strong Bayesian advocates but at one point in the discussion, Eric Bradlow mentioned that “frequenstist can answer similar questions […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[49,46,43,42],"tags":[],"class_list":["post-2670","post","type-post","status-publish","format-standard","hentry","category-bayesian-model-building","category-python-tips-tricks","category-sports-analytics","category-sports-science"],"_links":{"self":[{"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/posts\/2670","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=2670"}],"version-history":[{"count":3,"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/posts\/2670\/revisions"}],"predecessor-version":[{"id":2680,"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/posts\/2670\/revisions\/2680"}],"wp:attachment":[{"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/media?parent=2670"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/categories?post=2670"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/tags?post=2670"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}