---
title: "Introduction to bbreg"
author: "Wagner Barreto-Souza, Vinicius D. Mayrink, Alexandre B. Simas"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to bbreg}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(bbreg)
```

# Brief introduction

The *bbreg* package deals with regression models with response variables being continuous
and bounded. It currently provides implementation of two regression models: beta
regression <https://www.tandfonline.com/doi/abs/10.1080/0266476042000214501> and bessel regression <https://doi.org/10.1111/anzs.12354> and <https://arxiv.org/abs/2003.05157>. For both of these models, the estimation is 
done with the EM-algorithm. The EM-algorithm approach for beta regression
was developed in <https://www.tandfonline.com/doi/abs/10.1080/00949655.2017.1350679>.

# Choosing between bessel regression and beta regression

There are several manners to choose between bessel and beta regression, most of them using diagnostic analysis. 
One approach we would like to call attention is the **discrimination test** developed in the bessel regression paper <https://doi.org/10.1111/anzs.12354> and <https://arxiv.org/abs/2003.05157>. This discrimination test has a very nice performance on selecting
the "correct" distribution (in the sense that it performs well with synthetic data)
and is implemented in the *bbreg* package. 

Another approach to select between bessel and beta regression models is to use quantile-quantile plots
with simulated envelopes. The *bbreg* package also has a function to build such plots. 

We will show below how to use the discrimination test and how to construct the quantile-quantile plot
with simulated envelopes.

# Using the bbreg package

From version 2.0.0 onward, the usage of the *bbreg* package is standard in the
sense that it uses *formula* to handle data frames and that it has several
standard S3 methods implemented, 
such as *summary*, *plot*, *fitted*, *predict*, etc.

If the regression model is not specified (via the *model* parameter), the discrimination test
is used to select the regression model to be used.

## Fitting *bbreg* (bessel or beta regression) models

Let us fit a model:

```{r}
fit <- bbreg(agreement ~ priming + eliciting, data = WT)
```

Let us now see its output (notice that it provides a nice print when the fitted object is called):

```{r}
fit
```

It is also noteworthy that the bessel regression model was chosen based on the discrimination criteria. 
If we want to fit a beta regression model (with estimates obtained from the EM algorithm, which are
more robust than the usual maximum likelihood estimates due to the flatness of the likelihood function
with respect to the precision parameter), we pass the argument *model = "beta"*:
```{r}
fit_beta <- bbreg(agreement ~ priming + eliciting, data = WT, model = "beta")
fit_beta
```

Observe that the *bbreg* package makes it explicit that the discrimination test was ignored and the 
beta regression was perfomed.

Let us now add a precision covariate. Let us add the covariate *priming* as precision covariate.
To this end, we add "*| priming*" to the end of the formula:
```{r}
fit_priming <- bbreg(agreement ~ priming + eliciting | priming, data = WT)
fit_priming
```

To add both *priming* and *eliciting* as precision covariates, we add "*|priming + eliciting*"
to the end of the formula:
```{r}
fit_priming_eliciting <- bbreg(agreement ~ priming + eliciting | priming + eliciting, data = WT)
fit_priming_eliciting
```

## Viewing summaries

To view a summary of the fitted model, we simply call the method *summary*
to the fitted *bbreg* object:

```{r}
summary(fit)
```
For the beta regression fit:
```{r}
summary(fit_beta)
```
For the fit with priming as covariates:
```{r}
summary(fit_priming)
```

## Changing link functions

The *bbreg* package allows for several link functions for the mean and precision
parameters. If the link function for the mean is not specified, the *bbreg* package will
consider a *logit* link function for the mean parameter. If the link function
for the precision parameter is not specified and there are no covariates for the
precision parameter, then the *identity* link function will be used. If the link
function for the precision parameter is not specified and there are covariates
for the precision parameter, then the *log* link function will be used.

To change the link function for the mean parameter, just change the *link.mean* argument.
The currently implemented link functions for the mean parameter are
*logit*,*probit*, *cauchit*, *cloglog*.

To change the link function for the precision parameter, just change the *link.precision* argument.
The currently implemented link functions for the precision parameter are
*logit*,*probit*, *cauchit*, *cloglog*.

In the next example we will consider the regression model in the **fit** example
above, but choosing the *cloglog* as link function for the mean:

```{r}
fit_cloglog <- bbreg(agreement ~ priming + eliciting, data = WT, link.mean = "cloglog")
fit_cloglog
```

Now, we will also set the link function for the precision parameter to be *sqrt*:
```{r}
fit_cloglog_sqrt <- bbreg(agreement ~ priming + eliciting, data = WT, 
                          link.mean = "cloglog", link.precision = "sqrt")
fit_cloglog_sqrt
```

We can also only change the link function for the precision parameter.

For instance, we will set the link function for the precision parameter to be *sqrt*:
```{r}
fit_priming_sqrt <- bbreg(agreement ~ priming + eliciting| priming, data = WT, 
                          link.precision = "sqrt")
fit_priming_sqrt
```

## Getting fitted values

If we want to get the fitted values, we can use the *fitted* method on the
fitted *bbreg* object. We can have four types of fitted values: the *response*,
which provides the fitted means, the *link* which provides the the fitted 
linear predictors for the means, the *precision* which provides the fitted
precisions and *variance* which provides the fitted variances for the response
variables.

For example, consider the **fit** object in the first example. 
We can obtain the fitted means by simply calling the *fitted* method:
```{r}
fitted_means <- fitted(fit)
head(fitted_means)
```
and by passing the argument *type = "variance"*, we get:
```{r}
fitted_variances <- fitted(fit, type = "variance")
head(fitted_variances)
```

## Making predictions

To get predictions, we use the *predict* method on the fitted *bbreg* object.
We can pass an additional argument *newdata* equal to a data frame containing the values
of the covariates to get predictions based on these covariates. If one does not pass
the *newdata* argument, this functions returns the *fitted* values (that is, if
the *newdata* parameter is *NULL*, then this
method is identical to the *fitted* method).

Let us create a data frame containing new covariates for the fitted model **fit_priming_eliciting**:
```{r}
new_data_example <- data.frame(priming = c(0,0,1), eliciting = c(0,1,1))
```

Now, let us obtain the corresponding predicted mean values and predicted variances for the response variables:
```{r}
predict(fit_priming_eliciting, newdata = new_data_example)
predict(fit_priming_eliciting, newdata = new_data_example, type = "variance")
```

It is interesting to observe that the above fitted model has priming and eliciting as covariates for both
the mean and the precision parameters. We do not have to worry about that, the *bbreg* handles it automatically.

The same *new_data_example* can also be used for the fitted model **fit**, which does not have precision
covariates (*priming* and *eliciting* are covariates for the mean):
```{r}
predict(fit, newdata = new_data_example)
predict(fit, newdata = new_data_example, type = "variance")
```

In short, the data frame to be passed as *newdata* must contain all the mean and precision covariates as columns,
the *bbreg* package handles the rest.

Lastly, observe that if we do not pass the argument *newdata*, the *predict* method coincides with the *fitted* method:

```{r}
predict_without_newdata <- predict(fit)
identical(predict_without_newdata, fitted_means)
```

## Creating simulated envelopes

To create simulated envelopes for a *bbreg* object, we just need to set the number of random draws to be made
at the *envelope* argument of the *bbreg* function:
```{r}
fit_envelope <- bbreg(agreement ~ priming + eliciting, envelope = 300, data = WT)
```

Observe that the *bbreg* function (through the usage of the *pbapply* package) provides a nice progress bar 
with an estimate of the remaining time to complete the simulation.

The *summary* method provides the percentage of data within the bands formed by the simulated envelopes. This can
also be used as a criterion to select between the bessel and beta regressions.

```{r}
summary(fit_envelope)
```

The fitted *bbreg* object with simulated envelopes can be used to produce a quantile-quantile plot with simulated
envelopes. 

## Diagnostic plots for *bbreg* objects

The *bbreg* package comes with a *plot* method that currently has four kinds of plots implemented:
Residuals vs. Index; Q-Q Plot (if the fit contains simulated envelopes, the plot will be with the simulated envelopes); Fitted means vs. Response and Residuals vs. Fitted means.

To produce all four plots, one may simply apply the *plot* method to the *bbreg* object:

```{r}
plot(fit)
```

Notice that since there were no simulated envelopes for the **fit** object, the Q-Q plot does
not contain simulated envelopes. Nevertheless, it contains a line connecting the first and third
quartiles of the normal distribution (one can remove this line by setting the argument
*qqline* to *FALSE*).

To choose a specific plot, one can set the parameter *which* to a vector containing the
numbers of the plots to be shown. The plot numbers are given by: Plot 1: 
Residuals vs. Index; Plot 2: Q-Q Plot (if the fit contains simulated envelopes,
the plot will be with the simulated envelopes); Plot 3: Fitted means vs. Response;
Plot 4: Residuals vs. Fitted means.

For instance, if we only want to plot the Q-Q plot, 
we set the argument *which* to 2:

```{r}
plot(fit_envelope, which = 2)
```

Notice that since the **fit_envelope** object contains simulated envelopes, the
Q-Q plot is displayed with simulated envelopes.

Finally, to avoid being asked between plots, just set the *ask* parameter to *FALSE*:

```{r}
plot(fit, which = c(1,4), ask = FALSE)
```