---
title: "How to use ugcflss"
output:
  rmarkdown::html_vignette:
    toc: true
    fig_width: 6.5
    fig_height: 3.5
vignette: >
  %\VignetteIndexEntry{ugcflss_tutorial}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

## Setup

### Load the package

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

### Data analysis example

I use hypothetical data in this example. The sum scores are the sum of 10 items, each with response categories coded 1 to 4. Hence, the sum score can range from 10 $(10 \times 1)$ to 40 $(10 \times 4)$. These data were collected from 300 people who belong to 6 groups: A -- F.

```{r create_data, cache = FALSE}
set.seed(12345)
dt <- data.frame(grp = sample(LETTERS[1:6], 300, replace = TRUE))
dt$sum_score <- as.integer(cut(stats::rnorm(
  nrow(dt), -15 + 1.5 * as.integer(factor(dt$grp)), 15
), stats::qnorm(seq(0, 1, length.out = 32), 0, 15))) + 9
summary(dt)
```

We can visualize the data using a boxplot:

```{r data_boxplot, cache = FALSE}
boxplot(sum_score ~ grp, dt, horizontal = TRUE)
```

## Analysis

### Run the model

Once we have a dataset with a grouping variable (`grp` in this example), and a sum score (`sum_score` in this example), we are ready to use the ugcflss model fitting function. It will place a call to RStan so expect the model to run for several seconds. 

When running the model, we need to inform the package about the minimum and maximum response categories on the items that were summed to create the sum score and the number of items we summed, here: 1, 4 and 10 respectively.

```{r fit_model, cache = FALSE}
fitted_model <- ugcflss_fit_model(
  data = dt,
  grouping_variable = "grp", sum_score = "sum_score",
  minimum_item_response = 1, maximum_item_response = 4,
  chains = 1, # you can likely ignore this argument
  number_items = 10
)
```

We see the iteration printout from Stan. Stan sometimes produces warnings about divergent transitions (see `adapt_delta` argument in `?ugcflss_fit_model()`) or maximum treedepth (see `max_treedepth` argument in  `?ugcflss_fit_model()`). Setting both arguments to higher values than their defaults -- `adapt_delta = .99` instead of .9 or `max_treedepth = 15` instead of 10 -- when calling `ugcflss_fit_model(..., adapt_delta = .99, max_treedepth = 15)` can resolve these warnings. We have no such warnings here, so we proceed.

## Model checking via the posterior predictive distribution

### Built-in checks

One way to examine Bayesian models is to simulate data from the models then assess how well data from the models reflects the original data -- this is called posterior predictive checking. To assess how well the simulated data matches the distribution of the original data, try:

```{r ppd_dens, cache = FALSE}
# ppd_samples is the number of datasets we want to simulate
ugcflss_ppd_plot(fitted_model, ppd_samples = 500)
```

where $y$ is the density of the overall data and $y_{rep}$ is the density of the simulated datasets.

You could also do this by group:

```{r ppd_dens_group, cache = FALSE}
ugcflss_ppd_plot(fitted_model, ppd_samples = 500, by_group = TRUE)
```

### Using the bayesplot package

For more detailed assessment, it is possible to directly request simulated datasets:

```{r ppd_req, cache = FALSE}
ppd_draws <- ugcflss_ppd(fitted_model, ppd_samples = 500)
dim(ppd_draws) # should be 500 rows by 300 people
```

The `bayesplot` package has several features for comparing this distribution to the original data. For example, we can compare the observed count for each value to the distribution of counts based on the model:

```{r ppd_bars, cache = FALSE}
# the ppd draws are on the 1--4 scale, transform to sum scores:
ppd_sum_score <- ppd_draws * fitted_model$user_input$n_items
bayesplot::ppc_bars(dt$sum_score, ppd_sum_score)
```

Here, the distribution ($y_{rep}$ error bars) should ideally contain the actual bars.

You could also check specific statistics, with the expectation that the model-implied distribution of the statistic includes the sample statistic:

```{r ppd_stats, cache = FALSE}
# Checking the SD
bayesplot::ppc_stat(dt$sum_score, ppd_sum_score, stat = "sd")
# Or the 75% percentile
q75 <- function(y) quantile(y, .75)
bayesplot::ppc_stat(dt$sum_score, ppd_sum_score, stat = "q75")
```

We now turn to several outputs that help us understand the group differences.

## Simple descriptive statistics

Results are all on the 1 -- 4 scale, i.e. assuming the items were averaged rather than scored. I made this choice because this is the scale I see reported in manuscripts.

### For the overall sample

```{r describe_overall_1, cache = FALSE}
overall_median <- ugcflss_describe(fitted_model, stat = "median")
overall_median
```

The estimate of the median is `r scales::number(overall_median$mean, .01)`. By default, I report 89% quantile intervals i.e. the middle 89% of the distribution of the parameter. So one might say the median was `r scales::number(overall_median$mean, .01)` (SE = `r scales::number(overall_median$sd, .001)`), 89% CrI [`r scales::number(overall_median$q5.5, .01)`, `r scales::number(overall_median$q94.5, .01)`], where CrI stands for credible interval.

Statistics such as `rhat`, `ess_bulk` and `ess_tail` are parameter convergence statistics, i.e. can we trust that the parameter has converged to a final solution? `ess` stands for effective sample size, and we want this statistic in the hundreds. `rhat` under 1.01 or 1.05 (liberal) suggests parameter convergence. One way to make `rhat` smaller and `ess_` larger is to increase the number of iterations when calling `ugcflss_fit_model()`. For example, you can modify the model fitting call: `ugcflss_fit_model(..., warmup = 2000, sampling = 2000)`, both arguments are 750 by default.

Other potential statistics are `stat = "median"`, `stat = "mean"`, `stat = "sd"`, or `stat = "quantile"`. When requesting quantiles, also specify the `tau` argument to set the desired percentile. I might also want conventional 95% intervals:

```{r describe_overall_2, cache = FALSE}
overall_q78 <- ugcflss_describe(
  fitted_model,
  stat = "quantile", tau = .78, interval = .95
)
overall_q78
```

The 78% percentile was `r scales::number(overall_q78$mean, .01)` (SE = `r scales::number(overall_q78$sd, .001)`), 95% CrI [`r scales::number(overall_q78$q2.5, .01)`, `r scales::number(overall_q78$q97.5, .01)`].

For those more familiar with Bayesian analysis, you can request these computed samples by calling `ugcflss_describe(..., return_draws = TRUE)`.

### By group

We can also request any of the statistics above by group:

```{r describe_group_1, cache = FALSE}
ugcflss_describe(fitted_model, stat = "sd", by_group = TRUE)
```

These are the group standard deviations. Out of curiosity, I also show their sample estimates:

```{r describe_group_2, cache = FALSE}
aggregate(I(sum_score / 10) ~ grp, dt, sd)
```

The model-based results from the package produce group standard deviations that are more similar to each other than those based on sample statistics.

As with the entire sample, one can modify the default interval using `interval = .xx` or request draws using `return_draws = TRUE`.

### Plotting the descriptive statistics

To plot the descriptive statistics e.g. the `mean`, call:

```{r describe_plot_1, cache = FALSE}
ugcflss_describe_plot(fitted_model, stat = "mean")
```

## Pairwise comparisons

We can perform pairwise comparisons using:

```{r pairwise_1, cache = FALSE, out.width = "150%"}
pairwise_mean_diff <- ugcflss_pairwise(fitted_model, stat = "mean")
pairwise_mean_diff
```

We see pairwise mean differences, for example: the mean of group B was less than the mean of group A by `r scales::number(pairwise_mean_diff$mean[1], .01)` (SE = `r scales::number(pairwise_mean_diff$sd[1], .001)`), 89% CrI [`r scales::number(pairwise_mean_diff$q5.5[1], .01)`, `r scales::number(pairwise_mean_diff$q94.5[1], .01)`].

The `p_direction` column has the probability of direction. This is the probability that a coefficient whose estimate is positive (negative) is indeed positive (negative). Hence, there was a `r scales::percent(pairwise_mean_diff$p_direction[1], .1)` chance that the `r scales::number(pairwise_mean_diff$mean[1], .01)` coefficient was `r ifelse(pairwise_mean_diff$mean[1] > 0, "positive", "negative")`. This is a Bayesian analog to _p_-values, and quantifies the evidence in favour of a particular direction for an effect.

We can also perform ratio and log-ratio comparisons. For example, here is the pairwise log-ratio comparison of group standard deviations:

```{r pairwise_2, cache = FALSE, out.width = "150%"}
ugcflss_pairwise(fitted_model, stat = "sd", comparison = "log-ratio")
```

Ratios of positive numbers have an assymetric scale. When the numerator is less, the ratio will range between 0 and 1, otherwise, it will range between 1 and infinity. Using a log-ratio comparison transforms 0--1 to to a number less than 0, and numbers greater than 1 to a number above 0 -- thus making the comparison scale symmetric.

`comparison` must be one of `"difference"` (default), `"ratio"`, or `"log-ratio"`.

Again, users who are more familiar with Bayesian statistics can request the computed samples of the pairwise comparisons by calling `ugcflss_pairwise(..., return_draws = TRUE)`.

### Probability of superiority

In addition to the earlier statistics, we can also compute the pairwise probability of superiority (PS). This statistics communicates how likely it is that members of one group will have responses greater than members from another group. For example, a PS of 60% implies that if we randomly recruited 100 pairs of respondents from two groups, A and B, members of group A would have the higher value in 60 of the pairs, while members of group B would have the higher value in 40 pairs. Thus a PS of 50% implies neither group has the higher values.

To perform pairwise comparisons using the PS, call:

```{r eval=FALSE}
ugcflss_pairwise(fitted_model, stat = "ps")
```

The PS is both the `stat` and `comparison` when calling the PS, the value of `comparison` is ignored.

### Plotting pairwise comparisons

It is easier to parse several pairwise comparisons using plots. ugcflss can plot pairwise comparisons:

```{r pairwise_plot_1_prep, cache = FALSE, include = FALSE}
pp_median_comp <- ugcflss_pairwise(
  fitted_model,
  stat = "median", comparison = "ratio"
)
```

```{r pairwise_plot_1, cache = FALSE, fig.height = 2.5}
ugcflss_pairwise_plot(fitted_model, stat = "median", comparison = "ratio")
```

The topleft cell shows the result of median(B) / median(A) -- its value is `r ifelse(pp_median_comp$mean[1] > 1, "above 1", "less than 1")` meaning the median of group B was `r ifelse(pp_median_comp$mean[1] > 1, "greater than", "less than")` the median of group A -- there was a `r scales::percent(pp_median_comp$p_direction[1], 1)` chance this is true. The two values in brackets are the credible intervals.

```{r pairwise_plot_2, cache = FALSE, fig.height = 2.5}
ugcflss_pairwise_plot(fitted_model, stat = "ps")
```

We do not show the probability of direction in the plot for the PS -- it can be confusing.

## Exceedance probability

Sometimes, specific values on scale scores have practical significance, so we may be interested in how likely respondents were to have a score that exceeded some threshold. We can request this statistic:

```{r exceed_overall, cache = FALSE}
exceed_overall <- ugcflss_exceed_all(fitted_model)
exceed_overall
```

Based on the model, `r scales::percent(exceed_overall$md[1], .1)` of respondents had a score greater than `r scales::number(exceed_overall$points[1], .1)`. If you simply count using sample statistics, we get `r scales::percent(exceed_overall$count[1], .1)`. The `lo` and `hi` columns are the 89% CrI by default, though these can be changed using `ugcflss_exceed_all(fitted_model, interval = .xx)`.

We can also request this statistic by group (results not shown):

```{r exceed_group, eval = FALSE}
ugcflss_exceed_group(fitted_model)
```

One can also visualize these results:

```{r exceed_overall_plot, cache = FALSE}
ugcflss_exceed_all_plot(fitted_model)
```

```{r exceed_group_plot, cache = FALSE}
ugcflss_exceed_group_plot(fitted_model)
```

And there are options to show some desired interval (89% by default), e.g.

```{r exceed_overall_plot_int, cache = FALSE}
ugcflss_exceed_all_plot(fitted_model, interval = .99, show_intervals = TRUE)
```

The same is possible by group (results not shown):

```{r exceed_group_plot_int, eval = FALSE}
ugcflss_exceed_group_plot(fitted_model, show_intervals = TRUE)
```