gates_npi_temp_rr_v1.Rmd

---
title: "gates_data_rr_temp_npi"
author: "aj_mitchell"
date: "June 15, 2019"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r, echo=F, message=FALSE, warning=FALSE}
library(tidyselect)
library(ggplot2)
library(broom)
library(dplyr)
library(corrplot)
library(kableExtra)
library(knitr)
library(Hmisc)
library(infer)
library(plotly) 
library(skimr)
library(stringr)
library(reshape2)
library(som)
library(manipulate)
library(lavaan)
library(semPlot)
library(FactoMineR)
library(factoextra)
library(readxl)
library(formattable)
library(DT)
library(xtable)
library(wesanderson)
library(lessR)
library(RColorBrewer)
```


```{r plotting}
gates_temp_NPI_combined_v1 <- read_excel("gates_temp_NPI_combined_v1.xlsx") 
#####^do not use!! but you can still use the same datafile as long as the subjects with missing data are cleaned like below
#View(gates_temp_NPI_combined_v1)
gates_temp_NPI_combined_v1$subj_id <- as.factor(gates_temp_NPI_combined_v1$subj_id)
##the v2 will have the subjects removed that are causing the NA, also stereotypy behaviors were removed becuase they did not have any variance and they were causing their own error

gates_temp_npi_combined_v2 <- gates_temp_NPI_combined_v1 %>%
  filter(!subj_id %in% c(35408 , 34404 , 34412))

selected_gates_temp_NPI_combined_v2 <- gates_temp_npi_combined_v2 %>%
  select(- subj_id, - group)
##use selected_gates_temp_NPI_combined_v2 to look at group differences based on composites 
View(gates_temp_npi_combined_v2)

#selected_gates_temp_NPI_combined_v1 <- gates_temp_NPI_combined_v1 %>%
  #select(- subj_id, - group)

#View(selected_gates_temp_NPI_combined_v2)

#str(selected_gates_temp_NPI_combined_v2)

```


```{r}

####everything below should use gates_temp_npi_combined_v2! <<<---

corr_gates_temp_npi_v2 <- rcorr(as.matrix(selected_gates_temp_NPI_combined_v2))

#View(corr_gates_temp_npi_v1$P)

mat_1_v2 <- corr_gates_temp_npi_v2$r
ordered_mat_1_v2 <- corReorder(mat_1_v2)

npi_temp_structure_allcorrs <- corrplot(mat_1_v2,
         method = "color",
         col = brewer.pal(n=10, name ="RdYlBu"),
         type = "lower",
         order = "hclust", 
         tl.col = "black",
         tl.cex = .7)

npi_temp_structure_allcorrs

View(mat_1_v2)

cor_order <- corReorder(mat_p2, heat_map = T)

melted_corr_gates_temp_npi_v2 <- melt(corr_gates_temp_npi_v2$r)


melted_pval_gates_temp_npi_v2 <- melt(corr_gates_temp_npi_v2$P) 

datatable(formattable(melted_pval_gates_temp_npi_v2))

mat_p2_v2 <- corr_gates_temp_npi_v2$P
##this is how to move forward with correlation visuals 
pvalsonly_npi_temp_structure <- corrplot(corr_gates_temp_npi_v2$r,
         p.mat = corr_gates_temp_npi_v2$P,
         insig = "blank",
         method = "color",
         col = brewer.pal(n=10, name ="RdYlBu"),
         type = "lower",
         order = "hclust", 
         tl.col = "black",
         tl.cex = .7)

#View(melted_pval_gates_temp_npi_v1)

```


```{r with normalized structural data}
#this chunk will be doing the same steps above except with normalized imaging data
#i want to see if it changes the relationships of the output
#i will show both to chris if so to see which one is appropriate
######doesnt matter if structure is normalized or not


###the figures here can be used for the meeting tomorrow
###these figures contain no composites and are representations of raw relationships. 

##use this file for raw comparisons ((7/23))

gates_temp_NPI_combined_ALLNORM_v2 <- read_excel("gates_temp_NPI_combined_ALLNORM_v2.xlsx")
#View(gates_temp_NPI_combined_ALLNORM_v2)

selected_gates_temp_NPI_combined_ALLNORM_v2 <- gates_temp_NPI_combined_ALLNORM_v2 %>%
  filter(!subj_id %in% c(35408 , 34404 , 34412)) %>%
  select(-subj_id, -group, -Sex)

#View(selected_gates_temp_NPI_combined_ALLNORM_v2)

corr_gates_ALLNORM_v2 <- rcorr(as.matrix(selected_gates_temp_NPI_combined_ALLNORM_v2))

#View(corr_gates_ALLNORM_v2$r)

pval_gates_ALLNORM_v2 <- melt(corr_gates_ALLNORM_v2$P) %>%
  filter(value <= .05)

#View(pval_gates_ALLNORM_v2)

corrplot(corr_gates_ALLNORM_v2$r,
         method = "color",
         col = brewer.pal(n=10, name ="RdYlBu"),
         order = "hclust", 
         tl.col = "black",
         tl.cex = .5)

corrplot(corr_gates_ALLNORM_v2$r,
         p.mat = corr_gates_ALLNORM_v2$P,
         insig = "blank",
         method = "color",
         col = brewer.pal(n=10, name ="RdYlBu"),
         type = "lower",
         order = "hclust", 
         tl.col = "black",
         tl.srt = 25,
         tl.cex = .5)

datatable(formattable(corr_gates_ALLNORM_v2$n))
#png(filename)

#next steps are to incorporate the rr data with the structure
#standardize the structural data with temp and NPI data

```



```{r}
#only use the average just to get an idea of what it looks like with
#1) structure
#2) temperament
#3) novel peer performance


#i do the agitation data separate becuase it misses data from different subjects; it's more informative on its own
agitation_percent_RRdata_plus_structure <- read_excel("agitation_percent_RRdata_plus_structure.xlsx")
View(agitation_percent_RRdata_plus_structure)
agitation_percent_RRdata_plus_structure$Sex <- as.factor(agitation_percent_RRdata_plus_structure$Sex)
agitation_percent_RRdata_plus_structure$Diet <- as.factor(agitation_percent_RRdata_plus_structure$Diet)

selected_agitation <- agitation_percent_RRdata_plus_structure %>%
  select(-"Monkey_#", -match, - Sex, - Diet)

corr_agitation <- rcorr(as.matrix(selected_agitation))

corrplot(corr_agitation$r,
         method = "color",
         col = brewer.pal(n=10, name ="RdYlBu"),
         type = "lower",
         order = "hclust", 
         tl.col = "black",
         tl.srt = 25,
         tl.cex = .5)

corrplot(corr_agitation$r,
         p.mat =  corr_agitation$P,
         insig = "blank",
         method = "color",
         col = brewer.pal(n=10, name ="RdYlBu"),
         type = "lower",
         order = "hclust", 
         tl.col = "black",
         tl.srt = 25,
         tl.cex = .5)


#%time agitated is interesting
```

```{r}
#the rr data
####heavy grains of salt because i am using averaged classification scores. 
avg_rr_data_plus_structure <- read_excel("avg_rr_data_plus_structure.xlsx")
str(avg_rr_data_plus_structure)

avg_rr_data_plus_structure$Sex <- as.factor(avg_rr_data_plus_structure$Sex)
avg_rr_data_plus_structure$Diet <- as.factor(avg_rr_data_plus_structure$Diet)

selected_rr_and_structure <- avg_rr_data_plus_structure %>%
  select(-"Monkey_#", - Sex, - Diet, - DOB)

corr_rr_plus_structure <- rcorr(as.matrix(selected_rr_and_structure))

corrplot(corr_rr_plus_structure$r,
         method = "color",
         col = brewer.pal(n=10, name ="RdYlBu"),
         type = "lower",
         order = "hclust", 
         tl.col = "black",
         tl.srt = 25,
         tl.cex = .5)

corrplot(corr_rr_plus_structure$r,
         p.mat = corr_rr_plus_structure$P,
         insig = "blank",
         method = "color",
         col = brewer.pal(n=10, name ="RdYlBu"),
         order = "hclust", 
         tl.col = "black",
         tl.cex = .5)

melted_rr_pval <- melt(corr_rr_plus_structure$P) %>%
  filter(value <= .05)

melted_rr_corrvals <- melt(corr_rr_plus_structure$r)


datatable(formattable(melted_rr_corrvals))

write.csv(melted_rr_pval, "response_reactivity_pvals.csv")

```



```{r}
##next need to do PCAs and SEMs on new outputs using the calculated thresh holds 


```

========
new post meeting changes - 6/24
========

```{r}
#this chunk is for creating the interactive plots to put on box
Pvalue_matrix_temp_stuct_NPI_interactive <- ggplotly(ggplot(pval_gates_ALLNORM_v2, aes(x=Var2, y=Var1, fill=value)) +
         geom_tile(stat='identity') +
    scale_fill_gradientn(colours = topo.colors(5)) + 
  theme_minimal() +
  theme(axis.text.x = element_text(margin = margin(b = .5, t=.5),
                                   angle = 90,
                                   vjust = .16,
                                   size = 10,
                                   hjust = 1
                                    )) +
  theme(axis.text.y = element_text(margin = margin(b = .5, t=.5),
                                   vjust = .16,
                                   size = 10,
                                   hjust = 1
                                    )))
Pvalue_matrix_temp_stuct_NPI_interactive

View(pval_gates_ALLNORM_v2)

write.csv(pval_gates_ALLNORM_v2, "pvals_gates_ALLNORM_csv.csv")

ggplot(pval_gates_ALLNORM_v2, aes(x=Var2, y=Var1, fill=value)) +
  geom_tile(stat='identity') +
  scale_fill_gradientn(colours = topo.colors(5)) + 
  theme_minimal() +
  theme(axis.text.x = element_text(margin = margin(b = .5, t=.5),
                                   angle = 90,
                                   vjust = .16,
                                   size = 10,
                                   hjust = 1
                                    )) +
  theme(axis.text.y = element_text(margin = margin(b = .5, t=.5),
                                   vjust = .16,
                                   size = 10,
                                   hjust = 1
                                    ))

melt2gether <- melt(corr_gates_ALLNORM_v2$r, corr_gates_ALLNORM_v2$P)


library(PerformanceAnalytics)

chart.Correlation(selected_gates_temp_NPI_combined_v2, histogram = T, pch = 20)

cormat <- corr_gates_ALLNORM_v2$r
pmat <- corr_gates_ALLNORM_v2$P

flattenCorrMatrix <- function(cormat, pmat) {
  ut <- upper.tri(cormat)
  data.frame(
    row = rownames(cormat)[row(cormat)[ut]],
    column = rownames(cormat)[col(cormat)[ut]],
    cor  =(cormat)[ut],
    p = pmat[ut]
    )
}

flatmat_ALL <- flattenCorrMatrix(corr_gates_ALLNORM_v2$r, corr_gates_ALLNORM_v2$P)

View(flatmat_ALL_v2)

flatmat_ALL_v2 <- flatmat_ALL %>%
  filter(p <= .05)

write.csv(flatmat_ALL_v2, "pval_corrval_stucture_temp_npi.csv")

```


```{r}

cormat1 <- corr_rr_plus_structure$r
pmat1 <- corr_rr_plus_structure$P

flattenCorrMatrix1 <- function(cormat1, pmat1) {
  ut <- upper.tri(cormat1)
  data.frame(
    row = rownames(cormat1)[row(cormat1)[ut]],
    column = rownames(cormat1)[col(cormat1)[ut]],
    cor  =(cormat1)[ut],
    p = pmat1[ut]
    )
}

flatmat_rr_v1  <- flattenCorrMatrix1(corr_rr_plus_structure$r, corr_rr_plus_structure$P)
#v2 = pvals > .05 removed
flatmat_rr_v2 <- flatmat_rr_v1 %>%
  filter(p <= .05)

View(flatmat_rr_v2)

write.csv(flatmat_rr_v2, "pval_corrval_rr_structure.csv")
```


```{r}
###goals for 7/23
####look at group differences for year 1 in temperament
####look at group differences for cortisol tomorrow morning


skim(gates_temp_npi_combined_v2)

skim(gates_temp_NPI_combined_ALLNORM_v2)

gates_temp_NPI_combined_ALLNORM_v2 <- read_excel("gates_temp_NPI_combined_ALLNORM_v2.xlsx")
#View(gates_temp_NPI_combined_ALLNORM_v2)

subremoved_gates_temp_NPI_combined_ALLNORM_v2 <- gates_temp_NPI_combined_ALLNORM_v2 %>%
  filter(!subj_id %in% c(35408 , 34404 , 34412))

View(subremoved_gates_temp_NPI_combined_ALLNORM_v2)


pr50_gates_temp_NPI_combined_ALLNORM_v2 <- subremoved_gates_temp_NPI_combined_ALLNORM_v2 %>% 
  filter(subj_id < 35000) %>%
  select(- ZYawn, - ZPlay, - ZInspect, - ZGroom, - ZManipulate_Object, - ZInspect_Object, - ZFear_Grimace)

View(pr50_gates_temp_NPI_combined_ALLNORM_v2)
 
pr33_gates_temp_NPI_combined_ALLNORM_v2 <- subremoved_gates_temp_NPI_combined_ALLNORM_v2 %>%
  filter(subj_id > 35000)

#View(pr33_gates_temp_NPI_combined_ALLNORM_v2)

```

```{r}
year1_ttests <- data.frame(lapply(pr50_gates_temp_NPI_combined_ALLNORM_v2[,4:98], function(i) t.test(i ~ pr50_gates_temp_NPI_combined_ALLNORM_v2$group)$p.value))
year1_ttests$behavior <- rownames(year1_ttests) # edit: new column for protein_name 
rownames(year1_ttests) <- NULL           # edit: new column for protein_name
year1_ttests[year1_ttests$p.value < 0.05/105,]

```


```{r}
###creating composites from the temperament data
subremoved_gates_temp_NPI_combined_ALLNORM_v2$interact_novelty <- 
  (subremoved_gates_temp_NPI_combined_ALLNORM_v2$ZInspect_Snake_temperament +  
  subremoved_gates_temp_NPI_combined_ALLNORM_v2$ZInteract_Snake_temperament + 
  subremoved_gates_temp_NPI_combined_ALLNORM_v2$ZInteract_Novel_Food_temperament)
  
subremoved_gates_temp_NPI_combined_ALLNORM_v2$ZBark_temperament + subremoved_gates_temp_NPI_combined_ALLNORM_v2$ZShriek_temperament + subremoved_gates_temp_NPI_combined_ALLNORM_v2$temera


subremoved_gates_temp_NPI_combined_ALLNORM_v2$ <- subremoved_gates_temp_NPI_combined_ALLNORM_v2$ZScratch_temperament + subremoved_gates_temp_NPI_combined_ALLNORM_v2$ZTeeth_Grind_temperament + subremoved_gates_temp_NPI_combined_ALLNORM_v2$ZSelf_Groom_temperament + subremoved_gates_temp_NPI_combined_ALLNORM_v2$tempera
subremoved_gates_temp_NPI_combined_ALLNORM_v2$anxious_temp <- 
subremoved_gates_temp_NPI_combined_ALLNORM_v2$_temp <-  



```

```{r differences in FA}



##lateral ventricles omitted
FA_vals_withgroup <- subremoved_gates_temp_NPI_combined_ALLNORM_v2 %>%
  select(subj_id, group, Zfrontal_FA_y2_4mo, Zfrontal_wm_FA_y2_4mo, Zsensory_motor_FA_y2_4mo, Ztemporal_FA_y2_4mo, Ztemporal_wm_FA_y2_4mo, Zparietal_FA_y2_4mo, Zoccipital_FA_y2_4mo, Zoccipital_wm_FA_y2_4mo, Zallocortex_FA_y2_4mo, Zcortex_skeketon_FA_y2_4mo, Zhippocampus_FA_y2_4mo, Zinternal_capsule_wm_FA_y2_4mo, Zexternal_extreme_capsule_wm_FA_y2_4mo, Zslemium_wm_FA_y2_4mo, Zbody_wm_FA_y2_4mo, Zgenu_wm_FA_y2_4mo, Zcerebral_wm_FA_y2_4mo) 

FA_vals_withgroup$group <- as.factor(FA_vals_withgroup$group)
FA_vals_withgroup$subj_id <- as.factor(FA_vals_withgroup$subj_id)


View(FA_vals_withgroup)

occipital_ttest <- t.test(FA_vals_withgroup$Zoccipital_FA_y2_4mo ~ FA_vals_withgroup$group)


melted_FA_vals <- melt(FA_vals_withgroup)

View(melted_FA_vals)

ggplot(melted_FA_vals, aes(x=group, y=value, color = group)) +
geom_boxplot() + geom_jitter() +
  facet_wrap(~variable, ncol = 3) + 
  labs(title = "FA Differences; all years")

skim(FA_vals_withgroup)

shapiro.test(FA_vals_withgroup$Zfrontal_wm_FA_y2_4mo)

frontal_wm_FA_ttest <- t.test(FA_vals_withgroup$Zfrontal_wm_FA_y2_4mo~FA_vals_withgroup$group)

hippcampal_wm_FA_ttest <- t.test(FA_vals_withgroup$Zhippocampus_FA_y2_4mo~FA_vals_withgroup$group)

temporal_wm_FA_ttest <- t.test(FA_vals_withgroup$Ztemporal_wm_FA_y2_4mo~FA_vals_withgroup$group)

temporal_FA_ttest <- t.test(FA_vals_withgroup$Ztemporal_FA_y2_4mo~FA_vals_withgroup$group)

sensorymotor_wm_FA_ttest <- t.test(FA_vals_withgroup$Zsensory_motor_FA_y2_4mo~FA_vals_withgroup$group)

occipital_wm_FA_ttest <- t.test(FA_vals_withgroup$Zoccipital_wm_FA_y2_4mo~FA_vals_withgroup$group)

occipital_FA_ttest <- t.test(FA_vals_withgroup$Zoccipital_FA_y2_4mo~FA_vals_withgroup$group)


occipital_wm_FA_ttest$p.value
#occipital not significant

occipital_FA_ttest$p.value
#not significant

sensorymotor_wm_FA_ttest$p.value
#not significant

temporal_FA_ttest$p.value
#not significant

temporal_wm_FA_ttest$p.value
#.059

hippcampal_wm_FA_ttest$p.value


_wm_FA_ttest$p.value

frontal_wm_FA_ttest$p.value
```

```{r}
temperament_ttests_data <- subremoved_gates_temp_NPI_combined_ALLNORM_v2 %>%
  select(subj_id, group, matches("temperament")) 

temperament_ttests_data$subj_id <- as.factor(temperament_ttests_data$subj_id)
temperament_ttests_data$group <- as.factor(temperament_ttests_data$group)

melted_temp_ttest_data <- melt(temperament_ttests_data)

View(melted_temp_ttest_data)

skim(temperament_ttests_data)

ggplot(melted_temp_ttest_data, aes(x=group, y=value, color = group)) +
geom_boxplot() + geom_jitter() +
  facet_wrap(~variable, ncol = 3) + 
  labs(title = "Behavioral Differences related to Temperament Testing")
##from the output above, it appears that stationary front, grunt, bark, and movement are the behaviors 

##below then we will look at ttest with that data 
```

```{r}
stationary_front_ttest <- t.test()
grunt_ttest <- 
bark_ttest <- wilcox.test()
movement_ttest <- 
```


```{r}



```