Competencies.Rmd

---
title: "Competencies.Rmd" 
---

## Bioinf survey data analysis v1
## 2024-10-01

```{r}
library(dplyr)
library(ggplot2)
library(tidyr) 
library(plotly) 
library(data.table) # v1.16.0
library(ggrepel) # v0.9.6

citation("dplyr")
citation("ggplot2")
citation("tidyr") 
citation("plotly") 
citation("data.table") # v1.16.0
citation("ggrepel")

my_data <- read.csv("Results survey 1 October 2024 - Form responses 1.csv",
                    header=T)

row.names(my_data) <- my_data$IndividualID  # Replace 'IndividualID'  
my_data$IndividualID <- NULL  # Remove IndividualID column

# Use gsub() to remove "Level" and "(most advanced)" from cells
my_data_cleaned <- as.data.frame(lapply(my_data, function(x) {
                      gsub("Level |\\(most advanced\\)", "", x)}))
my_data_cleaned$Timestamp <- 1:nrow(my_data_cleaned) 

# Convert data to long form using pivot_longer() (modern replacement for gather)
my_data_long <- my_data_cleaned %>%
  pivot_longer( cols = -Timestamp,  # All columns except Timestamp
    names_to = "Original_Column",  # Store original column names for 
    values_to = "Competency_Score" )  # Values are the competency scores

# Assign levels based on the original column names
my_data_long <- my_data_long %>% mutate(Level = case_when(
      grepl("Profile..Core.facility.scientist.III..", Original_Column) ~ 3,
      grepl("Profile..Core.facility.scientist.II..", Original_Column) ~ 2,
      grepl("Profile..Core.facility.scientist.I..", Original_Column) ~ 1,
      grepl("Managerial.role", Original_Column) ~ 4 ))

# Extract competency labels (like A3, B3, C3) from the Original_Column names
my_data_long <- my_data_long %>% mutate(Competency = case_when(
    grepl("\\.A3\\.", Original_Column) ~ "A3",
    grepl("\\.B3\\.", Original_Column) ~ "B3",
    grepl("\\.C3\\.", Original_Column) ~ "C3",
    grepl("\\.D3\\.", Original_Column) ~ "D3",
    grepl("\\.E3\\.", Original_Column) ~ "E3",
    grepl("\\.F3\\.", Original_Column) ~ "F3",
    grepl("\\.G3\\.", Original_Column) ~ "G3",
    grepl("\\.H3\\.", Original_Column) ~ "H3",
    grepl("\\.I3\\.", Original_Column) ~ "I3",
    grepl("\\.J3\\.", Original_Column) ~ "J3",
    grepl("\\.K3\\.", Original_Column) ~ "K3",
    grepl("\\.L3\\.", Original_Column) ~ "L3",
    grepl("\\.M3\\.", Original_Column) ~ "M3",
    grepl("New..Identify.and.support.users..needs.", Original_Column) ~ "N1",
    grepl("New..Project.management.", Original_Column) ~ "O1",
    grepl("New..People.management..focusing.on.staff..",Original_Column)~"P1",
    grepl("New..Collaborator.engagement.", Original_Column) ~ "Q1",
    grepl("New..Training.", Original_Column) ~ "R1",
    grepl("New..Leadership.", Original_Column) ~ "S1" ))

my_data_long2 <- my_data_long
my_data_long2$Original_Column <- NULL
my_data_long2 <- my_data_long2[-1,]
# my_data_long2[my_data_long2 == "Not Applicable"] <- NA
my_data_long2[my_data_long2 == "Not Applicable"] <- "0" # implicit coercion!!
my_data_long2$Competency_Score <- as.numeric(my_data_long2$Competency_Score)
my_data_long2$Level <- as.numeric(my_data_long2$Level)
View(my_data_long2)
write.csv(my_data_long2, "cleaned_long_data.v2.csv", row.names = FALSE)

## heatmap

mode_stat <- function(x) {  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))] }

# Use the custom mode function in your summarise() call
my_data_summary <- my_data_long2 %>%
  group_by(Competency, Level) %>%
  summarise(mean_score = mode_stat(Competency_Score)) %>%
  ungroup() %>%
  mutate(Competency = factor(Competency, levels = sort(unique(Competency)))) 
my_data_summary <- my_data_summary[-77,]
my_data_summary <- my_data_summary %>% mutate(mean_score = factor(mean_score))

pdf("heatmap.mode.v1.pdf", width = 3, height = 6)
ggplot(my_data_summary, aes(x=factor(Level), y=Competency, fill=mean_score)) +
  geom_tile(color = "white") +  
  geom_text(aes(label = mean_score), color = "white", size = 3) +  
  labs(x = "Level", y = "Competency") + theme_minimal() +
  scale_y_discrete(limits = rev(levels(my_data_summary$Competency))) +   
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.text.y = element_text(size = 10)) + coord_fixed() +
  scale_fill_manual(values = c("0"="black", "1" = "grey20", "2" = "grey40",
                  "3" = "grey60"), name = "Score") 
dev.off()

# black to white
pdf("heatmap.mode.v2.pdf", width = 3, height = 6)
ggplot(my_data_summary, aes(x=factor(Level), y=Competency, fill=mean_score)) +
  geom_tile(color = "white") +  
  geom_text(aes(label = mean_score), color = "white", size = 3) +  
  labs(x = "Level", y = "Competency") + theme_minimal() +
  scale_y_discrete(limits = rev(levels(my_data_summary$Competency))) +   
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.text.y = element_text(size = 10)) + coord_fixed() +
  scale_fill_manual(values = c("0"="grey60", "1" = "grey40", "2" = "grey20",
                  "3" = "black"), name = "Score")  
dev.off()

my_data_summary <- my_data_long2 %>% group_by(Competency, Level) %>%
  summarise(mean_score = median(Competency_Score, na.rm =T)) %>% ungroup() %>%
  mutate(Competency = factor(Competency, levels = sort(unique(Competency))))
my_data_summary <- my_data_summary[-77,] 

pdf("heatmap.median.v1.pdf", width = 3, height = 6)
ggplot(my_data_summary, aes(x=factor(Level), y=Competency, fill=mean_score)) +
  geom_tile(color = "white") +  
  geom_text(aes(label = mean_score), color = "white", size = 3) +  
  labs(x = "Level", y = "Competency") + theme_minimal() +
  scale_y_discrete(limits = rev(levels(my_data_summary$Competency))) +   
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.text.y = element_text(size = 10)) + coord_fixed()  +
  scale_fill_gradient(low = "grey70", high = "black", limits = c(0, 3),
  name="Score")
dev.off()

# black to white
pdf("heatmap.median.v2.pdf", width = 3, height = 6)
ggplot(my_data_summary, aes(x=Level, y=Competency, fill=mean_score)) +
  geom_tile(color = "white") +  
  geom_text(aes(label = mean_score), color = "white", size = 3) +  
  labs(x = "Level", y = "Competency") + theme_minimal() +
  scale_y_discrete(limits = rev(levels(my_data_summary$Competency))) +   
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.text.y = element_text(size = 10)) + coord_fixed()  +
  scale_fill_gradient(low = "black", high = "grey70", limits = c(0, 3),
  name="Score")
dev.off()

# Get summary statistics  
summary_by_level <- my_data_long2 %>% group_by(Level, Competency) %>%  
  summarise(Mean_Score = mean(as.numeric(Competency_Score), na.rm=T),
            Median_Score =median(as.numeric(Competency_Score), na.rm=T),
            SD_Score = sd(as.numeric(Competency_Score), na.rm = T))
# print(summary_by_level)

# change to Spearman's
# keep colour scheme consistent

set1 <- na.omit(my_data_long2) %>% group_by(Competency) %>%
  summarise( Correlation = cor.test(Level, Competency_Score,
                                    method = c(  "spearman"))$estimate, 
    P_Value = cor.test(Level, Competency_Score)$p.value, .groups = 'drop')
set1$adj_P_Value <- p.adjust(set1$P_Value, "BH", length(set1$P_Value))
set1$Correlation <- round(set1$Correlation,2) 
print(set1)
write.csv(set1,  "correlations.v2.csv")

# Calculate means for each competency and level
mean_scores <- my_data_long2 %>% group_by(Competency, Level) %>%
  summarise(Mean_Score=mean(Competency_Score, na.rm=T), .groups='drop')
mean_scores <- mean_scores[-77,] # remove NA
print (mean_scores)

# Calculate global means for each  level as baseline for comparison
mean_level <- my_data_long2 %>% group_by(Level) %>%
  summarise(Mean_Score=mean(Competency_Score, na.rm=T), .groups='drop')
mean_level <- mean_level[-5,] # remove NA
mean_level2 <- rbind(mean_level, mean_level,mean_level, mean_level,
                     mean_level, mean_level,mean_level, mean_level,
                     mean_level, mean_level,mean_level, mean_level,
                     mean_level, mean_level,mean_level, mean_level,
                     mean_level, mean_level,mean_level) # 19 copies

# Create the plot with the global mean line with CIs
mean_scores_ci <- my_data_long2 %>%
  group_by(Competency, Level) %>%
  summarise(
    Mean_Score = mean(Competency_Score, na.rm = TRUE),
    Lower_CI=Mean_Score - qt(0.975, df=n()-1)*sd(Competency_Score,na.rm=T)/sqrt(n()),  
    Upper_CI=Mean_Score + qt(0.975, df=n()-1)*sd(Competency_Score,na.rm=T)/sqrt(n()),  
    .groups = 'drop' )
mean_scores_ci <- mean_scores_ci[-77,] # remove NA

set1 <- set1 %>% rename(Correlation = Correlation) # why this? 
set1$Correlation <- round(set1$Correlation,2)
mean_scores_ci <- mean_scores_ci %>%
  left_join(set1 %>% select(Competency, Correlation), by = "Competency")

mean_scores_ci$Correlation <- round(mean_scores_ci$Correlation, 2)
# for plotting
competency_colors <- c(
  "A3" = "#F28E82",  # Light pink/red
  "B3" = "#F4A261",  # Orange
  "C3" = "#E9C46A",  # Yellow
  "D3" = "#A8D08D",  # Green
  "E3" = "#A5A58D",  # Olive
  "F3" = "#8EC07C",  # Light green
  "G3" = "#89C4C8",  # Light cyan
  "H3" = "#52B788",  # Deep green
  "I3" = "#00AFFF",  # Light blue
  "J3" = "#8B9DC3",  # Light purple
  "K3" = "#C39BD3",  # Lavender
  "L3" = "#B587FF",  # Light purple
  "M3" = "pink",  # Light pink/purple
  "N1" = "black",  # Light blue/purple
  "O1" = "grey70",  # Light purple
  "P1" = "grey40",  # Pink
  "Q1" = "yellow",  # Pink/purple
  "R1" = "yellowgreen",  # Rose pink
  "S1" = "#E9C35B")   # Deep pink

pdf("Fig_2A_Competency_change_with_CI.pdf", width = 6, height = 6)
ggplot(mean_scores_ci, aes(x = Level, y = Mean_Score,
                           group = Competency, color = Competency)) + 
  geom_ribbon(aes(ymin = Lower_CI, ymax = Upper_CI), fill="grey", alpha=0.2) +
  geom_line(linewidth = 1.3, alpha = 0.7) +   
  geom_point(size = 1.2, alpha = 0.7) +   
  geom_line(data = mean_level2, aes(x = Level, y = Mean_Score), color="black",
            linewidth = 0.8, linetype = "dashed", inherit.aes = F) + 
  facet_wrap(~ Competency) + 
  geom_text(aes(x = 2, y = 0.3, label = paste0("r=", round(Correlation, 2))), 
            color = "black", size = 3, hjust = 0, vjust = 0) +   
  labs(x = "Stage", y = "Mean Score") + ylim(0, 3) + 
  theme_minimal() +  
  theme(legend.position = "bottom", 
        strip.text = element_text(size = 10, face = "bold")) +
  scale_color_manual(values =competency_colors)
dev.off()

mean_scores_ci <- mean_scores_ci %>%
  mutate(Competency_Group = case_when(
    Competency %in% c("A3", "B3", "C3") ~ "A3-C3: Bioscience",
    Competency %in% c("D3", "E3", "F3") ~ "D3-F3: Data science",
    Competency %in% c("G3", "H3", "I3") ~ "G3-I3: Computer science",
    Competency %in% c("J3", "K3", "L3", "M3") ~ "J3-M3: Professional conduct",
    Competency %in% c("N1", "O1", "P1", "Q1", "R1", "S1") ~
      "N1-S1: Core Facility-focused" )) # Strategy/Management

m2 <- mean_scores_ci %>% filter(Level == 2)
m3 <- mean_scores_ci %>% filter(Level == 3)
m4 <- mean_scores_ci %>% filter(Level == 4)
pdf("Fig_2C_Competency_means.pdf", width =10.2, height =5)
mean_scores_ci %>% filter(Level == 1) %>%
  ggplot(aes(x = Competency, y = Mean_Score, group = Competency)) + 
  geom_point(size = 5, alpha=0.5, aes(color = "I", shape = "I", fill = "I"),
             show.legend = TRUE) +
  geom_point(data = m2, aes(x = Competency, y = Mean_Score, color = "II",
         shape = "II", fill = "II"), size = 5, alpha=0.5, show.legend = T) +
  geom_point(data = m3, aes(x = Competency, y = Mean_Score, color = "III",
           shape = "III", fill = "III"), size = 5, alpha=0.5, show.legend =T)+
  geom_point(data = m4, alpha=0.7, aes(x=Competency, y=Mean_Score,
          color="M managerial", shape = "M managerial", fill ="M managerial"), 
             size = 5, stroke = 1.5, show.legend = TRUE) +  
  geom_text(data = m4, aes(x = Competency, y = Mean_Score, label = "M"), 
            color = "black", size =4, fontface = "bold", show.legend =F) +
  labs(x = "Competency", y = "Mean Score") +  ylim(0, 3) + theme_minimal() +  
  theme( legend.position = "bottom", legend.box = "vertical", 
    panel.spacing.x = unit(1, "lines"), legend.subtitle=element_text(size=4))+
  scale_color_manual( name = "Career",  guide = guide_legend(order = 1), 
    values = c("I" = "#8ecae6", "II" = "#219ebc", "III" = "#023047",
               "M managerial" = "#ee9b00") )  +
  scale_fill_manual( name = "Career",    guide = guide_legend(order = 1),
    values = c("I" = "#8ecae6", "II" = "#219ebc", "III" = "#023047",
               "M managerial" = "#ee9b00") ) +
  scale_shape_manual( name = "Career",  guide = guide_legend(order = 1),
    values = c("I" = 21, "II" = 24, "III" = 22, "M managerial" = 21)) +
  facet_grid(    ~Competency_Group, scales = "free_x",  space = "free_x" )
dev.off()

##### get outliers from above

# Merge the confidence intervals with the global mean for comparison
comparison_data <- mean_scores_ci %>%
  left_join(mean_level2, by="Level", suffix=c("_ci", "_global"))

# Identify cases where both Lower_CI and Upper_CI are either above or below the global mean
comparison_data <- comparison_data %>%
  mutate( Comparison = case_when(
      Lower_CI > Mean_Score_global & Upper_CI > Mean_Score_global ~ "Above",
      Lower_CI < Mean_Score_global & Upper_CI < Mean_Score_global ~ "Below",
      TRUE ~ "Within"  ) )

outliers <- comparison_data %>% filter(Comparison != "Within")

# View the output: possible 19x4 (76 instances in total)

# Perceived as less important at the given level
subset(unique(outliers), Comparison=="Below") # 11 instances
# O1, P1, S1 - initially

# Percevied as more important at the given level
subset(unique(outliers), Comparison=="Above") # 16 instances
# D3, E3, F3, H3, L3 - initially

## correlations

# Least correlated => r<0.24
# needed less at managerial stage:
#     B3 = Prepare life science data for computational analysis
#     H3 = Make appropriate and efficient use of scripting and programming languages

# more correlated =>  r > 0.5
# career stage progress 
#    S3 - correlated with later stages - "New..Leadership."
#    O3 - somewhat correlated with later stages - "Project.management."
#    P3 - correlated with later stages - "People.management..focusing.on.staff.
#    Q3 - "New..Collaborator.engagement."
#    K3 - Communicate meaningfully with a range of audiences - within ...

#### check ranks

mean_scores_ci %>% filter(Level==1) %>% arrange(desc(Mean_Score)) %>%
  select(Competency, Mean_Score)
mean_scores_ci %>% filter(Level==2) %>% arrange(desc(Mean_Score)) %>%
  select(Competency, Mean_Score)

### parallel plot

mean_scores_ci <- mean_scores_ci %>%
  group_by(Level) %>%
  mutate(Rank = rank(-Mean_Score))  # Rank descending by Mean_Score

# Ensure the data is ordered correctly
mean_scores_ci <- mean_scores_ci[order(mean_scores_ci$Level,
                                       mean_scores_ci$Rank), ]

### pairwise correlations 

my_data_wide <- my_data_long2 %>% 
  spread(key=Competency, value=Competency_Score) %>%
  filter(!is.na(Level)) 
my_data_wide2 <- my_data_wide[,-22]

correlation_test <- function(data) {
  n <- ncol(data)
  cor_matrix <- matrix(NA, n, n)
  p_matrix <- matrix(NA, n, n)
  for (i in 1:(n-1)) {
    for (j in (i+1):n) {
      res <- cor.test(unlist(data[,i]),unlist(data[,j]), method="spearman",
                      use = "pairwise.complete.obs")
      cor_matrix[i, j] <- res$estimate
      p_matrix[i, j] <- res$p.value  } }
  list(correlation = cor_matrix, p_value = p_matrix) }

result2 <- correlation_test(my_data_wide2[,c(-1,-2)])
result <- result2$correlation
print (result)
p_value_matrix_adj <- matrix(p.adjust(as.vector(result2$p_value), method="BH"), 
                      nrow=ncol(result2$p_value), ncol = ncol(result2$p_value))

# Fill lower triangles to make symmetric matrices 
#cor_melt <- melt(result, na.rm = T)   # old
cor_melt <- reshape2::melt(result, na.rm =T)

competency_names <- colnames(my_data_wide2)[3:21] 
cor_melt$Var1 <- competency_names[cor_melt$Var1]
cor_melt$Var2 <- competency_names[cor_melt$Var2] 
str(cor_melt) # Var1 Var2 value 

pdf("Figure_3_Pairwise_correlation.pdf", width=6, height=6)
cor_melt %>%  ggplot(aes(Var1, Var2, fill = value)) +   geom_tile() +
  geom_text(aes(label = round(value, 2)), color = "black", size=2.1)+
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", 
                    midpoint = 0.5, limit = c(0, 1), space = "Lab",
                       name = "Correlation") +
  theme_minimal() +  labs(x = "Competency", y = "Competency") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1),
        axis.text.y = element_text())
dev.off()

subset(cor_melt, value>0.75) # get top ones
subset(cor_melt, value<.3)  # get lower ones

## PCA

competency_data_wide2 <- t(my_data_wide2)[c(-1,-2),]# get constant cols
constant_cols <- apply(competency_data_wide2, 2, function(x) var(x) == 0) 
competency_data_wide2_clean <- competency_data_wide2[, !constant_cols]  
pca_result <- prcomp(competency_data_wide2_clean, scale = T) 

pdf("Figure_4_PCA.pdf", width=7, height=4)
plot(pca_result$x[, 1], pca_result$x[, 2], cex=1.4, pch = 19, xlim=c(-10,22),
     xlab = paste0("PC1 ", round(100*summary(pca_result)$importance[2,1],0), "%"),
     ylab = paste0("PC2 ", round(100*summary(pca_result)$importance[2,2],0), "%"), 
     col = competency_colors[rownames(pca_result$x)] )
grid(lty="dotted", lwd=2) 
text(pca_result$x[, 1], pca_result$x[, 2], labels = rownames(pca_result$x), 
     pos = c(3,1,1,1,4,1,3,1,3,1,2,1,1,4,1,1,1,2,3), cex = 1.3,
     col = competency_colors[rownames(pca_result$x)]) 
legend("topright", legend = rownames(pca_result$x), ncol=2,
       col = competency_colors[rownames(pca_result$x)], pch = 19, cex=1)
dev.off()

pdf("Figure_PCA_PC3_PC4.pdf", width=7, height=4)
plot(pca_result$x[, 3], pca_result$x[, 4], cex=1.4, pch = 19, xlim=c(-10,22),
     xlab = paste0("PC3 ", round(100*summary(pca_result)$importance[2,3],1), "%"),
     ylab = paste0("PC4 ", round(100*summary(pca_result)$importance[2,4],1), "%"), 
     col = competency_colors[rownames(pca_result$x)] )
grid(lty="dotted", lwd=2) 
text(pca_result$x[, 3], pca_result$x[, 4], labels = rownames(pca_result$x), 
     pos = c(3,1,1,1,4,1,3,1,3,1,2,1,1,4,1,1,1,2,3), cex = 1.3,
     col = competency_colors[rownames(pca_result$x)]) 
legend("topright", legend = rownames(pca_result$x), ncol=2,
       col = competency_colors[rownames(pca_result$x)], pch = 19, cex=1)
dev.off()

###
```
Profile: Core facility scientist I [A3: Work at depth in at least one technical area aligned with the life sciences]	Profile: Core facility scientist I [B3: Prepare life science data for computational analysis]	Profile: Core facility scientist I [C3: Have a positive impact on scientific discovery through bioinformatics]	Profile: Core facility scientist I [D3: Use data science methods suitable for the size and complexity of the data]	Profile: Core facility scientist I [E3: Manage own and others’ data according to community standards and principles]	Profile: Core facility scientist I [F3: Make appropriate use of bioinformatics tools and resources]	Profile: Core facility scientist I [G3: Contribute effectively to the design and development of user-centric bioinformatics tools and resources]	Profile: Core facility scientist I [H3: Make appropriate and efficient use of scripting and programming languages]	Profile: Core facility scientist I [I3: Construct, manage and maintain bioinformatics computing infrastructure of varying complexity]	Profile: Core facility scientist I [J3: Comply with professional, ethical, legal and social standards and codes of conduct relevant to computational biology]	Profile: Core facility scientist I [K3: Communicate meaningfully with a range of audiences - within and beyond your profession]	Profile: Core facility scientist I [L3: Work effectively in teams to accomplish a common goal]	Profile: Core facility scientist I [M3: Engage in continuing professional development in bioinformatics]