In the previous segment, we saw two outliers in the model relating the log of infant mortality to the log of income. Here, we're going to discuss the options available to you if you decide that these outliers really do belong in the dataset. You can't just delete them. The first approach is to look for additional covariants or explanatory variables that maybe able to explain the outliers. For example, there could be a number of variables that provide information about infant mortality above and beyond what income provides as an explanation. Looking back at our data, there are two variables we haven't used yet. The first one is region and the second is oil. The oil variable indicates oil exporting countries, both Saudi Arabia and Libya are oil exporting countries, so perhaps this might explain part of the anomaly. We've already coded up and run model 1 where log_income was the only explanatory variable. Here now is model 2 where we include the indicator variable is oil, meaning is this an oil exporting country along with a new coefficient for that term. We've added that coefficient, as well, to our prior and we've kept the same prior on our observation variance. To run model 2, we need to modify the data that we send to the model to include the is_oil indicator. To create that indicator, we need to transform the original data, which is in terms of yes and no. We want ones to represent yes and zeroes to represent no. So this says as.numeric turn the true or false statement into a number for whether dat$oil is a yes. If we run that, we get a set of 1's and 0's. Here is the batch of oil exporting countries right here. Everything else in the model setup is the same, except now we're going to draw three initial values for our beta vector, and we'll run mod2 with the new data and the new initial values. We'll run 1,000 iterations of burn-in, And run our samples to be saved. Let's go ahead and run all of this code. First, we'll run the rjags library, And then run model 2. Every time you ran a new model, you should always check the convergence diagnostics for the chain you just ran. We're going to skip that in the interest of time here but I will report that the chains look okay here. Let's take a look at the posterior summary from our chains. It looks like not much has changed since model 1. The intercept is a similar number, the effect related to income has not changed much, and we see a positive relationship between oil production and the log of infant mortality. Because this data are merely observational, we can't say that oil production causes an increase in infant mortality, indeed, that most certainly is not true. We can say that they are positively correlated. Now let's check the residuals. We'll follow the same steps we did for model 1 where we create the design matrix. We'll call it X2 in this case where we've now added the is_oil indicator. So we'll run that line and take a look at the top of X2. So now we have a third column for the oil indicator. We'll collect our posterior means for the parameters, as well as the predicted values based on those posterior means. And finally, we'll calculate the residuals as the difference between the value and the prediction. Now let's look at the residuals. We also want to compare these residuals to the residuals from model 1. So let's put both of these plots on the same screen with par(mfrow. We want two rows of plots and one column of plots. Start an empty plot and we'll fill it first with the residuals from the new model and then the residuals from the old model. It looks like we do have an improvement. The residuals in the new model are closer to the bulk of the data than the residuals were in the old model, where they are much further away. The standard deviation of the residuals in the second model here is 0.64. The strongest outlier up here is more than 2 away from the mean, which means it's more then 3 standard deviations, 3 times 0.64 away from the mean. So three standard deviations is still quite an extreme outlier. So although we've seen improvement, these outliers are pretty far away still. We might consider adding the other covariant region, but instead, let's look at another option when we are faced with strong outliers. That option is to change the distribution we use for the likelihood. The normal distribution or the normal likelihood has thin tails. Almost all of the probability in a normal distribution is concentrated within the first few standard deviations of the mean. This does not accommodate outliers well. Consequently, models with a normal likelihood might be overly influenced by outliers. Remember the t distribution is similar to the normal distribution. Here's a plot with the pdf of the normal in black and the pdf of a t distribution with one degree of freedom in red. The t distribution is similar to the normal distribution but it has thicker tails, which are better at accommodating outliers. The t model might look something like this. We're not going to run this model but we're going to step through it and talk about it a little bit. Notice that the t distribution, which we're using as our likelihood here, has three parameters. A location, in JAGS, tau is an inverse scale parameter and a degrees of freedom parameter. The smaller the degrees of freedom, the heavier the tails of the t distribution. We might fix the degrees of freedom to some number like we did in this plot or we can assign it a prior distribution, that's what we've done here. The degrees of freedom has to be a positive number, so we'll give it an exponential prior. The inverse scale, or the tau parameter, behaves similarly to the precision parameter in the normal model. It is not exactly the same as the precision parameter but it's close, so we'll give it the same dgamma prior that we used before. The inverse scale is related to the standard deviation of the errors by this equation here. Finally, we're going to relate the location parameter with our linear model form right here. As a side note, we should aware that the t distribution does not have a mean and a variance if the degrees of freedom is less than two. That can happen under this model. If we wanted to force the model, the likelihood, to have a mean and a variance, we could increase the degrees of freedom by two by adding a new variable. Let's call it nu and then the degrees of freedom would simply be nu + 2.0. This would guarantee that the degrees of freedom is greater than two. Again, we're not going to fit this model, we'll leave it up to you.
