4 Bayesian Regressions with Linear, Binomial Logistic, Poisson, and Exponential

 To set up a model with pymc, all you need to do is wrap the variables in a model object. Then the first thing you can do is create Data objects for each x and y. This isn't a nessessary step but it helps with predictions later on. You should specify that the x's are mutable and the y is not. I've seen many people use the x1, x2, x3,... styling but I prefer descriptive variable names as they'll be handy later when we make inference. 

 Next you can set variables for each beta coefficient you'll need. Here for brevity I typically shorten the variable names, but continue to give the descriptive names as a parameter. These will show up in the summary inference table later on, and it will be easier to have the descriptive names handy at that point. These are uninformative priors, normal distributions with very wide sigmas, with the mu centered at 0. The idea here is that since we will be running many iterations, and generating many samples, the uninformative priors will eventually update themselves to be near the true beta values. This is the flexibility in Bayesian modeling. You'll also need ot include a beta coefficent for the intercept of course. 

 Now all we have to do is put these varaibles together with the x covariate variables in a linear fashion, similar to what we would do in any linear regression. Notice how creating interaction terms is exactly how you would expect it to be. We can simply add, multiply, square, etc... the variables together. The likelihood, represents our y variable, in our case spend, which we can model was a flexible Normal distribution. For this variable we will use the linear combination we created as it's mean or mu. For the sigma of the likelihood, I'm actually choosing to model it with an uninformative Gamma hyper prior. Sometimes pymc can have divergence issues with really really small values. So using a Gamma hyper prior will help avoid these issues. We also pass the observed spend (y) values into the likelihood, into the observed parameter. 

 Finally we create the trace varaible which will hold all of our generated samples from the Monte Carlo process. By setting the tune parameter we are essentially throwing out the first N smaples. The purpose of these burn samples is to allow the uninformative priors time to learn, thereby creating more accurate estimates.  Typically I set it to about 1% of however many samples I'm generating. So here for 10,000 samples, I'll burn 1,000 samples.

 After the model has finsished generating it's samples, you can call a summary table through the arviz package in order to view the coefficients and their respective 95% Bayesian credible intervals, the Bayesian equivalent of a confidence interval. HDIs or Highest Density Intervals  are sometimes also refered to as HPD or Highest Posterior Distributions. Here I'm using the var_names and filer_vars parameters to get only the coefficients that start with "b".  

Looking into the variables whose 95% intervals do not contain 0, we find that:

• b1 - for every $1 increase in average_spend we can expect the spend to increase by $0.89 on average, keeping all other variables constant, with a 95% probability of falling into the interval of [$0.27, $1.53].  

• b3 - for every increase (of one time) in times_visited_this_week we can expect the spend to increase by $35.49 on average, keeping all other variables constant, with a 95% probability of falling into the interval of [$32.21, $38.79].  

• b4 - for every increase of one minute in view_time we can expect the spend to increase by $15.93 on average, keeping all other variables constant, with a 95% probability of falling into the interval of [$10.02, $22.03].  

It's also really simple to add in interaction terms. Just create new uninformative priors for these terms, and add them into the linear model.

Although adding in the interaction terms can confuscate the interpretation of coefficients, we can still  gain some valuable insights. 

Interestingly, you can see that view_time * times_visited_this_week resulted in a negative coefficient, while it's parts both resulted in positive coefficients. The interaction being negative here alludes to the fact that perhaps users who visit often and view too long are getting cold feet. This could be something to look further into, in order to try to retain these sales. Maybe these users would be a good candidate for our coupons. 

 Another inntersting thing we can do with our model, is that we can make predictions. For instance let's say that we want to know the average spend for a particular type of user. Let's say that this user typically spends about $30, it's not their first time on the site, and they've already visted the site 2 times this week, and looked at a product for 2 minutes.  Maybe we want to targetthis user with a coupon. We can then simply update our model with the set_data function, and regenerate the summary table. 

Now we can say that we expect this user to spend about $190.19 with a 95% proability of spending in the interval of [75.43, 304.87]. Isn't that a much nicer way to make inference, rather than say with a p-value? With the Bayesian credible intervals you can communicate your finding on the parameter in terms of direct probability. 

 We can plot our posterior to vizualy see the distrubtuions of each beta coefficient, along with their 95% credible intervals. 

 The very normal shape of each of those beta coefficients is a good check that the model converged well. To evaluate the model you can generate samples posterior predictive distribution, then compare them to the actual known spend distribution. An familiar metric here is R2, which you can call from the arviz package. 

 Another check for you model would be to plot the trace. You basically want to see lot's of fuzzy caterpillars, which look exactly like what we've got here.

 And here's a DAG of our simple Linear regression model. 

 The Binomial distribution models the number of sucesses in a certain number of trials and takes two paramters n (the number of trials) and p (the probability of success). We are also given the length of the melanoma (measured in mm) in our data which we will use as our x explanatory variable.  Here we set up the pymc model the same way as before. Even though we are only using a singular x, we will also create Data objects for s, the number of survivals, since we can pass this into the Binnomial distribution as an observed value. Same thing with the total number of experiments, n. Here, as before, I'm using uninformative priors for all beta coefficients. 

 Now the likelihood here will be the Binomial distribution, and we can pass the number of experiemnts to n and the number of observed survivals to the observed parameter. It also requires p, the probability of success (survival). For this we can  use the logistic function, aka  the sigmoid function! And within the logistic function we can pass the beta coefficients and expalanatory variable, length. These will update and generate the approprite probabilities of survival. Wow! 

With pymc we also have the option of modeling the invlogit() function if we wish to track this varaible, or we can use the logit_p parameter in the Binomial function and just simply pass in the linear equation. So flexible! By using the invlogit() fucntion we will be modeling the direct probability, not the log odds, so when we make a predicition there wll be no need to apply exp().

 Calling the summary function we can see that the length of the melanoma negatively effects the patients chances for survival. We can see that since the HDI doesn't contain 0 this is statistically significant. So for every unit increase in length (1 mm), the posterior probability of survival decreases by 46.8%.  

 And we can vizualize that interval, with the plot posterior function. 

  Now what if we wanted to make a prediction based off of our model? Let's say that a patient walks into our office and has a possible melanoma growth of 2.75mm. We can use the coefficients of our model to make that prediction. But to do this we'll need to wrap our linear function in the sigmoid function. We can use scipy.special.expit() for our logistic sigmoid function. And that patient will have a 91.8% probability of survival. 

This time I can squeeze the entire model into one frame. Similar to the Logistic regression from before we can use the linear function wrapped into the appropriate link function. For the Poisson distribution, this is the exponential function. We can call it from the pymc.math library. We're also using uninformative priors here, as they will update pretty quickly. 

 Calling the summary function with a 95% HDI, we can see that both chlorine and radium are significant predictors of cancer here, but in opposite dierctions. It seems that for every unit increase in cancer (1 person) there is an increase of 0.024 ppb int he water in that city. 

 And let's say that we'd like to test out our own cities water and we find that there is 3.24 ppb chlorine and 8.16 ppb radium. Well as it turn out we could expecet 26 people to develop this particular cancer this year. And we can say that there is a 95% proability that the number of people who develop cancer will be in the interal of [16.068,  35.947].

 Now I'll set some variables to help the model later on. Min, max, mean, and standard deviation. Then I'm going to create numpy masks to help the model understand that there is missing data. 

The model follows a similar structure as the previous models. For the variable that is not the likelihood (days), I will model it with a Normal distribution, specifying the sample mean and standard deviation to set it to about the right place. The pymc model will impute the missing data automatically. Yes, that's right - automatically. I can do this easily because we've alreay specified the distribution of the parameter, and Bayesian modeling relies heavily on sampling. All of the other coefficients are set to be uninformative Normal priors. The mu variable is set to an exponential equation, rather than linear, to model population decline. And the likelihood (number of salmon) are modeled with a Truncated Normal distribution. For this distributions you can specify cutoffs for the distribution. And in our case it deosn't make contextual sense to have negative number of salmon in the lakes. For the maximum, I couldv'e chosen to use the None argument in the paramter, but instead I've set it to the maximum number of salmon + 10% for wiggle room. 

 Here's the summary. And rather than focus on the coefficients, let's focus on these imputed values. Each missing value is returned for each varaible as a mean with a 95% HDI. Amazing. 

 Now I'm creating an new column to flag the vales as imputed, and also using some simple for loops to impute the missing values that were retunred in the summary. 

 And finally we can vizuallize our imputed data points along side our actuals. The imputed values seem to be right in line with what I would've expected their actual values to be. Notably, you can see that they don't fall on an unrealistic straight line, but rather are scattered according to the distribution they were sampled from. It's much more natural. You can even see that the imputation method can not only handle the interpolations, but also exterpolations as seen in the point to the far right.  That point was 50 days away from any other actual data point and the model imputed it perfectly.