https://github.com/pymc-devs/pymc-examples/blob/main/examples/bart/bart_categorical_hawks.ipynb
(bart_categorical)=
:::{post} May, 2024 :tags: BART, regression :category: beginner, reference :author: Pablo Garay, Osvaldo Martin :::
In this example, we will model outcomes with more than two categories.
:::{include} ../extra_installs.md
:::
Here we will use a dataset that contains information about 3 species of hawks (CH=Cooper's, RT=Red-tailed, SS=Sharp-Shinned). This dataset has information for 908 individuals in total, each one containing 16 variables, in addition to the species. To simplify the example, we will use the following 5 covariables:
Also we are going to eliminate the NaNs in the dataset. With these we will predict the "Species" of hawks, in other words, these are our dependent variables, the classes we want to predict.
The following compares covariables to get a rapid data visualization for the 3 species.
It can be seen that the RT species have distributions more differentiated than the other two in almost all covariables, and the covariables wing, weight, and culmen present certain separations between species. Still, none of the variables have a marked separation among the species distributions such that they can cleanly separate them. It is possible to make a combination of covariables, probably wing, weight, and culmen, to achieve the classification. These are the principal reasons for realizing the regression.
First, we are going to prepare the data for the model, using "Species" as response and, "Wing", "Weight", "Culmen", "Hallux" and "Tail" as predictors. Using pd.Categorical(Hawks['Species']).codes we can codify the name species into integers between 0 and 2, being 0="CH", 1="RT" and 2="SS".
We only can have an instance of {class}~pymc_bart.BART() in each pymc model (for now), so to model 3 species we can use coordinate and dimension names to specify the shapes of variables, indicating that there are 891 rows of information for 3 species. This step facilite the later selection of groups from the InferenceData.
In this model we use the pm.math.softmax() function, for $\mu$ from pmb.BART(), because guarantees that the vector sums to 1 along the axis=0 in this case.
Now fit the model and get samples from the posterior.
It may be that some of the input variables are not informative for classifying by species, so in the interest of parsimony and in reducing the computational cost of model estimation, it is useful to quantify the importance of each variable in the dataset. PyMC-BART provides the function {func}~pymc_bart.plot_variable_importance(), which generates a plot that shows on his x-axis the number of covariables and on the y-axis the R$^2$ (the square of the Pearson correlation coefficient) between the predictions made for the full model (all variables included) and the restricted models, those with only a subset of the variables. The error bars represent the 94 % HDI from the posterior predictive distribution.
It can be observed that with the covariables Hallux, Culmen, and Wing we achieve the same $R^2$ value that we obtained with all the covariables, this is that the last two covariables contribute less than the other three to the classification. One thing we have to take into account in this is that the HDI is quite wide, which gives us less precision on the results; later we are going to see a way to reduce this.
We can also plot the scatter plot of the submodels' predictions to the full model's predictions to get an idea of how each new covariate improves the submodel's predictions.
Let's check the behavior of each covariable for each species with pmb.plot_pdp(), which shows the marginal effect a covariate has on the predicted variable, while we average over all the other covariates. Since our response variable is categorical, we'll pass softmax as the inverse link function to plot_pdp.
You can see we have to be careful with the softmax function, because it's not vectorized: it considers relationships between elements, so the specific axis along which we apply it matters. By default, scipy applies to all axes, but we want to apply it to the last axis, since that's where the categories are. To make sure of that, we use np.apply_along_axis and pass it in a lambda function.
The Partial Dependence Plot, together with the Variable Importance plot, confirms that Tail is the covariable with the smaller effect over the predicted variable: in the Variable Importance plot, Tail is the last covariate to be added and does not improve the result; in the PDP plot Tail has the flattest response.
For the rest of the covariate in this plot, it's hard to see which of them have more effect over the predicted variable, because they have great variability, showed in the HDI width.
Finally, some variability depends on the amount of data for each species, which we can see in the counts of each covariable for each species:
Now we are going to compare the predicted data with the observed data to evaluate the fit of the model, we do this with the Arviz function az.plot_ppc().
We can see a good agreement between the observed data (black line) and those predicted by the model (blue and orange lines). As we mentioned before, the difference in the values between species is influenced by the amount of data for each one. Here there is no observed dispersion in the predicted data as we saw in the previous two plots.
Below we see that the in-sample predictions provide very good agreement with the observations.
So far we have a very good result concerning the classification of the species based on the 5 covariables. However, if we want to select a subset of covariable to perform future classifications is not very clear which of them to select. Maybe something sure is that Tail could be eliminated. At the beginning when we plot the distribution of each covariable we said that the most important variables to make the classification could be Wing, Weight and, Culmen, nevertheless after running the model we saw that Hallux, Culmen and, Wing, proved to be the most important ones.
Unfortunately, the partial dependence plots show a very wide dispersion, making results look suspicious. One way to reduce this variability is adjusting independent trees, below we will see how to do this and get a more accurate result.
The option to fit independent trees with pymc-bart is set with the parameter pmb.BART(..., separate_trees=True, ...). As we will see, for this example, using this option doesn't give a big difference in the predictions, but helps us to reduce the variability in the ppc and get a small improvement in the in-sample comparison. In case this option is used with bigger datasets you have to take into account that the model fits more slowly, so you can obtain a better result at the expense of computational cost. The following code runs the same model and analysis as before, but fitting independent trees. Compare the time to run this model with the previous one.
Now we are going to reproduce the same analyses as before.
Comparing these two plots with the previous ones shows a marked reduction in the variance for each one. In the case of pmb.plot_variable_importance() there are smallers error bands with an $R^{2}$ value closer to 1. And for pmb.plot_pdp() we can see thinner HDI bands. This is a representation of the reduction of the uncertainty due to adjusting the trees separately. Another benefit of this is that the behavior of each covariable for each one of the species is more visible.
With all these together, we can select Hallux, Culmen, and, Wing as covariables to make the classification.
Concerning the comparison between observed and predicted data, we obtain the same good result with less uncertainty for the predicted values (blue lines). And the same counts for the in-sample comparison.
:::{bibliography} :filter: docname in docnames :::
:::{include} ../page_footer.md :::