R Notebook

Loading data and setting up.

install.packages('tidyverse') # only need to do this first time - can delete or comment using "#" after
install.packages('abind') # also need to install abind
install.packages('MASS') # same
install.packages('readr') # last one - this one reads CSV files faster and has some useful defaults, though it is not necessary

Will only use one function from abind, readr, and MASS, so don’t have to load them using library

library(tidyverse) 

setwd('~/dropbox/1_research/generality') # change this to folder that has the data

data <- readr::read_csv("all_observations.csv") 
data # print the data, equivalent to using print(data)

data <- rbind(data, c(NA, '3.1', '3.1', NA)) # this is because `7C_8-3` and `7C_12-2` have no 3.1 obs - will remove this later

data[, 2:4] <- sapply(data[, 2:4], function(x) floor(as.numeric(x))) # this recodes our 10 variables to just six (0, 1, 2, 3, 4, and 5)

Look at likelihood of profile (using all data & six codes).

to_prep <- tidyr::gather(data, key, val, -id)

to_plot <- to_prep %>% 
    group_by(key, val) %>% 
    summarize(n = n()) %>% 
    filter(!is.na(val)) %>% 
    tidyr::spread(val, n, fill = 0) %>% 
    tidyr::gather(Code, the_val, -key) %>% 
    group_by(key) %>% 
    mutate(the_val = as.numeric(the_val),
           the_val = the_val / sum(the_val))

to_plot$key <- factor(to_plot$key,
                      levels = c("7C_8-3", "7C_12-2", "8B_6-2"))

to_plot$the_key <- factor(to_plot$Code)

to_plot <- to_plot %>% arrange(key, Code)

t1_freq <- as.vector(table(data$`7C_8-3`))
t2_freq <- as.vector(table(data$`7C_12-2`))
t3_freq <- as.vector(table(data$`8B_6-2`))
the_mat <- cbind(t1_freq, t2_freq, t3_freq)
row.names(the_mat) <- sort(unique(data$`8B_6-2`))
colnames(the_mat) <- names(data)[2:4] 

This throws a warning (not an error) because some cells have few (< 5) cases; it’s probably okay as long as we interpret those cells with caution.

chi_sq_output <- chisq.test(the_mat) 

Here are the z-scores (> 1.96 means cell is more likely than expected; < 1.96 means cell is less likely than expected). This finds cells (i.e., shifts) which are more (+) or less (-) expected, or as frequent as expected (=).

round(chi_sq_output$stdres, 3)  

Here is the plot.

to_plot$sig <- ifelse(as.vector(chi_sq_output$stdres) >= 1.96, "+",
                      ifelse(as.vector(chi_sq_output$stdres) <= -1.96, "-",
                             "="))

ggplot(to_plot, aes(x = key, y = the_val, color = Code, group = Code, label = sig)) +
    geom_point(size = 1.5) +
    geom_line() +
    ggrepel::geom_label_repel(show.legend = F) +
    ylab("Proportion of Responses") +
    xlab(NULL) +
    theme(text = element_text(size = 14)) +
    theme(legend.title = element_blank()) +
    labs(
        title = "Number of Responses by Code for All Observations",
        caption = "Note. +, -, and = labels indicate code is more likely than expected evaluated using a chi-square test of proportions."
    ) +
    theme(text = element_text(family = "Times")) +
    theme(plot.caption = element_text(size = 8))

# save figure in working directory
ggsave("bar_chart_all_observations.png", width = 8, height = 6)