{"id":2872,"date":"2023-01-25T14:23:29","date_gmt":"2023-01-25T14:23:29","guid":{"rendered":"http:\/\/optimumsportsperformance.com\/blog\/?p=2872"},"modified":"2023-01-25T14:23:29","modified_gmt":"2023-01-25T14:23:29","slug":"bayesian-linear-regression-getting-started-with-pymc3","status":"publish","type":"post","link":"https:\/\/optimumsportsperformance.com\/blog\/bayesian-linear-regression-getting-started-with-pymc3\/","title":{"rendered":"Bayesian Linear Regression: Getting started with PyMC3"},"content":{"rendered":"

Previously I’ve used {rstanarm<\/strong>}, {brms<\/strong>}, and Stan<\/strong> for fitting Bayesian models. However, as I continue to work on improving my Python skills, I figured I’d try and delve into the PyMC3<\/strong> framework for fitting such models. This article will go through the following steps:<\/p>\n

    \n
  1. Fitting the model<\/li>\n
  2. Making a point estimate prediction<\/li>\n
  3. Making a point estimate prediction with uncertainty<\/li>\n
  4. Calculating a posterior predictive distribution<\/li>\n<\/ol>\n

    I’ve covered the last three steps in a prior blog on making predictions with a Bayesian model<\/a><\/strong><\/span>. I know there are probably functions available in PyMC3<\/strong> that can do these things automatically (just as there are in {rstanarm<\/strong>}) but instead of falling back on those, I create the posterior distributions here using numpy<\/strong> and build them myself.<\/p>\n

    The entire code and data are available on my GITHUB page<\/a><\/span><\/strong>, where I also have the model coded in {rstanarm<\/strong>}, for anyone interested in seeing the steps in a different code language.<\/p>\n

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

    The data I’ll be using is the {mtcars<\/strong>} data set, which is available in R. I’ve saved a copy in .csv format so that I can load it into my Jupyter notebook.<\/p>\n

    \r\nimport pandas as pd\r\nimport numpy as np\r\nimport matplotlib.pyplot as plt \r\nimport seaborn as sns\r\nimport pymc3 as pm\r\nimport arviz as az \r\nimport os\r\n\r\n# load mtcars\r\nd = pd.read_csv('mtcars.csv')\r\nd.head()\r\n<\/pre>\n
    \"\"<\/a><\/p>\n
    Exploratory Data Analysis<\/strong><\/span><\/p>\n
    The model will regress mpg on engine weight (wt). Let’s plot and describe those two variables so that we have a sense for what we might be working with.<\/p>\n
    \"\"<\/a> \"\"<\/a><\/p>\n
    Linear Regression<\/strong><\/span><\/p>\n
    Before fitting the Bayesian model, I want to fit a simple regression model to see what the coefficients look like.<\/p>\n
    \r\nimport statsmodels.api as sm\r\n\r\nx = d['wt']\r\ny = d['mpg']\r\n\r\nx = sm.add_constant(x)\r\n\r\nfit = sm.OLS(y, x).fit()\r\n\r\nfit.summary()\r\n<\/pre>\n
    \"\"<\/a><\/p>\n
    We can see that for every one unit increase in engine weight the miles per gallon decrease, on average, by about 5.3.<\/p>\n
    Bayesian Regression (PyMC3)<\/strong><\/span><\/p>\n
    Fitting a Bayesian regression model in PyMC3 requires us to specify some priors. For this model I’ll use a prior intercept of 40 \u00b1 10 and a prior beta for the wt<\/strong> variable of 0 \u00b1 10. The thing to note here is that the priors I’m specifying priors were created by first looking at the data that has been collected (which is technically cheating). Normally we would have priors BEFORE<\/strong><\/em> collecting our data (using prior published research, data from a pilot study, prior intuition, etc) and then combine the prior with the observations to obtain a posterior distribution. However, the aim here is to understand how to code the model, so I’ll use these priors. I’ll write a longer blog on priors and what they do to a model in the coming weeks.<\/p>\n
    Some notes on fitting the model in PyMC3<\/strong>:<\/p>\n