Likert data analyses

This document shows several different analyses of a single Likert-type item. It is by no means exhaustive.

If you want to skip the details, there is a summary of the results at the bottom of the document.

Setup

library(gamlss)
library(gamlss.tr)
library(simpleboot)
library(tidyverse)
library(broom)
library(ggstance)
library(tidybayes)
library(ordinal)
library(betareg)
library(brms)
library(rstanarm)
library(rstan)
library(emmeans)
library(modelr)
library(forcats)

options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
theme_set(theme_light())

Data

Read in likert-type data from the “blinded with science” study. The study has 4 experiments. In each experiment, each participant is either shown a text about a hypothetical drug (“no graph” condition) or the same text with a chart (“graph” condition). The participant is then asked to assess to what extent the drug is effective, on a scale from 1 to 7. We focus on this data here. The study has another dependent variable that’s more important (error of the response to a comprehension question) but it’s not Likert-type data.

Let’s read in the data:

blinded <- read_csv("blinded.csv",
  col_types = cols(
    experiment = col_integer(),
    condition = col_character(),
    effectiveness = col_integer()
  )) %>%
  mutate(
    experiment = paste0("e", experiment),
    condition = fct_recode(condition, `no graph` = "no_graph")
  )
head(blinded)

It looks like this:

blinded %>%
  ggplot(aes(x = ordered(effectiveness))) +
  stat_count() +
  facet_grid(rows = vars(condition), cols = vars(experiment)) +
  xlab("effectiveness")

Analyses

Bootstrap CI

For each experiment, compute the difference between the two means and its 95% bootstrap CI.

First, we’ll create a helper function to compute bootstrap CIs in a tidy way (with broom-compatible output):

tidy_bootstrap_mean_diff <- function(group1, group2, conf.level = 0.95, R = 1000, seed = 98432098) {
  set.seed(seed) # make deterministic
  bootstrap_samples <- two.boot(sample1 = group1, sample2 = group2, FUN = function(x, d) mean(x[d]), R = R)
  bootci <- boot.ci(bootstrap_samples, type = "bca", conf = conf.level)
  data.frame(
    estimate = mean(group1) - mean(group2),
    conf.low = bootci$bca[4],
    conf.high = bootci$bca[5]
  )
}

Now we’ll apply that function within each experiment:

result_boot = blinded %>%
  group_by(experiment) %>%
  do(tidy_bootstrap_mean_diff(
    filter(., condition == "no graph")$effectiveness,
    filter(., condition == "graph")$effectiveness
  )) %>%
  mutate(
    method = "bootstrap",
    estimate_type = "mean"
  )

result_boot

t-test

For each experiment, compute the difference between the two means, its 95% t-based CI and the p-value.

We’ll use the broom::tidy function, which cleans up test and model outputs into a common tidy format:

result_t = blinded %>%
  group_by(experiment) %>%
  do(tidy(t.test(effectiveness ~ fct_rev(condition), data = .))) %>%
  mutate(
    method = "t-test",
    estimate_type = "mean"
  )

result_t

Wilcoxon

Compute the Wilcoxon signed-rank test and p-value. The estimate (and 95% CIs) in this are for the median of the difference between samples (not the difference in the medians).

result_wilcoxon = blinded %>%
  group_by(experiment) %>%
  do(tidy(wilcox.test(effectiveness ~ fct_rev(condition), data = ., conf.int = TRUE, exact = FALSE, correct = FALSE))) %>%
  mutate(
    method = "wilcoxon",
    estimate_type = "Hodges–Lehmann"
  )

result_wilcoxon

Beta regression

For an argument for the use of Beta regression on data like this (in a frequentist mode), see here. There are two tricky bits: (1) the data needs to be rescaled to be in ((0,1)); in the case of a Likert item response in ([1,k]) a natural approach might be to divide the values by (k + 1), in this case, 10. (2) Getting estimates of mean differeneces on the original scale is a pain—for the frequentist version I won’t do this, but for the Bayesian version I will (because I actually know how to do that).

First, we’ll fit a beta regression within each experiment, using default priors:

models_beta_freq = blinded %>%
  group_by(experiment) %>%
  do(model = betareg(I(effectiveness/10) ~ condition, data = .))

Then we’ll get an estimate for the log odds ratio between conditions:

result_beta_freq = models_beta_freq %>% 
  mutate(fits = list(tidy(model, conf.int = TRUE))) %>%
  unnest(fits) %>%
  filter(term == "conditionno graph", component == "mean") %>%
  mutate(
    method = "beta regression",
    estimate_type = "log odds ratio"
  )

result_beta_freq

Bayesian Beta regression

This is a Bayesian formulation of Beta regression. I show the Bayesian approach here because it makes it easier to generate marginal means on the original scale (and intervals on them).

First, we’ll fit a beta regression within each experiment, using default priors:

set.seed(31818) # for reproducibility
models_beta = blinded %>%
  nest(-experiment) %>%
  mutate(model = map(data, ~ stan_betareg(I(effectiveness/10) ~ condition, data = .)))

Then we’ll get an estimate for the predicted marginal means in each condition, and then find the posterior distribution for the difference in those means:

result_beta = models_beta %>% 
  mutate(draws = map(model, fitted_draws, data_grid(blinded, condition))) %>%
  unnest(draws) %>%
  # .value contains draws from marginal means on 0-1 scale,
  # transform back to original scale
  mutate(.value = .value * 10) %>%
  group_by(experiment) %>%
  # get draws from distribution of difference in marginal means
  compare_levels(.value, by = condition) %>%
  median_qi(.value) %>%
  to_broom_names() %>%
  mutate(
    method = "beta reg (Bayes)",
    estimate_type = "mean"
  )

result_beta

Ordinal logistic regression

Another model often applied to Likert-type data is ordinal models. Ordinal models come in plenty of flavors; we’ll use a cumulative link logisitic regression model (also often called ordinal logisitic regression), which is a common variety.

models_ord_freq = blinded %>%
  nest(-experiment) %>%
  mutate(model = map(data, ~ clm(ordered(effectiveness, levels = as.character(1:9)) ~ condition, data = .)))

Then we’ll get an estimate for the log odds ratio between conditions:

result_ord_freq = models_ord_freq %>% 
  mutate(fit = map(model, ~ tidy(pairs(emmeans(., ~ condition), reverse = TRUE), infer = TRUE))) %>%
  unnest(fit) %>%
  rename(
    conf.low = asymp.LCL,
    conf.high = asymp.UCL
  ) %>%
  select(-data, -model) %>%
  mutate(
    method = "ordinal reg",
    estimate_type = "log odds ratio"
  )

result_ord_freq

Bayesian Ordinal logistic regression

We’ll do an ordinal logistic regression in a Bayesian manner because it will make it easier to get estimates of marginal means on the original scale afterwards:

set.seed(221818) # for reproducibility
models_ord = blinded %>%
  nest(-experiment) %>%
  mutate(model = map(data, ~ brm(
    effectiveness ~ condition, data = ., family = cumulative(link = "logit"),
    prior = c(prior(normal(0, 1), class = b), prior(normal(0, 5), class = Intercept)),
    control = list(adapt_delta = .99)
  )))

## Compiling the C++ model

## Start sampling

## Warning: There were 17 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

## Warning: Examine the pairs() plot to diagnose sampling problems

## Compiling the C++ model

## recompiling to avoid crashing R session

## Start sampling

## Compiling the C++ model

## recompiling to avoid crashing R session

## Start sampling

## Compiling the C++ model

## recompiling to avoid crashing R session

## Start sampling

Then pull out the estimates:

result_ord = models_ord %>%
  mutate(draws = map(model, fitted_draws, data_grid(blinded, condition), value = "P(effectiveness|condition)")) %>%
  unnest(draws) %>%
  # calculate marginal means
  group_by(experiment, condition, .draw) %>%
  mutate(effectiveness = as.numeric(as.character(.category))) %>%
  summarise(.value = sum(effectiveness * `P(effectiveness|condition)`)) %>%
  # distribution of difference in marginal means
  group_by(experiment) %>%
  compare_levels(.value, by = condition) %>%
  median_qi(.value) %>%
  to_broom_names() %>%
  mutate(
    method = "ordinal reg (Bayes)",
    estimate_type = "mean"
  )

result_ord

Robust regression

We’ll use a robust, heteroskedastic linear regression: a Student t error distribution instead of Gaussian error distribution, and estimate a different variance parameter for each group. This is essentially Kruschke’s BEST test (a Bayesian, robust, hetereoskedastic alternative to the t-test), but estimated using a frequentist procedure instead of a Bayesian one:

gamlss_fixed = function(..., data) {
  # gamlss has a bug where it does not work if `data` is not in the global environment (!!!!)
  temp_df <<- data
  gamlss(..., data = temp_df)
}

models_robust = blinded %>%
  group_by(experiment) %>%
  do(model = gamlss_fixed(effectiveness ~ condition, sigma.formula = ~ condition, family = "TF", data = .))

## GAMLSS-RS iteration 1: Global Deviance = 431.1339 
## GAMLSS-RS iteration 2: Global Deviance = 431.0836 
## GAMLSS-RS iteration 3: Global Deviance = 431.0688 
## GAMLSS-RS iteration 4: Global Deviance = 431.0642 
## GAMLSS-RS iteration 5: Global Deviance = 431.0627 
## GAMLSS-RS iteration 6: Global Deviance = 431.0623 
## GAMLSS-RS iteration 1: Global Deviance = 613.3834 
## GAMLSS-RS iteration 2: Global Deviance = 613.0218 
## GAMLSS-RS iteration 3: Global Deviance = 612.8633 
## GAMLSS-RS iteration 4: Global Deviance = 612.7906 
## GAMLSS-RS iteration 5: Global Deviance = 612.7573 
## GAMLSS-RS iteration 6: Global Deviance = 612.7422 
## GAMLSS-RS iteration 7: Global Deviance = 612.7352 
## GAMLSS-RS iteration 8: Global Deviance = 612.7321 
## GAMLSS-RS iteration 9: Global Deviance = 612.7307 
## GAMLSS-RS iteration 10: Global Deviance = 612.7301 
## GAMLSS-RS iteration 1: Global Deviance = 650.9131 
## GAMLSS-RS iteration 2: Global Deviance = 650.8548 
## GAMLSS-RS iteration 3: Global Deviance = 650.834 
## GAMLSS-RS iteration 4: Global Deviance = 650.8244 
## GAMLSS-RS iteration 5: Global Deviance = 650.82 
## GAMLSS-RS iteration 6: Global Deviance = 650.818 
## GAMLSS-RS iteration 7: Global Deviance = 650.8171 
## GAMLSS-RS iteration 1: Global Deviance = 566.1012 
## GAMLSS-RS iteration 2: Global Deviance = 566.0954 
## GAMLSS-RS iteration 3: Global Deviance = 566.0953

Then pull out the estimates:

result_robust = models_robust %>%
  mutate(estimates = list(tidy(pairs(emmeans(model, ~ condition), reverse = TRUE), infer = TRUE))) %>%
  unnest(estimates) %>%
  rename(conf.low = asymp.LCL, conf.high = asymp.UCL) %>%
  select(-model) %>%
  mutate(
    method = "robust",
    estimate_type = "mean"
  )

result_robust

Truncated regression

Or how about a truncated normal regression? As before, we’ll also include hetereoskedasticity.

gen.trun(c(min(blinded$effectiveness - 0.5), max(blinded$effectiveness) + 0.5), type = "both")

## A truncated family of distributions from NO has been generated 
##  and saved under the names:  
##  dNOtr pNOtr qNOtr rNOtr NOtr 
## The type of truncation is both 
##  and the truncation parameter is 0.5 9.5

models_trun = blinded %>%
  group_by(experiment) %>%
  do(model = gamlss_fixed(effectiveness ~ condition, sigma.formula = ~ condition, family = "NOtr", data = .))

## GAMLSS-RS iteration 1: Global Deviance = 427.0338 
## GAMLSS-RS iteration 2: Global Deviance = 426.9323 
## GAMLSS-RS iteration 3: Global Deviance = 426.9306 
## GAMLSS-RS iteration 4: Global Deviance = 426.9305 
## GAMLSS-RS iteration 1: Global Deviance = 606.3642 
## GAMLSS-RS iteration 2: Global Deviance = 606.1256 
## GAMLSS-RS iteration 3: Global Deviance = 606.12 
## GAMLSS-RS iteration 4: Global Deviance = 606.1197 
## GAMLSS-RS iteration 1: Global Deviance = 643.449 
## GAMLSS-RS iteration 2: Global Deviance = 643.361 
## GAMLSS-RS iteration 3: Global Deviance = 643.3602 
## GAMLSS-RS iteration 1: Global Deviance = 561.0226 
## GAMLSS-RS iteration 2: Global Deviance = 560.0491 
## GAMLSS-RS iteration 3: Global Deviance = 560.044 
## GAMLSS-RS iteration 4: Global Deviance = 560.0439

Then pull out the estimates:

result_trun = models_trun %>%
  mutate(estimates = list(tidy(pairs(emmeans(model, ~ condition), reverse = TRUE), infer = TRUE))) %>%
  unnest(estimates) %>%
  rename(conf.low = asymp.LCL, conf.high = asymp.UCL) %>%
  select(-model) %>%
  mutate(
    method = "truncated",
    estimate_type = "mean"
  )

result_trun

Summary of results

First, we’ll stick all the result data frames together into a single tidy data frame:

results = bind_rows(
  result_boot,
  result_t,
  result_wilcoxon,
  result_beta,
  result_beta_freq,
  result_ord_freq,
  result_robust,
  result_trun,
  result_ord
) %>%
  ungroup() %>%
  mutate(
    method = fct_reorder(method, estimate, max, .desc = TRUE),
    `estimate type` = estimate_type,
    `95% interval contains 0` = conf.low < 0 & 0 < conf.high
  )

Here are the estimates:

results_plot = results %>%
  ggplot(aes(x = estimate, xmin = conf.low, xmax = conf.high, y = method, shape = `estimate type`, color = `95% interval contains 0`)) +
  geom_pointrangeh() +
  facet_grid(cols = vars(experiment)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_y_discrete(labels = . %>% map(function(x) bquote(underline(.(x))))) +
  theme(axis.text.y = element_text(color = "#4466ff"))

results_plot_original_units = results_plot %+% 
  filter(results, estimate_type != "log odds ratio") +
  xlab("no graph - graph, original units") +
  scale_x_continuous(breaks = c(-1,0,1))

results_plot_original_units

Pretty much however you do it the results come out the same. The above plots show estimates that are in the original units. The Wilcoxon estimates look weird, but it is also estimating a slightly different quantity (median of the differences), so that probably isn’t too surprising.

From two of the models we have also output estimates on a log-odds scale:

results_plot_log_odds = results_plot %+% 
  filter(results, estimate_type == "log odds ratio") +
  guides(shape = FALSE) +
  xlab("graph - no graph, log odds ratio")

results_plot_log_odds

While these are both on the log odds scale, they are measuring log odds of different things, so it is hard to compare the values directly. It does look like if they were scale by their variances they might be similar.

Finally, the p values from all methods (where p values are provided):

p_value_plot = results %>%
  ggplot(aes(x = p.value, y = method, color = p.value < 0.05)) +
  geom_point() +
  facet_grid(cols = vars(experiment)) +
  geom_vline(xintercept = 0.05, linetype = "dashed")

p_value_plot

## Warning: Removed 12 rows containing missing values (geom_point).

Or the same chart on a log scale:

p_value_plot + 
  scale_x_log10()

## Warning: Removed 12 rows containing missing values (geom_point).

And as a table:

results %>%
  arrange(experiment) %>%
  select(experiment, method, statistic, p.value, estimate, estimate_type, conf.low, conf.high)