Analysis of Principal Components x Demographic Factors
In this notebook, principal components are generated from question topic prevalences and disctionary-based estimates of linguistic patterns — and their relationships to various demographic data are explored through statistical tests.
We used a Principle Components Analysis (PCA) in R (‘prcomp’ function) to integrate multiple variables and conduct an analysis for each question and set of responses. PCA is a multivariate statistical technique that extracts important information about the inter-correlations of multiple variables and re-represents that information as new variables called principle components. These new principle components can then be used to show the similarity and graphic relationships between variables as a way of exploring the data. Because PCA condenses multiple variables into fewer dimensions, the new principle components may also be used to integrate new descriptions of the combined variables. These reified principle components may then be used as covariates to estimate relationships with other variables such as demographic factors.
Twenty-seven dictionary-based variables (LIWC, AIC, and Docuscope) were not relevant and were excluded from the analysis, leaving a total of 221 dictionary-based variables in the PCA. This was in addition to the measures of topic prevalence, which each had different numbers of topic variables depending on the question/response type. In order to remove variables with high numbers of 0s (no observations), we set a 0.3 maximum threshold, leaving variable columns with 70% or more observations in the analysis. The number of responses and variables included for each question in the analysis are provided in Table 1. Because of the different scales used by the different dictionaries and topic measures (for example, topic prevalence is 0-1 while AIC is measured on a scale from 1-7), the PCA was based on a correlation matrix. The data were also centered around zero and rotated.
Load libraries and data
library(tidyverse)
library(broom)
library(FactoMineR)
library(magrittr)
library(emmeans)
library(ggplot2)
library(ggExtra)
load('../data/merged_artsengagement.rda')
load('../data/tidy_questions_best.Rda')
source('../scripts/select_helpers.R')
load('../data/pca_labels.rda')
Subset demographic variables
motiv_demos <- df %>% select(contains('sr_motivation')) %>% names
identity_demos <- df %>% select(contains('sr_identity')) %>% names
participation_demos <- df %>% select(all_arts('participation')) %>% select(starts_with('sr')) %>% names
real_demos <- c('ethnic_group', 'sex', 'school', 'parented', 'income', 'artsincollege', 'hstype', 'hssize',
'hslocation', 'hs_arts_freq', 'hs_encouragement', 'hs_required', 'hs_fees',
'so_childhood1', 'so_childhood3', 'so_childhood5', 'sr_participated',
'sr_highestdegreeplanned')
demo_groups <- c('real_demos','motiv_demos','identity_demos','participation_demos')
real_demo_prompts <- c(
'Ethnic Group',
'Sex',
'School',
'Highest Parent\nEducation',
'Family Income',
'Prior to beginning college,\ndid you expect to be involved in the arts during college?',
'High School Type',
'High School Size',
'High School Location',
'How would you rate the frequency\nof your overall involvement\nin arts events and activities during high school?',
'Were you encouraged by your high school to be involved in the arts?',
'Were you required by your high school to be involved in the arts?',
'Was there a fee associated with participation in arts activities at your high school?',
'Were you involved in the arts during your childhood?',
'During your childhood, did your parents sign you up for art classes of any kind?',
'During your childhood, did your parents talk to you about the arts?',
'Arts Participation Senior Year',
'What is the\nhighest academic degree\nthat you intend on obtaining?'
)
Define helper functions
# Used within the mergeQuestion function to merge variable names
stripnames <- function(data, question) {
data <- data %>% select('rownums', starts_with(question))
varnames <- names(data)
newnames <- varnames
for (i in seq(3, length(varnames))) {
newnames[i] <- substr(varnames[i], gregexpr("\\.\\.", varnames[i])[[1]][1]+2, nchar(varnames[i]))
}
names(data) <- newnames
return(data)
}
# Remove irrelavant variables
stripVars <- function(data) {
removenames <- names(data %>% select(matches("AIC..Words|AIC..IC_Differentiation|AIC..IC_Integration|AIC..DIAL_Differentiation|AIC..DIAL_Integration|AIC..ELAB_Differentiation|AIC..ELAB_Integration|docuscope..OralCues|docuscope..DialogCues|docuscope..XWordTokens|docuscope..XPunctuationTokens|docuscope..XTokens|LIWC..WC|LIWC..Sixltr|LIWC..Dic|LIWC..WPS|LIWC..swear|LIWC..netspeak|LIWC..assent|LIWC..nonflu|LIWC..filler|LIWC..AllPunc|LIWC..Period|LIWC..Comma|LIWC..Colon|LIWC..SemiC|LIWC..QMark|LIWC..Exclam|LIWC..Dash|LIWC..Quote|LIWC..Apostro|LIWC..Parenth|LIWC..OtherP")))
data <- data[ , !names(data) %in% removenames]
return(data)
}
# Merges responses between years
mergeQuestion <- function(data, question) {
if (substr(question, 3,3) == '_') {
question <- substr(question, 4, nchar(question))
}
# if you're aware of other prefixes, add them here
yrs <- c(NA,'fl','f1','f2','so','jr','sr','sl')
data$rownums <- seq(nrow(data)) # add rownums so resulting data can be re-integrated
mergeddf <- data.frame()
for (yr in yrs) {
if (question == 'definition' && is.na(yr))
qyr <- question
else if (is.na(yr))
next
else {
qyr <- paste0(yr, '_', question)
}
if (data %>% select(starts_with(qyr)) %>% ncol == 0) { next }
tempdf <- stripnames(data, qyr)
names(tempdf)[2] <- question
mergeddf <- bind_rows(mergeddf, tempdf)
}
mergeddf <- mergeddf[complete.cases(mergeddf),]
mergeddf <- mergeddf %>% stripVars
rownames(mergeddf) <- c()
return(mergeddf)
}
Density Plots
# Filter demographics to those within min & max levels (inclusive) into small_demos variable
min = 2
max = 5
temp <- demodfpcs %>% sapply(function(x)(x %>% levels %>% length)) %>% as.data.frame
temp$demo <- rownames(temp)
temp <- temp %>% filter(. >= min & max >= .)
small_demos <- temp %>% pull(demo)
small_demos
dir_to_save_results = './results'
plot_type = 'density_plots' # plural
output_format = '.png' # with dot
dir.create(dir_to_save_results, showWarnings = F)
# for all stm'd questions
for(question in question_names[-c(4:6,9:14)]) {
dir.create(paste(dir_to_save_results, question, sep = '/'), showWarnings = F)
# data subsetting, PC creation
data_subset <- mergeQuestion(df, question)
data_subset_rows <- data_subset[[1]]
data_subset <- data_subset[-1]
data_subset <- data_subset[-1]
percent <- 0.4
data_subset <- data_subset[ lapply( data_subset, function(x) sum(x==0) / length(x) ) < percent ]
# principal component generation
pd <- prcomp(data_subset, retx= TRUE, center=TRUE, scale=TRUE)
# merge PCs with demographic cols
pcs <- data.frame(rownum = data_subset_rows, pd$x[,seq(5)] %>% as.data.frame)
rownums <- pcs$rownum %>% sort
pcs <- pcs %>% arrange(rownum)
pcs <- pcs[-1]
demodfpcs <- bind_cols(df[rownums,c('key',real_demos)], pcs)
# weighting by # of responses per respondent
resps_by_key <- demodfpcs$key %>% table %>% as.data.frame
names(resps_by_key) <- c('key','n')
resps_by_key %<>% filter(n!=0)
for(r in seq(nrow(demodfpcs))) {
nresps <- resps_by_key %>% filter(key == demodfpcs[r,]$key) %>% .$n
demodfpcs[r, names(select(demodfpcs,matches('PC[12345]')))] <- (1/nresps) *
(demodfpcs[r, names(select(demodfpcs,matches('PC[12345]')))])
}
# separate PCs and demographics
demodfpcs <- demodfpcs %>% select(-'key')
demodf <- demodfpcs %>% select(-matches('PC[12345]'))
pcs <- demodfpcs %>% select(matches('PC[12345]'))
# plotting (by PC)
for(p in seq(pcs)) {
# combine individual PC with demographics
temp <- bind_cols(pcs[p], demodf)
names(temp) <- c("pc", names(temp)[-1])
# for each demographic...
for(d in seq(small_demos)) {
mu <- temp %>% na.omit %>% group_by_at(vars(small_demos[d])) %>%
summarise(grp.median=median(pc))
ggplot(na.omit(demodfpcs), aes_string(x=names(pcs)[p], fill=small_demos[d])) +
geom_density(alpha=0.4) +
xlim(round(na.omit(pcs[[p]]) %>% summary %>% tidy %>% .$minimum)-0.5,
round(na.omit(pcs[[p]]) %>% summary %>% tidy %>% .$maximum)+.5) +
geom_vline(data=mu, aes_string(xintercept='grp.median', color=names(mu)[1]),
linetype="dashed") +
labs(subtitle=paste0('pos-x: ', pca_labels %>%
filter(question_name==question) %>%
select(matches(paste0(p,'_pos'))) %>%
pull %>% as.character,
'\nneg-x: ', pca_labels %>%
filter(question_name==question) %>%
select(matches(paste0(p,'_neg'))) %>%
pull %>% as.character)) +
theme(plot.subtitle = element_text(size=10))
# saving
dir.create(paste(dir_to_save_results, question, names(pcs)[p], sep='/'), showWarnings = F)
dir.create(paste(dir_to_save_results, question, names(pcs)[p], plot_type,sep='/'), showWarnings = F)
ggsave(paste(dir_to_save_results, question, names(pcs)[p], plot_type,
paste0(small_demos[d], output_format),sep='/'),
width = 5, height=4)
}
}
}
Marginal Plots
dir.create('./results', showWarnings = F)
for(question in question_names[-c(4:6,9:14,16)]) {
# directories
dir.create(paste('./results', question, sep = '/'), showWarnings = F)
# data subsetting, PC
data_subset <- mergeQuestion(df, question)
data_subset_rows <- data_subset[[1]]
data_subset <- data_subset[-1]
data_subset <- data_subset[-1]
percent <- 0.4
data_subset <- data_subset[ lapply( data_subset, function(x) sum(x==0) / length(x) ) < percent ]
pd <- prcomp(data_subset, retx= TRUE, center=TRUE, scale=TRUE)
# merge with demographic cols
pcs <- data.frame(rownum = data_subset_rows, pd$x[,seq(5)] %>% as.data.frame)
rownums <- pcs$rownum %>% sort
pcs <- pcs %>% arrange(rownum)
pcs <- pcs[-1]
demodfpcs <- bind_cols(df[rownums,c('key',real_demos)], pcs)
# weighting by # of responses per respondent
resps_by_key <- demodfpcs$key %>% table %>% as.data.frame
names(resps_by_key) <- c('key','n')
resps_by_key %<>% filter(n!=0)
for(r in seq(nrow(demodfpcs))) {
nresps <- resps_by_key %>% filter(key == demodfpcs[r,]$key) %>% .$n
demodfpcs[r, names(select(demodfpcs,matches('PC[12345]')))] <- (1/nresps) *
(demodfpcs[r, names(select(demodfpcs,matches('PC[12345]')))])
}
demodfpcs <- demodfpcs %>% select(-'key')
demodf <- demodfpcs %>% select(-matches('PC[12345]'))
pcs <- demodfpcs %>% select(matches('PC[12345]'))
# analysis
for(p in c(1,3)) {
# directories
cur_dir <- paste('./results', question, paste0('pc', p, '_pc', p+1, '_marginal_plots'), sep = '/')
dir.create(cur_dir, showWarnings = F)
for(d in seq(demodf)) {
temp <- bind_cols(pcs[p], pcs[p+1], demodf[d])
names(temp) <- c(paste0('pc0'),paste0('pc00'),'demo') # for easier ggplot column selection
if((temp$demo %>% table %>% as.data.frame %>% filter(Freq != 0) %>% nrow) == 1) next
pca_idx <- 0
if (p == 1) {
pca_idx <- 2
}
if (p == 3) {
pca_idx <- 6
}
# hexbins
plot <- ggplot(temp %>% na.omit, aes(x=pc0, y=pc00, color=demo)) + geom_point(alpha=4/10) +
labs(color=real_demo_prompts[d]) +
theme(legend.position = "bottom", text = element_text(size=22)) +
xlab(paste0(as.character(filter(pca_labels, question_name == question)[[pca_idx+1]]),
' (-) v. ',
as.character(filter(pca_labels,question_name == question)[[pca_idx]]), ' (+)')) +
ylab(paste0(as.character(filter(pca_labels,question_name == question)[[pca_idx+3]]),
' (-) v. ',
as.character(filter(pca_labels,question_name == question)[[pca_idx+2]]), ' (+)'))
if ((question == 'sr_othergrowth' || question == 'sr_society') && p == 1) {
plot <- plot + xlab(paste0(as.character(filter(pca_labels,question_name == question)$pc1_neg),
' (-) v.\n',
as.character(filter(pca_labels,question_name == question)$pc1_pos), ' (+)'))
}
plot <- ggMarginal(plot, type='density', groupColour=TRUE)
size <- 1248
if (d == 3) { # school index
size <- size * 1.5
}
png(paste0(cur_dir, '/pc', p, '_pc', p+1, '_', demodf[d] %>% names, '.png'), width = size, height = size, units = "px")
print(plot)
dev.off()
}
}
}
ANOVA, EMMeans
sig_contrasts.demos <- c()
sig_aov_posthoc <- c()
dir.create('./results', showWarnings = F)
for(question in question_names[-c(4:6,9:14)]) {
# directories
dir.create(paste('./results', question, sep = '/'), showWarnings = F)
# data subsetting
data_subset <- mergeQuestion(df, question)
data_subset_rows <- data_subset[[1]]
data_subset <- data_subset[-1]
data_subset <- data_subset[-1]
percent <- 0.4
data_subset <- data_subset[ lapply( data_subset, function(x) sum(x==0) / length(x) ) < percent ]
# generate principal components
pd <- prcomp(data_subset, retx= TRUE, center=TRUE, scale=TRUE)
# merge with demographic cols
pcs <- data.frame(rownum = data_subset_rows, pd$x[,seq(5)] %>% as.data.frame)
rownums <- pcs$rownum %>% sort
pcs <- pcs %>% arrange(rownum)
pcs <- pcs[-1]
demodfpcs <- bind_cols(df[rownums,c('key',real_demos)], pcs)
# weighting by # of responses per respondent
# ( each respondent has up to 4 responses per question, and each response is given a weight between 1/4 and 1 )
resps_by_key <- demodfpcs$key %>% table %>% as.data.frame
names(resps_by_key) <- c('key','n')
resps_by_key %<>% filter(n!=0)
for(r in seq(nrow(demodfpcs))) {
nresps <- resps_by_key %>% filter(key == demodfpcs[r,]$key) %>% .$n
demodfpcs[r, names(select(demodfpcs,matches('PC[12345]')))] <- (1/nresps) *
(demodfpcs[r, names(select(demodfpcs,matches('PC[12345]')))])
}
demodfpcs <- demodfpcs %>% select(-'key')
demodf <- demodfpcs %>% select(-matches('PC[12345]'))
pcs <- demodfpcs %>% select(matches('PC[12345]'))
# analysis
for(p in seq(pcs)) {
# creating directories
dir.create(paste('./results', question, names(pcs)[p], sep = '/'), showWarnings = F)
dir.create(paste('./results', question, names(pcs)[p], 'emmeans', sep = '/'), showWarnings = F)
dir.create(paste('./results', question, names(pcs)[p], 'anova_tukey', sep = '/'), showWarnings = F)
for(d in seq(demodf)) {
temp <- bind_cols(pcs[p], demodf[d])
names(temp) <- c('pc','demo')
if((temp$demo %>% table %>% as.data.frame %>% filter(Freq != 0) %>% nrow) == 1) next
# emmeans contrasts
lm.out <- lm(pc ~ demo, data=temp, na.action = na.exclude)
em.out <- emmeans(lm.out, ~ demo, data=temp, weights = 'proportional')
contrast.out <- contrast(em.out, method="del.eff") %>% tidy
contrast.out$p.value %<>% p.adjust(method='holm')
contrast.out %<>% mutate(question_name=question,
PC_num = paste0('PC',p),
PC_pos = pca_labels %>%
filter(question_name==question) %>%
select(matches(paste0(p,'_pos'))) %>%
pull %>% as.character,
PC_neg = pca_labels %>%
filter(question_name==question) %>%
select(matches(paste0(p,'_neg'))) %>%
pull %>% as.character,
demographic = names(demodf)[d]
)
# anova tukey
aov.out <- aov(pc ~ demo, data=temp) %>% TukeyHSD %>% tidy
aov.out %<>% select(-term) %>% mutate(question_name=question,
PC_num = paste0('PC',p),
PC_pos = pca_labels %>%
filter(question_name==question) %>%
select(matches(paste0(p,'_pos'))) %>%
pull %>% as.character,
PC_neg = pca_labels %>%
filter(question_name==question) %>%
select(matches(paste0(p,'_neg'))) %>%
pull %>% as.character,
demographic = names(demodf)[d]
)
# bind rows
options(warn=-1)
sig_contrasts.demos %<>% bind_rows(contrast.out)
sig_aov_posthoc %<>% bind_rows(aov.out)
options(warn=0)
# write csvs
write.csv(contrast.out, paste('./results',question,names(pcs)[p],'emmeans',paste0(names(demodf)[d],'.csv'),sep='/'))
write.csv(aov.out, paste('./results',question,names(pcs)[p],'anova_tukey',paste0(names(demodf)[d],'.csv'),sep='/'))
}
}
}
# # write rdas
save(sig_contrasts.demos, file='./results/emmeans_contrasts.rda')
save(sig_aov_posthoc, file='./results/anova_tukey.rda')
MANOVA
# skip aftercollege bc it causes errors (index 7)
for(question in (question_names[-c(4:6,9:14)])[-7]) {
# directories
dir.create(paste('./results', question, sep = '/'), showWarnings = F)
# data subsetting, PC
data_subset <- mergeQuestion(df, question)
data_subset_rows <- data_subset[[1]]
data_subset <- data_subset[-1]
data_subset <- data_subset[-1]
percent <- 0.4
data_subset <- data_subset[ lapply( data_subset, function(x) sum(x==0) / length(x) ) < percent ]
pd <- prcomp(data_subset, retx= TRUE, center=TRUE, scale=TRUE)
# merge with demographic cols
pcs <- data.frame(rownum = data_subset_rows, pd$x[,seq(5)] %>% as.data.frame)
rownums <- pcs$rownum %>% sort
pcs <- pcs %>% arrange(rownum)
pcs <- pcs[-1]
demodfpcs <- bind_cols(df[rownums,c('key',real_demos)], pcs)
# weighting by # of responses per respondent
resps_by_key <- demodfpcs$key %>% table %>% as.data.frame
names(resps_by_key) <- c('key','n')
resps_by_key %<>% filter(n!=0)
for(r in seq(nrow(demodfpcs))) {
nresps <- resps_by_key %>% filter(key == demodfpcs[r,]$key) %>% .$n
demodfpcs[r, names(select(demodfpcs,matches('PC[12345]')))] <- (1/nresps) *
(demodfpcs[r, names(select(demodfpcs,matches('PC[12345]')))])
}
demodfpcs <- demodfpcs %>% select(-'key')
demodf <- demodfpcs %>% select(-matches('PC[12345]'))
pcs <- demodfpcs %>% select(matches('PC[12345]'))
# manova
demo_level_nums <- demodf %>% na.omit %>% sapply(function(x)(x %>% table %>% as.data.frame %>% filter(Freq != 0) %>% nrow)) %>% as.data.frame
demo_level_nums$demo <- rownames(demo_level_nums)
names(demo_level_nums) <- c('n','demo')
filter_list <- demo_level_nums %>% filter(n < 2) %>% pull(demo)
man_demodf <- demodf %>% select(-filter_list)
print(nrow(man_demodf))
man.out <- manova(pcs %>% as.matrix ~ .,
data=man_demodf)
# write csv
write.csv(tidy(man.out), paste('./results',question,'manova_table.csv',sep='/'))
sink(paste('./results',question,'manova_output.txt',sep='/'))
print(man.out)
sink()
}