---
title: "A Multimodal Symphony: Data analysis"
date: 2025-02-12
authors:
  - name: Matteo Spanio
    email: spanio@dei.unipd.it
    id: ms
    url: https://matteospanio.github.io
    orcid: 0000-0002-2436-7208
    affiliation:
      - name: University of Padova
        city: Padua
        url: https://www.unipd.it
    roles: writing
    corresponding: true
  - name: Massimiliano Zampini
    email: massimiliano.zampini@unitn.it
    id: mz
    url: https://webapps.unitn.it/du/en/Persona/PER0000123/Curriculum
    orcid: 0000-0001-5950-7365
    affiliation:
      - name: University of Trento
        city: Rovereto
        url: https://www.unitn.it/
    roles: writing
    corresponding: false
  - name: Antonio Rodà
    email: roda@dei.unipd.it
    id: ar
    url: https://www.dei.unipd.it/~roda/
    orcid: 0000-0001-9921-0590
    affiliation:
      - name: University of Padova
        city: Padua
        url: https://www.unipd.it
    roles: writing
    corresponding: false
abstract: >
  In recent decades, neuroscientific and psychological research has traced direct relationships between taste and auditory perceptions. This article explores multimodal generative models capable of converting taste information into music, building on this foundational research. We provide a brief review of the state of the art in this field, highlighting key findings and methodologies. Additionally, we present an experiment in which we fine-tuned a generative music model (MusicGEN) to generate music based on detailed taste descriptions provided for each musical piece. The results are promising: the fine-tuned model produces music that more coherently reflects the input taste descriptions compared to the non-fine-tuned model. This study represents a significant step towards understanding and developing embodied interactions between AI, sound, and taste, opening new possibilities in the field of generative AI.
keywords:
  - Generative AI
  - Crossmodal correspondences
  - Statistical analysis
  - Synesthesia
license: "CC BY"
citation:
    type: software
    container-title: "Zenodo"
    doi: 10.5281/zenodo.14864652
    url: https://zenodo.org/records/14864652
funding: "We thank SoundFood s.r.l. for the economic support of this research."
bibliography: references.bib
csl: ieee.csl
---

## Demographic Analysis

The data used for this study were collected through an [online survey](https://www.psytoolkit.org/c/3.4.6/survey?s=YmxcY) via PsyToolkit's web platform [@stoet_psytoolkit_2017]. Inclusion criteria were that participants had to be over eighteen years old and have access to a device capable of playing audio files.

```{r}
#| include: false

library(tidyverse)
library(cluster)
library(psych)

source("src/load_data.R")
```

Of `{r} nrow(data_raw)` people reached, only `{r} nrow(data_pinfo)` completed the questionnaire (`{r} sum(data_pinfo$sex == 'Male')` males, `{r} sum(data_pinfo$sex == 'Female')` females, `{r} sum(data_pinfo$sex == 'Other')` other, `{r} sum(data_pinfo$sex == 'Not specified')` prefer not to say, see @fig-demographic-info-1), null or partial answers were considered as withdrawal from the questionnaire, therefore only complete answers were taken into consideration for the following analysis.

:::: {#fig-demographic-info layout="[[1, 1],[1, 1]]"}
::: {#fig-demographic-info-1}
```{r}
library(ggplot2)

# Assuming your dataframe is named 'data_pinfo' and the gender column is 'sex'
gender_distribution <- data_pinfo |>
  count(sex) |>
  mutate(percentage = n / sum(n) * 100)

ggplot(gender_distribution, aes(x = "", y = percentage, fill = sex)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar("y") +
  theme_void()
```
Gender distribution of participants.
:::

::: {#fig-demographic-info-2}
```{r}
ethnicity_distribution <- data_pinfo |>
  count(ethnicity) |>
  mutate(percentage = n / sum(n) * 100)

ggplot(ethnicity_distribution, aes(x = "", y = percentage, fill = ethnicity)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar("y") +
  theme_void()
```
Distribution of participants’ ethnic backgrounds.
:::

::: {#fig-demographic-info-3}
```{r}
hearing_experience_distribution <- data_pinfo |>
  count(hearing_experience) |>
  mutate(percentage = n / sum(n) * 100)

ggplot(hearing_experience_distribution, aes(x = "", y = percentage, fill = hearing_experience)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar("y") +
  theme_void()
```
Distribution of participants’ hearing experiences.
:::

::: {#fig-demographic-info-4}
```{r}
eating_experience_distribution <- data_pinfo |>
  count(eating_experience) |>
  mutate(percentage = n / sum(n) * 100)

ggplot(eating_experience_distribution, aes(x = "", y = percentage, fill = eating_experience)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar("y") +
  theme_void()
```
Distribution of participants’ eating experiences.
:::
Demographic characteristics of the study’s participants.
::::

Overall the reached population has a mean age of `{r} round(mean(data_pinfo$age, na.rm = TRUE))` years, whit a maximum of `{r} max(data_pinfo$age, na.rm = TRUE)` and a minimum of `{r} min(data_pinfo$age, na.rm = TRUE)`. The mean time spent on the survey was `{r} round(mean(data_pinfo$time, na.rm = TRUE))` minutes, with a standard deviation of `{r} round(sd(data_pinfo$time, na.rm = TRUE), 1)`. Along with age, gender and execution time also data about musical and eating experience have been collected: @fig-demographic-info-2 displays the ethnicity distribution of the population, the majority of participants recognize themself as *White/European American*, participants have an almost equally distributed experience in listening to music (see @fig-demographic-info-3), while just one participant recognized himself as an experienced eater and the major part of the sample population declared to be not-experienced in tasting food (@fig-demographic-info-4).

## Model Preference Analysis

The first task in the survey involved participants listening to two audio clips, each corresponding to a taste description chosen randomly from *sweet*, *sour*, *bitter*, and *salty*. The goal was to determine which audio sample best matched the given taste description. The two clips were generated by different versions of the MusicGEN model [@copet2024simplecontrollablemusicgeneration]: a fine-tuned version and the original [^1], base model, released by Meta. Participants were asked to express their preference for the first or second clip by moving a slider ranging from 0 to 10, where 0 indicated a strong preference for the first clip and 10 indicated a strong preference for the second, @fig-task1 shows the survey's first question interface.

[^1]: A full description of the model and its finetuning process is available at the publication related to this analysis.

![Screenshot of the survey's first task interface.](assets/img/model_preference_woh.png){#fig-task1 .lightbox}

To ensure randomization and avoid any bias, the taste descriptions and audio clips were presented in a random order. In the analysis the scores are normalized as follows: scores from 0 to 4 are interpreted as a preference for the base model, scores between 6 and 10 indicate a preference for the fine-tuned model, and scores of 5 are treated as neutral.

The underlying research question guiding this task was to assess if the fine tuned model output better matches taste descriptions than the sounds generated by the base model.

### Data Visualization

The distribution of scores across all participants is presented in @fig-model-pref-1. The histogram and density plot show the spread of scores, allowing us to visually assess the preference for one model over the other. The base model and fine-tuned model preferences are expected to manifest as peaks around the lower and higher end of the score range, respectively.

:::: {#fig-model-pref layout="[[1, 1], [1]]"}

::: {#fig-model-pref-1}
```{r}
ggplot(data_task1, aes(x = score)) +
  geom_histogram(binwidth = 1, fill = "skyblue", color = "black", alpha = 0.7) +
  geom_density(aes(y = after_stat(count)), fill = "lightgreen", alpha = 0.4) +
  labs(x = "Score", y = "Count") +
  theme_minimal()
```

Overall evalutation of the models.
:::

::: {#fig-model-pref-2}

```{r}
ggplot(data_task1, aes(x = prompt, y = score, fill = prompt)) +
  geom_violin() +
  geom_hline(yintercept = 5, color = "red", linetype = "dashed", linewidth = 1) +
  labs(
    x = "Taste",
    y = "Score"
  ) +
  scale_fill_brewer(palette = "Set2") +
  theme_minimal() +
  theme(legend.position = "none")
```

Score boxplot by taste category.
:::

::: {#fig-model-pref-3}
```{r}
# Creazione dei grafici per ciascun valore della variabile 'prompt'
ggplot(data_task1, aes(x = score)) +
  geom_histogram(binwidth = 1, fill = "skyblue", color = "black", alpha = 0.7) +
  geom_density(aes(y = after_stat(count)), fill = "lightgreen", alpha = 0.4) +
  labs(x = "Score", y = "Count") +
  theme_minimal() +
  facet_wrap(~ prompt, scales = "free_y")
```
Score distribution by taste category.

:::

Score distribution between the two models.
::::

@fig-model-pref-2 goes further by breaking down the preferences based on the taste category described in the audio sample. Each taste category is visualized using boxplots, where the median score for each taste can be assessed. This enables us to examine whether the model preference varies depending on the taste label, with the red dashed line at a score of 5 acting as the neutral threshold.




### Hypothesis Testing

Next, we assess whether the average score for the audio samples significantly differs from a neutral score of 5, under the assumption that the preference for the fine-tuned model should be greater than this threshold. To do this, we need to verify whether the data follows a normal distribution.
In order to assess normality of data we applied both visual and computational methods, then firstly a Q-Q plot was generated to visually inspect the normality of the score distribution, see @fig-qqplot. The resulting plot shows deviations from the expected straight line, indicating that the scores do not follow a normal distribution.

```{r}
#| label: fig-qqplot
#| fig-cap: Q-Q plot to assess the normality distribution of the collected data.
library(qqplotr)

ggplot(data_task1, aes(sample = score)) +
  stat_qq_band(bandType="pointwise", fill="#8DA0CB",alpha=0.4) +
  stat_qq_line(color = "#8DA0CB") +
  stat_qq_point() +
  labs(
    x = "Theoretical Quantiles",
    y = "Sample Quantiles"
  ) +
  theme_minimal()
```


```{r}
# Test Shapiro-Wilk
shapiro.pvalue <- (shapiro.test(data_task1$score))$p.value
```
In addition the Shapiro-Wilk test confirmed that the data significantly deviate from a normal distribution (with a resulting $p$-value equals to $`{r} signif(shapiro.pvalue, digits = 4)`$). Therefore, we decide to apply the non-parametric Wilcoxon signed-rank test to see if the models preference score expressed by participants has a mean major than the null preference (score equals to five), more formally we are testing the hypothesis:

$$
H_0: \mu = 5 \quad H_1: \mu > 5
$$
where $H_0$ means that there is no preference between the two models while $H_1$ means that the fine-tuned model is preferred over the other one with a $95%$ confidence interval.

```{r}
alpha <- 0.05
hyp_mu <- 5

wilcox.pvalue <- (wilcox.test(data_task1$score, mu = hyp_mu, alternative = "greater", conf.level = 1-alpha))$p.value
```

The result of the Wilcoxon test shows a $p$-value of $`{r} signif(wilcox.pvalue, digits = 4)`$, which is less than $0.05$, indicating that we can reject the null hypothesis and conclude that the median score is indeed significantly greater than $5$. This supports the hypothesis that the participants prefer the fine-tuned model overall.

### Post-Hoc Analysis by Taste

While the Wilcoxon test confirms that the overall preference goes to the fine-tuned model the boxplots reveal a variation of the score by taste, to confirm the variation we perform separate Wilcoxon tests for each taste group (*sweet*, *sour*, *bitter*, *salty*). We use a Bonferroni correction to adjust for the multiple comparisons and control the family-wise error rate. The results of the post-hoc tests are shown below, in @tbl-wilcoxon-taste.

```{r}
#| label: tbl-wilcoxon-taste
#| tbl-cap: Results of the Wilcoxon test performed on data filtered by taste.

library(knitr)

levels_taste <- unique(data_task1$prompt)

# List to stpre the results
test_results <- list()

# Execute the Wicoxon test for each category
for (level in levels_taste) {
  # Filter data
  subset_data <- data_task1[data_task1$prompt == level, ]

  test_result <- wilcox.test(subset_data$score, mu = 5, alternative = "greater", conf.level = 1 - alpha)

  # Add the category result
  test_results[[level]] <- test_result
}

# extract the pvalues
p_values <- sapply(test_results, function(res) res$p.value)

# apply the Bonferroni correction
adjusted_p_values <- p.adjust(p_values, method = "bonferroni")

# Visualizza i risultati con i p-value corretti
result_table <- data.frame(
  Taste = names(test_results),
  p.value = p_values,
  adjusted.p.value = adjusted_p_values
) |> select(c("p.value", "adjusted.p.value"))

kable(result_table)
```

The analysis reveals that the fine-tuned model was significantly preferred for the sweet taste category, with an adjusted $p$-value of $`{r} adjusted_p_values[4]`$, well below the conventional threshold of $0.05$. This suggests a strong alignment between the musical outputs and participants’ expectations of sweetness.
Conversely, the bitter and sour categories also exhibited significant preferences, with adjusted $p$-values of $`{r} adjusted_p_values[1]`$ and $`{r} adjusted_p_values[3]`$, respectively. However, these results, while statistically significant, indicate a less robust preference compared to the sweet category.
Notably, the salty taste group did not demonstrate a significant preference for the fine-tuned model, as indicated by an adjusted $p$-value near to 1. This lack of significance suggests that the model’s performance may not align with participants’ expectations for salty flavors, warranting further investigation into the underlying factors influencing this outcome.

```{r}
subset_data <- data_task1[data_task1$prompt == "salty", ]
salty_base <- wilcox.test(subset_data$score, mu = 5, alternative = "less", conf.level =1 - alpha)$p.value
```

Since the finetuned model did not show to perform well on the salty group, we performed a Wilcoxon test to test if its mean is significally lower than the tie situation ($H_0: \mu_{\text{salty}} = 5, H_1: \mu_{\text{salty}} < 5$) the result is that the first model is statistically preferred over the finetuned variant according to the Wilcoxon test with a $p$-value equal to $`{r} salty_base`$ without Bonferroni correction[^2].

[^2]: The Bonferroni correction has not been applied due to the non indepenent nature of the test, in fact the test has been performed after the results of the previous Wilcoxon test, which, instead, was testing independently 4 groups.

## Recognisability of Tastes

In the second task of the survey, participants were asked to listen to musical pieces generated exclusively by the fine-tuned model to better investigate the intrinsic qualities carried by the generated music. Following each listening session, participants were required to quantify the flavors they perceived in the music using a graduated scale, ranging from 1 to 5 (where 1 means *not at all* and 5 means *very much*), for each of the four primary taste categories: salty, sweet, bitter, and sour. Unlike the first task, there were no imposed labels for specific flavors, allowing participants the freedom to associate values with each taste based on their personal interpretations of the musical experience.
Additionally, to enrich the assessment, participants had to evaluate their emotional responses to the music by rating various non-gustatory parameters, including happiness, sadness, anger, disgust, fear, surprise, hot, and cold. @fig-task2 displays the web interface used to collect participants' responses.

![Screenshot of the survey's second task interface.](assets/img/music_synesthesia.png){#fig-task2 .lightbox}

The underlying research questions guiding this task were:

1. Can the music generated by the model induce sensory-gustatory responses?
2. What correlations exist between music and taste?
3. Which emotions mediate these sensory responses?

### ANOVA test

To address the first research question, we performed an Analysis of Variance (ANOVA) to evaluate whether the participants' ratings of stimuli varied according to distinct stimulus characteristics. The dependent variable was the value assigned by participants to each stimulus, while the independent variables included stimulus-related factors. The dataset was filtered to include only participants identifying as *Male* or *Female*, excluding participants classified as *Professional Eaters* due to insufficient representation of this category. This preprocessing step ensured robust and meaningful comparisons between groups.

```{r}
#| label: tbl-anova-value
#| tbl-cap: Results of the ANOVA test.
library(broom)

df <- data_aov |>
  filter(sex == "Male" | sex == "Female") |>  # too few people recognize themself as non male/female
  filter(eating_experience != "Professional") # the only professional eater is filtered out

res <- aov(value ~ prompt * adjective + hearing_experience + eating_experience + sex, data = df)
kable(tidy(res), digits = c(0, 0, 3, 3, 3, 3),
      col.names = c("", "Df", "Sum Sq", "Mean Sq", "F value", "Pr(>F)"))
```

The results of the ANOVA are summarized in @tbl-anova-value, which presents the degrees of freedom (Df), sum of squares (Sum Sq), mean squares (Mean Sq), F-statistics, and $p$-values for each factor and interaction.

Prior to interpreting the results, the homoskedasticity assumption was assessed by examining the residuals. A Shapiro-Wilk test indicated evidence of heteroskedasticity ($p < 0.05$). Despite this violation, the ANOVA analysis proceeded, following recommendations from prior research [@glass_1972_consequences_of_failure; @harwell_1992_meta-analytic_methods; @lix_1996_presence_of_variance] suggesting that ANOVA is robust to deviations from normality under moderate violations, particularly with large sample sizes such as the one in this study.
The results show that different prompts and adjectives lead to significantly different adjectives quantified by participants, similarly the adjectives used influence the participants’ feelings. Furthermore the significant interaction effect implies that the effect of one variable depends on the level of the other. In other words, the way a prompt influences feelings may vary depending on the adjective used. This result highlights that the participants to the survey deliberately operated consistent choices while evaluating the stimuli.

### Post-Hoc Analysis

To further explore the results of the ANOVA, we conducted Tukey’s Honest Significant Difference (HSD) test. This post-hoc analysis is particularly useful for identifying which specific group means are significantly different from each other after finding a significant overall effect in the ANOVA. Given that our analysis revealed significant main effects for prompt, adjective, their interaction and the hearing experience, it is essential to determine the nature of these differences.

The Tukey test compares all possible pairs of group means while controlling for the family-wise error rate, thus providing a robust method for multiple comparisons. This is crucial in our context, as we aim to understand how different prompts, adjectives and hearing experience levels influence participants’ evaluations of the stimuli.
Upon executing the Tukey test, we examined the adjusted $p$-values for each comparison. The results indicated several significant differences between specific combinations of prompts and adjectives, as summarized in the table below:

```{r}
res <- aov(value ~ prompt + adjective + hearing_experience, data = data_aov)
tukey <- TukeyHSD(res)

# For each Tukey result, filter where p adj is less than 0.05
significant_prompts <- as.data.frame(tukey$prompt) |> filter(`p adj` < 0.05)
significant_adjectives <- as.data.frame(tukey$adjective) |> filter(`p adj` < 0.05)
significant_hearing_experience <- as.data.frame(tukey$hearing_experience) |> filter(`p adj` < 0.05)
```

```{r}
#| tbl-cap: Tukey test between different prompts with a $p$-value lower 0.05.
kable(significant_prompts)
```

```{r}
#| tbl-cap: Tukey test between different adjectives with a $p$-value lower 0.05.
kable(significant_adjectives)
```

```{r}
#| tbl-cap: Tukey test between different hearing experience groups with a $p$-value lower 0.05.
kable(significant_hearing_experience)
```

These significant comparisons illuminate the subtleties in participants’ responses to different stimuli. Certain prompt-adjective combinations elicited stronger emotional reactions than others, indicating that the interaction between prompts and adjectives significantly shapes participants’ perceptions.
Notably, some combinations yielded adjusted $p$-values below the conventional threshold of $0.05$, signifying statistically significant differences. This finding reinforces the ANOVA results, confirming that the presentation of prompts and adjectives can meaningfully impact emotional responses.

### Interaction between *prompt* and *adjective*

The ANOVA analysis evidenced also a significant interaction between prompt and adjectives used to evaluate the sounds, as we know, the design of the experiment fixed the prompt before generating the audio files, therefore adjectives has to be intended as dependent variables while the prompts are independent; in other words, participants assigned different values to the adjectives to the sounds on the basis of their generation prompt. This interaction can be seen in @fig-heatmap. In particular @fig-heatmap-taste shows the mean value assigned to each taste adjective by their prompt, we can clearly seen the major diagonal emerge by the $4 \times 4$ matrix, this means that, the mean value assigned to the adjective that matches the prompt of each sound is the highest. The rest of the interaction between adjectives and prompt can be seen in @fig-heatmap-emotions , a deeper analysis of emotional aspect assigned to the sounds is presented in the next section.

:::: {#fig-heatmap layout="[[1, 1]]"}
::: {#fig-heatmap-taste}
```{r}
library(reshape2)

# Calcolo della media per ciascun prompt
heatmap_data <- data_task2 |>
  group_by(prompt) |>
  summarise(
    sour = mean(sour, na.rm = TRUE),
    sweet = mean(sweet, na.rm = TRUE),
    bitter = mean(bitter, na.rm = TRUE),
    salty = mean(salty, na.rm = TRUE)
  )

# Riorganizzazione per la heatmap
heatmap_matrix <- heatmap_data |>
  pivot_longer(cols = -prompt, names_to = "adjective", values_to = "mean_value") |>
  pivot_wider(names_from = adjective, values_from = mean_value) |>
  column_to_rownames("prompt") |>
  as.matrix()

axis_order <- c("bitter", "salty", "sour", "sweet")

# Creazione della heatmap con ggplot2
heatmap_plot <- heatmap_data |>
  pivot_longer(cols = -prompt, names_to = "adjective", values_to = "mean_value") |>
  mutate(
    prompt = factor(prompt, levels = rev(axis_order)),
    adjective = factor(adjective, levels = axis_order)
  ) |>
  ggplot(aes(x = adjective, y = prompt, fill = mean_value)) +
  geom_tile() +
  geom_text(aes(label = round(mean_value, 2)), color = "black", size = 4) +
  scale_fill_gradient(low = "white", high = "blue") +
  labs(
    x = "Adjective",
    y = "Prompt",
    fill = "Mean value"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Mostrare il grafico
print(heatmap_plot)
```
Heatmap of perceived taste in correspondence of the intended one.
:::

::: {#fig-heatmap-emotions}
```{r}
# Calcolo della media per ciascun prompt
heatmap_data <- data_task2 |>
  group_by(prompt) |>
  summarise(
    happy = mean(happy, na.rm = TRUE),
    sad = mean(sad, na.rm = TRUE),
    anger = mean(anger, na.rm = TRUE),
    disgust = mean(disgust, na.rm = TRUE),
    fear = mean(fear, na.rm = TRUE),
    surprise = mean(surprise, na.rm = TRUE),
    hot = mean(hot, na.rm = TRUE),
    cold = mean(cold, na.rm = TRUE)
  )

# Riorganizzazione per la heatmap
heatmap_matrix <- heatmap_data |>
  pivot_longer(cols = -prompt, names_to = "adjective", values_to = "mean_value") |>
  pivot_wider(names_from = adjective, values_from = mean_value) |>
  column_to_rownames("prompt") |>
  as.matrix()

# Definire l'ordine degli assi
axis_order <- c("sad", "anger", "disgust", "fear", "surprise", "happy", "hot", "cold")

# Creazione della heatmap con ggplot2
heatmap_plot <- heatmap_data |>
  pivot_longer(cols = -prompt, names_to = "adjective", values_to = "mean_value") |>
  mutate(
    prompt = factor(prompt, levels = rev(unique(prompt))),
    adjective = factor(adjective, levels = axis_order)
  ) |>
  ggplot(aes(x = adjective, y = prompt, fill = mean_value)) +
  geom_tile() +
  geom_text(aes(label = round(mean_value, 2)), color = "black", size = 4) +
  scale_fill_gradient(low = "white", high = "blue") +
  labs(
    x = "Adjective",
    y = "Prompt",
    fill = "Mean Value"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Mostrare il grafico
print(heatmap_plot)
```
Heatmap of perceived emotional response in correspondence of the suggested taste.
:::
Heatmaps
::::

### Factorial analysis

```{r}
# Standardize the data excluding the 'taste' column
data <- data_task2 |> select(-prompt)
```

To explore the underlying relationships between sensory qualities and emotional states, we conducted a factor analysis on the standardized data, excluding the ‘taste’ column. The initial step involved calculating the correlation matrix, which revealed notable relationships among the adjectives. Specifically, we observed that negative emotions were positively correlated, while the pair happy-sad exhibited a negative correlation. Furthermore, sweetness demonstrated a strong correlation with happiness and warmth, whereas bitterness was associated with anger and fear. Sourness, on the other hand, was evidently correlated with disgust and fear. Other variables did not show strong correlations at first glance, prompting us to proceed with the factor analysis.

```{r}
#| label: fig-corr-matrix
#| fig-cap: Correlation matrix
library("ggcorrplot")

cor_matrix <- cor(data)

ggcorrplot(cor_matrix,
           lab = TRUE,
           lab_size = 2.9,
           type="lower")
```
The correlation matrix is illustrated in @fig-corr-matrix, showcasing these relationships clearly.

```{r message=FALSE}
#| label: fig-scree-plots
#| fig-cap: Parallel Analysis Scree Plots
#| message: false
# Perform factor analysis
# We'll start by estimating the number of factors, then perform the analysis
fa_parallel <- fa.parallel(data, fm = "ml", fa = "fa", main="")
num_factors <- fa_parallel$nfact # Optimal number of factors from parallel analysis
```

To determine the optimal number of factors for our analysis, we employed parallel analysis, which indicated an optimal number of `{r} num_factors` factors. This estimation serves as a foundation for our subsequent factor analysis.
Following this, we executed the factor analysis using the identified number of factors, applying an oblique rotation (oblimin) to allow for potential correlations among the factors. The results of the factor analysis, including the factor loadings, are presented in the output. The factor loadings indicate how strongly each variable contributes to the identified factors, providing insights into the underlying structure of the data.

```{r}
#| label: fig-factor-analysis
#| fig-cap: Factor analysis graph.

# Run the factor analysis with the identified number of factors
factor_analysis <- fa(data, nfactors = num_factors, rotate = "oblimin", fm = "ml", scores=TRUE)

fa.diagram(factor_analysis, main="")
```

```{r}
# Display factor loadings
print(factor_analysis$loadings)
```

To visualize the relationships among the factors, we generated biplots for various factor combinations. The biplots, shown in the subsequent figures, illustrate the distribution of variables across the identified factors, highlighting the clustering of adjectives associated with similar emotional states.

:::: {#fig-fa-biplot layout="[[1, 1], [1, 1]]"}

::: {#fig-fa-biplot-1}
```{r}
biplot.psych(
  factor_analysis,
  choose=c(1,2),
  smoother = TRUE,
  main = NULL,
  )
```
:::

::: {#fig-fa-biplot-2}
```{r}
biplot.psych(
  factor_analysis,
  choose=c(3,4),
  main = NULL,
  smoother = TRUE)
```
:::

::: {#fig-fa-biplot-3}
```{r}
biplot(
  factor_analysis,
  choose=c(1,3),
  main = NULL,
  smoother=TRUE
)
```
:::

::: {#fig-fa-biplot-4}
```{r}
biplot(
  factor_analysis,
  main = NULL,
  smoother = TRUE,
  choose=c(1,4),
)
```
:::

Loadings biplot over the four factors.
::::

```{r}
library(dplyr)
library(knitr)
library(kableExtra)

# Creiamo un dataframe con gli scores e il prompt
scores_df <- cbind(data_task2["prompt"], as.data.frame(factor_analysis$scores))

# Rinominare le colonne dei fattori se necessario
colnames(scores_df)[2:5] <- paste0("Factor", 1:4)

# Calcoliamo la media e la deviazione standard per ogni fattore, raggruppati per prompt
summary_table <- scores_df |>
  group_by(prompt) |>
  summarise(across(starts_with("Factor"),
                   list(mu = mean, sigma = sd),
                   .names = "{.col}_{.fn}"))

# Riorganizziamo i dati in un formato tabellare con le coppie (mu, sigma)
summary_table_formatted <- summary_table |>
  mutate(across(contains("_mu"),
                ~ sprintf("$\\mu=%.2f$, $\\sigma=%.2f$", .x, get(sub("_mu", "_sigma", cur_column())))))

# Selezioniamo solo le colonne formattate
summary_table_kable <- summary_table_formatted |>
  select(prompt, contains("_mu"))

# Rinominiamo le colonne per renderle più leggibili
colnames(summary_table_kable) <- c("Prompt", paste0("Fattore ", 1:4))

# Creiamo la tabella con kableExtra
kable(summary_table_kable)
```

Lastly, we performed a multi-factor analysis to further explore the dimensions of negativity and positivity within the data. The results of this analysis are depicted in the multi-factor diagram, which categorizes the factors into two overarching themes: Negativity and Positivity.

::: {#fig-multilevel-fa}
```{r}
fam <- fa.multi(data, nfact2=2, nfactors=4)
fa.multi.diagram(fam, f2labels = c("Negativity", "Positivity"), main="")
```

Hierarchical (multilevel) factors' structure.
:::

## References

::: {#refs}
:::
