---
tags: Rscripts
---

# Heat maps of t-tests X sex
### Author: Jocelyn P. Colella

**Hypothesis:** Males and Females will have significantly different physiological responses (e.g., dependent variables) during similar experiments

**Plots**
Heatmaps of bonferonni adjusted p-values for t-tests between males and females across each experiment (in total) and across each dependent variable:

- _males and females are MOSTLY significantly different_
- _M and F are NOT different in VO2 or EE under HOT experimental conditions_

Heatmaps of Bonferonni adjusted p-values for t-tests between males and females across each experiment (BL, hot, cold), day and night, and across each dependent variable:
- _Males and females significantly different for MOST comparisons_
- _Hot conditions are freq not different btwn M and F (but not universal):
H2O hot_night_M/hot_night_F
VCO2 hot_day_M/hot_day_F_
- All other NON-signif comparisons are between different treatment groups (e.g., hot_night_M vs. cold_day_M)

# R Code, by JPC
dependent_variables = c("EE", "RQ", "VO2", "VCO2", "H2Omg")

### M and F (all)
MF_all_list = c("all_noOL_F", "all_noOL_M",
                "BL_noOL_F", "BL_noOL_M",
                "hot_noOL_F", "hot_noOL_M",
                "cold_noOL_F", "cold_noOL_M")
               
#Write header to mean/sd file
meanSD_header <- paste('DV', 'Dataset', 'mean', 'sd', sep=',')
write.table(meanSD_header, "mean_sd_eachTreatment_Xsex.csv", sep=',', col.names = FALSE, row.names = FALSE, quote = FALSE)

#write head to t test results file
t_header <- paste('Data1', 'Data2', 'DV', 'p-value', 'unadj. p', sep=',')
write.table(t_header, "ttest_results_Xsex.csv", sep=',', col.names = FALSE, row.names = FALSE, quote = FALSE)


#For each experiment
#Calcualte the average and SD of each dependent variable
count = 1
for(df in MF_all_list){
  print(df)
  for (DV in dependent_variables){
    print(DV)
    this_mean = mean(get(df)[[DV]])
    this_sd = sd(get(df)[[DV]])
    line = paste(DV, df, this_mean, this_sd, sep=',')
    print(line)
    print('\n')
    write.table(line, "mean_sd_eachTreatment_Xsex.csv", sep = ',', append=TRUE, col.names = FALSE, row.names = FALSE, quote = FALSE)
    
    #two sample t test - test whether the mean differs from other exp
    templist <- list("all_noOL_F", "all_noOL_M",
                     "BL_noOL_F", "BL_noOL_M",
                     "hot_noOL_F", "hot_noOL_M",
                     "cold_noOL_F", "cold_noOL_M")
    
    # Adjust for multiple comparisons (Bonferonni)
    k <- length(templist)
    num_pw_comparisons <- k*(k-1) / 2 #where k is the number of conditions being compared
    for(data in templist){
      ts <- t.test(get(df)[[DV]], get(data)[[DV]], conf.level = 0.95)
      new_p <- p.adjust(ts$p.value, method = "bonferroni", num_pw_comparisons)
      this_line=paste(df, data, DV, new_p, ts$p.value, sep = ',')
      print(this_line)
      write.table(this_line, "ttest_results_Xsex.csv", sep=',', append = TRUE, col.names = FALSE, row.names = FALSE, quote = FALSE)
    }
    count = count + 1
  }
}

##REPEAT for each experiment split into night and day

To make paneled heat maps:

#Read in ttest results produced above
ttest_xSex <- read.csv("ttest_results_Xsex.csv")

#replace _noOL in the variable names for plotting
ttest_xSex_clean <- ttest_xSex %>% 
      mutate(Data1 = str_replace(Data1, "_noOL", "")) %>%
      mutate(Data2 = str_replace(Data2, "_noOL", ""))

# GENERATE PLOT FOR EACH DEPENDENT VARIABLE
ee_xSex = ttest_xSex_clean[ttest_xSex_clean$DV == 'EE', ]  
ee_mat <- dcast(ee_xSex, Data1~Data2, value.var="p.value")

#Remove female (F) columns and male (M) rows)
col2drop <- c("all_F", "BL_F", "cold_F", "hot_F")
eemat2 <- ee_mat[ , !(names(ee_mat) %in% col2drop)]
eemat3 <- eemat2[-c(1,3,5,7), ]
eemat3_m<- melt(eemat3)

MF_ee <- ggplot(data = eemat3_m, aes(Data1, variable, fill = value))+
  geom_tile(color = "white")+
  scale_fill_viridis(option="viridis", alpha = 0.75, limit = c(0, 0.8), na.value="white") +
  labs(y="Males", x = "Females") +
  theme_minimal()+ 
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        legend.justification = c(0.9, 0),
        legend.position = "none",
        legend.direction = "horizontal",
        legend.title = element_blank())+#,
  coord_fixed() +
  geom_text(aes(x=Data1,y=variable,label=round(value,2)),size=3,color="black") #text IN heatmap boxes
MF_ee