06-poly-lca-plot.Rmd

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, 
                      warning = FALSE,
                      message = FALSE) #Here, I have made it so that when you knit your .rmd, warnings and messages will not show up in the html markdown. 
library(extrafont)
loadfonts()
```


# Item Probability Plots

Polytomous LCA deals with variables that have more than two categories, such as survey questions with responses like `never`, `sometimes`, and `always`. The workflow of a polytomous LCA model is similar to that of an LCA model with binary indicators. However, polytomous LCA captures more complex response patterns, which can make interpretation more complex. The following code demonstrates an example, along with a visualization of the model.

------------------------------------------------------------------------

## Example: Elections

"Two sets of six questions with four responses each, asking respondents’ opinions of how well various traits describe presidential candidates Al Gore and George W. Bush. Also potential covariates vote choice, age, education, gender, and party ID. Source: The National Election Studies (2000)." (poLCA, 2016) [See documentation here](https://cran.r-project.org/web/packages/poLCA/poLCA.pdf)

Two sets of six questions with four responses each, asking respondents’ opinions of how well various traits describe presidential candidates Al Gore and George W. Bush. In the election data set, respondents to the 2000 American National Election Study public opinion poll were asked to evaluate how well a series of traits—moral, caring, knowledgeable, good leader, dishonest, and intelligent—described presidential candidates Al Gore and George W. Bush. Each question had four possible choices: (1) extremely well; (2) quite well; (3) not too well; and (4) not well at all.

------------------------------------------------------------------------

Load packages

```{r, cache = FALSE}

library(poLCA)
library(tidyverse)
library(janitor)
library(gt)
library(MplusAutomation)
library(here)
library(glue)
```

------------------------------------------------------------------------

## Prepare Data

```{r, eval=TRUE}
data(election)

# Detaching packages that mask the dpylr functions 
detach(package:poLCA, unload = TRUE)
detach(package:MASS, unload = TRUE)

df_election <-  election %>% 
  clean_names() %>% 
  select(moralb:intelb) %>% 
  mutate(across(everything(), 
                ~ as.factor(as.numeric(gsub("\\D", "", .))), 
                .names = "{.col}1")) 

# Quick summary
summary(df_election)
```

------------------------------------------------------------------------

## Descriptive Statistics

```{r}
ds <- df_election %>% 
  pivot_longer(moralb1:intelb1, names_to = "variable") %>% 
  count(variable, value) %>%  # Count occurrences of each value for each variable
  group_by(variable) %>%
  mutate(prop = n / sum(n)) %>% 
  arrange(desc(variable))

# Create the table
prop_table <- ds %>% 
  gt() %>% 
  tab_header(title = md("**Descriptive Summary**")) %>%
  cols_label(
    variable = "Variable",
    n = md("*N*"),
    prop = md("Proportion")
  ) %>%
  fmt_number(c("n", "prop"), decimals = 2) %>%  # Format both n and prop columns
  cols_align(
    align = "center",
    columns = c(prop, n)
  ) 

# View the table
prop_table

# Save as a Word doc
#gtsave(prop_table, here("figures", "prop_table.docx"))
```


------------------------------------------------------------------------

## Item Probability Plot 

The functions `poLCA_stacked` and `poLCA_grouped` create visualizations of class probabilities for LCA with polytomous indicators. Each function takes the following arguments:

- **`model_name`**: The LCA model read into R using the `readModels` function from the `MplusAutomation` package.
- **`category_labels`**: A character vector of category labels for the response options (e.g., survey answers). 

Note: Double check that the labels are in the correct order!


```{r, fig.width= 18, fig.height= 11}
source(here("functions","poLCA_plot.R"))

# Read in models
output_election <- readModels(here("poLCA"), filefilter = "election", quiet = TRUE)


poLCA_stacked(output_election$c5_election.out, category_labels = c("1" = "1: Extremely well", 
                                                                   "2" = "2: Quite Well", 
                                                                   "3" = "3: Not Too Well", 
                                                                   "4" = "4: Not Well at All"))
```


Alternative plot

```{r,  fig.width= 18, fig.height= 11}
poLCA_grouped(output_election$c5_election.out, category_labels = c("1" = "1: Extremely well", 
                                                                   "2" = "2: Quite Well", 
                                                                   "3" = "3: Not Too Well", 
                                                                   "4" = "4: Not Well at All"))
```

------------------------------------------------------------------------

## APA-formatted Plot

```{r}
# Model 
model <- output_election$c5_election.out


# Title
title <- "2000 Descriptions of Presidential Candidate George W. Bush; Item Probabilities by Class"

# Item names
item_labels <- c("CARESB1" = "Caring",
                 "DISHONB1" = "Dishonest",
                 "INTELB1" = "Intelligent",
                 "KNOWB1" = "Knowledgeable",
                 "LEADB1" = "Good Leader",
                 "MORALB1" = "Moral")

# Item Category
category_labels <- c("1" = "1: Extremely well", 
                     "2" = "2: Quite Well", 
                     "3" = "3: Not Too Well", 
                     "4" = "4: Not Well at All")

# Class labels
class_labels <- c("1" = "Poor Decsription (9.95%)",
                  "2" = "Mostly Poor Description (22.40%)",
                  "3" = "In-Between (24.06%)",
                  "4" = "Mostly Well-Described But Not Intelligent (28.29%)",
                  "5" = "Well-Described But Not Intelligent (15.30%)")


#### END EDIT ####
```

### Extract data needed for plotting

```{r}
# Extract data needed for plotting 
plot_data <- data.frame(model$parameters$probability.scale) %>%
  dplyr::select(est, LatentClass, param, category) %>%
  mutate(
    items = factor(param, labels = item_labels),
    class = factor(LatentClass, labels = class_labels),
    cat = factor(category, labels = category_labels)
  ) %>% 
  mutate(class = factor(class, levels = rev(levels(factor(class))))) 
```

------------------------------------------------------------------------

### Final grouped bar plot

```{r, fig.height=9, fig.width=15}

## Plot data
plot_data %>%
  ggplot(aes(
    x = items,
    y = est,
    fill = cat,
    group = cat
  )) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(label = sub("^0\\.", ".", sprintf("%.2f", est))), 
            position = position_dodge(width = 0.9), 
            vjust = -0.5, size = 3) +
  facet_wrap(~ class) + 
  ylim(0, 1) +
  scale_x_discrete(
    "",
    labels = function(x)
      str_wrap(x, width = 10)
  ) +
  labs(title = "Figure 1",
       subtitle = title,
       y = "Probability") +
  theme_bw(12) +
  scale_fill_grey(start = 0.8, end = 0.2) + # Gives different shades
  theme(
    text = element_text(family = "sans", size = 12),
    legend.text = element_text(family = "sans", size = 12, color = "black"),
    legend.title = element_blank(),
    legend.position = "bottom",
    legend.justification = "center", 
    axis.text.x = element_text(vjust = 1),
    plot.subtitle = element_text(face = "italic", size = 15),
    plot.title = element_text(size = 15),
    strip.background = element_rect(fill = "grey90", color = "black", size = 1),
    strip.text = element_text(size = 12)
  ) 
```

------------------------------------------------------------------------

Save figure:

```{r, eval = FALSE}
ggsave(here("figures", "APA_plot1.png"), dpi="retina", bg = "white", height=9, width=15, units="in")
```

------------------------------------------------------------------------

### Alternative

```{r, fig.height=10, fig.width=17}
## Plot data
plot_data %>%
  ggplot(aes(
    x = items,
    y = est,
    fill = cat,
    group = cat
  )) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(label = sub("^0\\.", ".", sprintf("%.2f", est))), 
            position = position_dodge(width = 0.9), 
            vjust = -0.5, size = 3) +
  facet_wrap(~ class) + 
  ylim(0, 1) +
  scale_x_discrete(
    "",
    labels = function(x)
      str_wrap(x, width = 10)
  ) +
  labs(title = "Figure 1",
       subtitle = title,
       y = "Probability") +
  theme_cowplot(12) +
  scale_fill_grey(start = 0.8, end = 0.2) + # Gives different shades
  theme(
    text = element_text(family = "sans", size = 12),
    legend.text = element_text(family = "sans", size = 12, color = "black"),
    legend.title = element_blank(),
    legend.position = "bottom",
    legend.justification = "center", 
    axis.text.x = element_text(vjust = 1),
    plot.subtitle = element_text(face = "italic", size = 15),
    plot.title = element_text(size = 15),
    strip.background = element_rect(fill = "grey90", color = "black", size = 1),
    strip.text = element_text(size = 12)
  )
```

------------------------------------------------------------------------

Save figure:

```{r, eval = FALSE}
ggsave(here("figures", "APA_plot2.png"), dpi="retina", bg = "white", height=10, width=17, units="in")
```

<div style="text-align: center;"><img src="images/ucsb_logo.png" width="75%" /></div>