.Rhistory

#---------------------------------------------------------------------
# Chunk 8
all_musical_indv_theory_model <- lmer(scale(rating) ~  scale(lerdhal_tension) +
scale(parncut_roughness) + scale(semitone_voice_move) +
as.character(rootmotion) + scale(number_tendency_tones4_6) +
scale(i_have_had_formal_training_in_music_theory_for_years) +
(1+scale( i_have_had_formal_training_in_music_theory_for_years)|participant),
data = model_data)
r.squaredGLMM(all_musical_indv_theory_model)
# Chunk 9
(re.effects <- plot_model(all_musical_indv_theory_model, type = "est", show.values = TRUE,colors = c("black","black"),value.offset = .25))
levels(re.effects$data$term) <- c("Music Theory Training", "^4, ^#4, b^6, ^6" , "Descending-Fifths Root Motion",
"Semitonal Voice Leading","Chord Roughness","Chord Distance")
re.effects +
labs(title = "Musical Feature and Musical Training Model", subtitle = "Fixed Effects") +
theme_minimal() -> model_3_figure
model_3_figure
ggsave(plot = model_3_figure,
filename = "img/Figures/Figure_Model_Three_fixed_effects_.png",
height = 6, width = 9, units = "in", dpi = 300)
# Chunk 10
(analysis_a_plot_1 + descriptive_theory_plot )/ (rm_musical_plot + model_3_figure) +
plot_annotation( tag_levels = "A") &
theme(plot.tag.position = c(0, 1),
plot.tag = element_text(size = 15, hjust = 0, vjust = 0)) -> compound_plot
compound_plot
ggsave(filename = "img/Figures/Figure2.tiff", height = 9, width = 13, units = "in", dpi = 300)
ggsave(filename = "img/Figures/Figure2.eps", height = 9, width = 13, units = "in", dpi = 300)
# ALL MODELS HERE
tab_model(chord_cat_rm, musical_training_rm, all_musical_indv_theory_model)
# Chunk 11
# Re-Import Tables
# Get only rating data
attraction_ratings <- response_data %>%
dplyr::select(I:`viio7/V`, group, participant)
# Demographic Data
gmsi <- read_csv("tidy_tables/gmsi_table.csv")
gmsi <- janitor::clean_names(gmsi)
gmsi <- gmsi %>%
rename(experimental_group = group,
participant = participant_no)
gmsi %>%
dplyr::select(participant, i_have_had_formal_training_in_music_theory_for_years) -> gmsi_2
attraction_ratings <- attraction_ratings %>%
left_join(gmsi_2)
chord_stimuli_data <- read_csv("tidy_tables/chord_experiment_stimuli_data.csv")
chord_stimuli_data <- chord_stimuli_data %>%
rename(chord = chord_symbol,
chord_family = group)
# Chunk 12
# Make Exploratory Plot
pca_rec <- recipe(~., data = attraction_ratings ) %>%
update_role(group, participant, i_have_had_formal_training_in_music_theory_for_years, new_role = "id") %>%
step_normalize(all_predictors()) %>%
step_pca(all_predictors())
pca_prep <- prep(pca_rec)
pca_prep
tidied_pca <- tidy(pca_prep, 2)
chord_stimuli_data_table <- chord_stimuli_data %>%
rename(terms = chord)
tidied_pca %>%
left_join(chord_stimuli_data_table) -> tidied_pca2
# Chunk 13
# Extract PCA Variance
sdev <- pca_prep$steps[[2]]$res$sdev
percent_variation <- sdev^2 / sum(sdev^2)
sum(percent_variation[1:3])
sum(percent_variation[1])
sum(percent_variation[2])
sum(percent_variation[3])
tibble(
component = unique(tidied_pca$component),
percent_var = percent_variation ## use cumsum() to find cumulative, if you prefer
) %>%
mutate(component = fct_inorder(component)) %>%
ggplot(aes(component, percent_var)) +
geom_col() +
scale_y_continuous(labels = scales::percent_format()) +
labs(x = NULL, y = "Percent Variance",
title = "Cumulative Variance Explained") +
theme_minimal() +
coord_flip() -> var_explained_plot
var_explained_plot
ggsave(filename = "img/Figures/Figure_Variation_Explained.png",
height = 6, width = 9, dpi = 300, plot = var_explained_plot)
# Chunk 14
# PLOTS HERE NOT IN MANUSCRIPT
# EXCLUDE
tidied_pca2 %>%
filter(component %in% paste0("PC", 1:5)) %>%
mutate(component = fct_inorder(component)) %>%
ggplot(aes(value, terms, fill = chord_family)) +
geom_col() +
facet_wrap(~component, nrow = 1) +
labs(y = "Chord", x = "Loading", title = "First Five Principal Components") -> pc15
pc15
ggsave(filename = "img/Figures/Figure_PC1to5.png", height = 6, width = 9, dpi = 300, plot = pc15)
# Chunk 15
tidied_pca2 %>%
filter(component %in% paste0("PC", 1)) %>%
mutate(component = fct_inorder(component)) %>%
mutate(value = round(value, 3)) %>%
mutate(chord_family = str_replace_all(chord_family, "_"," ")) %>%
mutate(chord_family = str_to_title(chord_family)) %>%
mutate(chord_family = str_replace_all(chord_family,"Pd","PDs")) %>%
mutate(chord_family = str_replace_all(chord_family,"Not","Non")) %>%
mutate(chord_family = str_replace_all(chord_family,"Common","Diatonic")) %>%
mutate(chord_family = factor(chord_family, levels = c("Chromatic PDs","Diatonic PDs","Bridge Chords", "Non PDs"))) %>%
ggplot(aes(x = value, y = reorder(terms, abs(value)))) +
geom_point(aes(color = chord_family), size = 3) + # swapped from shape
theme_minimal() +
scale_color_viridis(option = "D", discrete = TRUE) +
theme(legend.position="bottom") +
labs(y = "Chord", x = "Loading",
title = "Response Strategy: Experimental Prompt", shape = "Chord Category",
color = "Chord Category") -> pc1
# chromatic - purple
# diantoinc - blue
# bridge -green
# nonpd - yellow
pc1
ggsave(filename = "img/Figures/PC1.png",
height = 9, width = 5, dpi = 300, plot = pc1)
# Chunk 16
tidied_pca2 %>%
filter(component %in% paste0("PC", 2)) %>%
mutate(component = fct_inorder(component)) %>%
ggplot(aes(x = value, y = reorder(terms, abs(value)), fill = chord_family)) +
geom_col() +
theme_minimal() +
scale_fill_manual(labels = c("Bridge Chords","Chromatic PDs","Common PDs","Non PDs"),
values = c("green","black","blue","red")) +
labs(y = "Chord", title = "PC2: Sensitivity to Training", fill = "Chord Family") -> pc2
pc2
ggsave(filename = "img/Figures/PC2.png", height = 6, width = 9, dpi = 300, plot = pc2)
# Chunk 17
# PANEL 2
juice(pca_prep) %>%
mutate(group = str_replace_all(group, "upperclass","Juniors/Seniors")) %>%
mutate(group = str_replace_all(group, "freshman","freshmen")) %>%
mutate(group = str_to_title(group)) %>%
mutate(`Experimental Group` = str_to_title(group)) %>%
mutate(`Experimental Group` = factor(`Experimental Group`, levels = c("Juniors/Seniors","Freshmen","Prolific"))) %>%
ggplot(aes(PC1, PC2,
# shape = `Experimental Group`,
color = `Experimental Group`)) +
geom_point(alpha = 1) +
scale_shape_manual(values=c( 8, 13, 4))+
scale_color_viridis(discrete = TRUE, option = "A", begin = .5, end = .9, direction = -1) +
theme_minimal() +
theme(legend.position="bottom") +
labs(title = "Scores of Items from PCA Loadings",
x = "PC1: Predominant Attraction",
y = "PC2: Sensitivity to Music Theory Training",
shape = "Participant Group") -> pc1vpc2
pc1vpc2
ggsave(filename = "img/Figures/PC1vsPC2.png", height = 6, width = 9, dpi = 300, plot = pc1vpc2)
# Chunk 18
# Make Panel
library(patchwork)
(pc1 + (pc1vpc2)) + plot_annotation(tag_levels = "A") -> pca_panel
library(cowplot)
cowplot::plot_grid(pc1, pc1vpc2,labels = "AUTO", label_size = 24) -> pca_panel_2
ggsave(plot = pca_panel, filename = "img/Figures/Figure3.png", height = 9, width = 12, dpi = 300)
ggsave(plot = pca_panel_2, filename = "img/Figures/Figure3cow.png", height = 9, width = 12, dpi = 300)
ggsave(plot = pca_panel, filename = "img/Figures/Figure3.tiff", height = 9, width = 12, dpi = 300)
ggsave(plot = pca_panel_2, filename = "img/Figures/Figure3cow.tiff", height = 9, width = 12, dpi = 300)
ggsave(plot = pca_panel, filename = "img/Figures/Figure3.eps", height = 9, width = 12, dpi = 300)
ggsave(plot = pca_panel_2, filename = "img/Figures/Figure3cow.eps", height = 9, width = 12, dpi = 300)
# Chunk 19
# Old Exploratory Plots, Not in Manuscript
juice(pca_prep) %>%
ggplot(aes(PC1, PC3, color = group)) +
geom_point(alpha = 0.7, size = 2) +
labs() -> pc1vpc3
ggsave(filename = "img/unsupervised/pc1vpc3.png", height = 6, width = 9, dpi = 300, plot = pc1vpc3)
juice(pca_prep) %>%
ggplot(aes(PC2, PC3, color = group)) +
geom_point(alpha = 0.7, size = 2) +
labs() -> pc2vpc3
ggsave(filename = "img/unsupervised/pc2vpc3.png", height = 6, width = 9, dpi = 300, plot = pc2vpc3)
# Chunk 20
# Correlation with general average
# tidied_pca2 %>%
#   filter(component == "PC1") %>%
#   dplyr::select(terms, value) %>%
#   rename(chord = terms) -> pca_terms
#
#
# feature_model_data %>%
#   dplyr::select(mean_chord_ratings, chord) %>%
#   left_join(pca_terms)-> cor_pca
#
# cor_pca %>%
#   ggplot(aes(x = value, y = mean_chord_ratings, label = chord)) +
#   geom_point() +
#   geom_text(nudge_y = .1) +
#   geom_smooth()
#
# cor.test(cor_pca$mean_chord_ratings, cor_pca$value)
# Chunk 21
chord_stimuli_data <- read_csv("tidy_tables/chord_experiment_stimuli_data.csv")
chord_stimuli_data <- chord_stimuli_data %>%
rename(chord = chord_symbol,
chord_family = group)
# Chord Features
chord_features <- read_csv("tidy_tables/chord_features.csv")
chord_features <- janitor::clean_names(chord_features)
chord_features <- chord_features %>%
separate(col = stimuli, into = c("chord_id","temp"), "\\.") %>%
separate(col = temp, into = c("chord","notes"), sep = "\\:") %>%
mutate(chord = str_trim(chord)) %>%
mutate(notes = str_trim(notes))
chord_features %>%
select(parncut_roughness:number_4_b6) -> cluster_features
rownames(cluster_features) <- chord_features$chord
cluster_features %>% scale() -> scaled_cluster_features
dist_cluster <- dist(scaled_cluster_features)
dist_cluster_object <- hclust(dist_cluster, method = "complete")
par(mfrow=c(1,1))
jpeg("img/Figures/Figure4.jpeg", height = 500, width = 900)
plot(dist_cluster_object , xlab="", sub="", cex=1.5, main = "HCA of Chord Features")
dev.off()
plot(dist_cluster_object , xlab="", sub="", cex=1.5, main = "HCA of Chord Features")
# ADD COLOR =====================================================================
#install.packages("dendextend")
library(dendextend)
# dont judge me for doing this manually
#scale_color_manual(values = c("#55C667FF","#440154FF","#39568CFF","#FDE725FF"))
# chromatic - purple
# diantoinc - blue
# bridge -green
# nonpd - yellow
chord_stimuli_data %>%
select(chord, chord_family) %>%
print(n = 100)
#jpeg("img/Figures/Figure4_color.jpeg", height = 400, width = 800)
#tiff("img/Figures/Figure4_color.tiff", height = 400, width = 800)
par(mar=c(1,1,1,1))
tiff("img/Figures/Figure4_color_rescale.tiff", width = 10, height = 6, units = 'in', res = 300)
dist_cluster_object %>%
as.dendrogram %>%
set("labels_col", c(
"#440154FF",  # Purple
"#440154FF",  # Purple
"#440154FF",  # Purple
"#440154FF",  # Purple
"#39568CFF",  # Blue
"#440154FF",  # Purple
"#440154FF",  # Purple
"#440154FF",  # Purple
"#39568CFF",  # Blue
"#39568CFF",  # Blue
"#FDE725FF",  # #55C667FF --> #FDE725FF
"#39568CFF",  # Blue
"#39568CFF",  # Blue
"#39568CFF",  # Blue
"#440154FF",  # Purple (N6)
"#FDE725FF",  # #55C667FF (vii07)--> #FDE725FF
"#39568CFF",  # Blue
"#39568CFF",  # Blue
"#FDE725FF",  # #55C667FF (i7)--> #FDE725FF
"#55C667FF",  # --  > #55C667FF
"#55C667FF",  # --> #55C667FF
"#55C667FF",  # --> #55C667FF
"#55C667FF",  # --> #55C667FF
"#55C667FF",  # --> #55C667FF
"#55C667FF",  # --> #55C667FF
"#FDE725FF",  # #55C667FF (I7)--> #FDE725FF
"#FDE725FF",  # #55C667FF --> #FDE725FF
"#FDE725FF",  # #55C667FF --> #FDE725FF
"#55C667FF",  # --> #55C667FF
"#55C667FF",  # --> #55C667FF
"#55C667FF"   # --> #55C667FF
)) %>%
set("labels_cex", 1.5) %>%
plot(main = "HCA Chord Features")
legend("topright",
legend = c("Chromatic PDs" , "Diatonic PDs" , "Bridge Chords" , "Non PDs"),
col = c(
"#440154FF",  # Purple
"#39568CFF",  # Blue
"#55C667FF",
"#FDE725FF" ), # Yellow
pch = c(20), bty = "n", cex = 1.5 ,
text.col = "black", horiz = FALSE, inset = c(0, 0))
dev.off()
#=================================================================================
# Linage does not seem to affect too much in terms of general categories
dist_cluster_object_complete <- hclust(dist_cluster,method = "complete")
dist_cluster_object_average <- hclust(dist_cluster,method = "average")
dist_cluster_object_single <- hclust(dist_cluster, method = "single")
par(mfrow=c(1,3))
plot(dist_cluster_object_complete ,main="Complete Linkage", xlab="", sub="",
cex=.9)
plot(dist_cluster_object_average , main="Average Linkage", xlab="", sub="",
cex=.9)
plot(dist_cluster_object_single, main="Single Linkage", xlab="", sub="",
cex=.9)
library(readr)
library(stringr)
library(ggplot2)
library(tidyr)
library(GGally)
library(MuMIn)
library(janitor)
library(lme4)
library(lmerTest)
library(patchwork)
library(gt)
library(psych)
library(Hmisc)
library(forcats)
library(multcomp)
library(ggpubr)
library(dplyr)
library(sjPlot)
library(sjlabelled)
library(forcats)
library(janitor)
library(tidymodels)
library(viridis)
# Importing Data
## Experimental Data
response_data <- read_csv("tidy_tables/response_data.csv")
chord_stimuli_data <- read_csv("tidy_tables/chord_experiment_stimuli_data.csv")
chord_stimuli_data <- chord_stimuli_data %>%
rename(chord = chord_symbol,
chord_family = group)
# Chord Features
chord_features <- read_csv("tidy_tables/chord_features.csv")
chord_features <- janitor::clean_names(chord_features)
chord_features <- chord_features %>%
separate(col = stimuli, into = c("chord_id","temp"), "\\.") %>%
separate(col = temp, into = c("chord","notes"), sep = "\\:") %>%
mutate(chord = str_trim(chord)) %>%
mutate(notes = str_trim(notes))
chord_features[chord_features$chord=="Gr6",]$chord <- "Ger6"
# Match Encoding of iiø7 and viiø7
chord_features[chord_features$chord=="iiØ7",]$chord <- "iiø7"
chord_features[chord_features$chord=="viiØ7",]$chord <- "viiø7"
# Demographic Data
gmsi <- read_csv("tidy_tables/gmsi_table.csv")
gmsi <- janitor::clean_names(gmsi)
gmsi <- gmsi %>%
rename(experimental_group = group,
participant = participant_no)
#View(chord_features)
# Build Big Data
df <- response_data %>%
pivot_longer(cols = I:`viio7/V`,names_to = "chord",values_to = "rating") %>%
select(participant, chord, rating, block, group, `How attracted`) %>%
left_join(chord_stimuli_data) %>%
left_join(gmsi) %>%
left_join(chord_features)
# Table 2
df %>%
select(chord_family, chord, parncut_roughness:number_tendency_tones4_6) %>%
distinct() -> table_2
#write_csv(table_2,"img/table2.csv")
```
## Descriptives
```{r}
df %>%
select(participant, group, i_have_had_formal_training_in_music_theory_for_years) %>%
distinct() %>%
group_by(group) %>%
summarise(mean_theory = mean(i_have_had_formal_training_in_music_theory_for_years),
sd_theory = sd(i_have_had_formal_training_in_music_theory_for_years))
df %>%
select(participant, group, i_have_had_formal_training_in_music_theory_for_years) %>%
distinct() %>%
group_by(group) %>%
summarise(mean_theory = mean(i_have_had_formal_training_in_music_theory_for_years),
sd_theory = sd(i_have_had_formal_training_in_music_theory_for_years))
## Results
### Supervised
First we check for test-restest reliability then collapse data based on R1F  =  0.92 Reliability of average of all items for one  time (Random time effects) and RkF  =  0.96 Reliability of average of all items and both times (Fixed time effects).
```{r}
# Check test-retest with psych package
response_data %>%
select(block,I:`viio7/V`) %>%
rename(time = block) %>%
data.frame() %>%
testRetest()
# Reduce df to variables for modeling
model_data <- df %>%
select(rating,
chord,
chord_family,
lerdhal_tension,
parncut_roughness,
semitone_voice_move,
rootmotion,
number_tendency_tones4_6,
participant,
block,
i_have_had_formal_training_in_music_theory_for_years,
f3_musical_training,
experimental_group) %>%
mutate(experimental_group = str_replace_all(experimental_group, pattern = "freshman",replacement = "freshmen")) %>%
mutate(experimental_group = str_replace_all(experimental_group, pattern = "upperclass",replacement = "Juniors/Seniors")) %>%
mutate(experimental_group = factor(experimental_group, c("Juniors/Seniors","freshmen", "prolific"))) %>%
mutate(chord_family_f = factor("non_pd","bridge","common_pd","chromatic"))
# Test retest, OK, collapse to one average rating per participant
model_data %>%
group_by(chord, participant) %>%
summarise(mean_rating = mean(rating)) -> summary_ratings
model_data %>%
select(chord, chord_family, lerdhal_tension:participant, i_have_had_formal_training_in_music_theory_for_years, f3_musical_training, experimental_group) %>%
left_join(summary_ratings) -> collapsed_model
# Check test-retest with psych package
response_data %>%
select(block,I:`viio7/V`) %>%
rename(time = block) %>%
data.frame() %>%
testRetest()
# Reduce df to variables for modeling
model_data <- df %>%
select(rating,
chord,
chord_family,
lerdhal_tension,
parncut_roughness,
semitone_voice_move,
rootmotion,
number_tendency_tones4_6,
participant,
block,
i_have_had_formal_training_in_music_theory_for_years,
f3_musical_training,
experimental_group) %>%
mutate(experimental_group = str_replace_all(experimental_group, pattern = "freshman",replacement = "freshmen")) %>%
mutate(experimental_group = str_replace_all(experimental_group, pattern = "upperclass",replacement = "Juniors/Seniors")) %>%
mutate(experimental_group = factor(experimental_group, c("Juniors/Seniors","freshmen", "prolific"))) %>%
mutate(chord_family_f = factor("non_pd","bridge","common_pd","chromatic"))
# Test retest, OK, collapse to one average rating per participant
model_data %>%
group_by(chord, participant) %>%
summarise(mean_rating = mean(rating)) -> summary_ratings
model_data %>%
select(chord, chord_family, lerdhal_tension:participant, i_have_had_formal_training_in_music_theory_for_years, f3_musical_training, experimental_group) %>%
left_join(summary_ratings) -> collapsed_model
## Regression Analyses
We now run each hypothesis as linear, mixed effects regression model.
We note here that [ANOVA models are a special case of regression](https://lindeloev.github.io/tests-as-linear/) and that running a linear mixed effects model with a random effect of participant is [analgous to running a within subjects repeated measures ANOVA in R](https://m-clark.github.io/docs/mixedModels/anovamixed.html).
### Chord Category Model
```{r}
# CHORD CATEGORY
# RM-ANOVA as LINEAR MIXED EFFECTS
chord_cat_rm <- lmer(mean_rating ~ chord_family + (1|participant), data = collapsed_model)
# Model Summary Output
summary(chord_cat_rm)
anova(chord_cat_rm)
r.squaredGLMM(chord_cat_rm)
options(scipen = 999)
summary(chord_cat_rm, ddf = "Satterthwaite")
# Post Hoc Comparisions
glht(chord_cat_rm, mcp(chord_family="Tukey")) %>% summary()
# CHORD CATEGORY
# RM-ANOVA as LINEAR MIXED EFFECTS
chord_cat_rm <- lmer(mean_rating ~ chord_family + (1|participant), data = collapsed_model)
# Model Summary Output
summary(chord_cat_rm)
anova(chord_cat_rm)
r.squaredGLMM(chord_cat_rm)
options(scipen = 999)
summary(chord_cat_rm, ddf = "Satterthwaite")
# Post Hoc Comparisions
glht(chord_cat_rm, mcp(chord_family="Tukey")) %>% summary()
#### Figure For Chord Category Model
```{r}
means_a <- tibble(means = c(4.39518,4.39518+.41810,4.39518+.52897,4.39518+.03238),
std_errors = c(.09757,.05051,.05051, .05229),
labz = c("Bridge Chords","Chromatic PDs","Diatonic PDs","Non PDs"))
means_a %>%
mutate(labz = factor(labz,
levels = c("Chromatic PDs","Diatonic PDs","Bridge Chords", "Non PDs"))) %>%
ggplot(aes(y = means, x = labz, fill = labz)) +
geom_bar(stat = "identity") +
#  geom_errorbar(aes(ymin = means - std_errors, ymax = means+std_errors)) +
labs(x = "Chord Category", y = "Adjusted Means", title = "Chord Category Model", fill = "Chord Category") +
scale_fill_viridis(discrete = TRUE) +
geom_bracket(inherit.aes = FALSE, xmin = c("Bridge Chords","Bridge Chords", "Non PDs","Non PDs"),
xmax = c("Chromatic PDs","Diatonic PDs","Chromatic PDs","Diatonic PDs"),
y.position = c(7, 6, 9 , 8),
label = c("***"),
tip.length = 0.1,
) +
theme_minimal() -> analysis_a_plot_1
analysis_a_plot_1
# Same analysis, run with base R calculations
#chord_cat_rm_anova <- aov(mean_rating ~ chord_family + Error(participant), data = collapsed_model)
#summary(chord_cat_rm_anova)
```
### Musical Training Model
```{r}
# MUSICAL TRAINING
# Run as Linear Mixed Effects
# Main Effects: Theory Training and Chord Category
# Interaction
# Participant as Random Effects, within subject
musical_training_rm <- lmer(mean_rating ~ chord_family*i_have_had_formal_training_in_music_theory_for_years  + (1|participant),
data = collapsed_model)
summary(musical_training_rm)
anova(musical_training_rm)
r.squaredGLMM(musical_training_rm)
summary(musical_training_rm, ddf = "Satterthwaite")
anova(musical_training_rm)
r.squaredGLMM(musical_training_rm)
anova(musical_training_rm)
sessionInfo()