npi_visuals.Rmd

---
title: "npi_firstpass"
author: "aj_mitchell"
date: "April 3, 2019"
output: 
  html_document:
    toc: true
    toc_float: true
    theme: yeti
---

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


```{r, echo=FALSE, warning=FALSE, message=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(readxl)
library(manipulate)
library(lavaan)
library(semPlot)
library(readxl)
library(formattable)
library(DT)
library(ggthemes)
library(ggcorrplot)
library(wesanderson)
library(RColorBrewer)
#library(ggthemes)
```

**This document is for explaining the pilot NPI behavior analysis.**

##Description##

A total of 21 subjects were included in the analysis. 

30 individual behaviors were included in the analysis and subsequently divided up into 4 composite groups to investigate the effects of PR. These composite groups were made based on aggregates of **anxious behaviors**, **pro-social behavior (social affiliation)**, **non-affiliative social behaviors**, and **object oriented behaviors**.  These scores are outlined in later tables below.  

The table below provides correlational values about the individual behaviors.

```{r, echo=FALSE, warning=FALSE, message=FALSE}
    
NPI_znorm_v1 <- read_excel("norm_correct_avgdsubj_NPI_v1.xlsx")
NPI_znorm_v1$diet_group <- as.factor(NPI_znorm_v1$diet_group)

#skim(selected_NPIgroup_znorm_v1)

selected_NPI_znorm_v1 <- NPI_znorm_v1 %>%
  select(- Result_Containers, - Statistics, - Observations, - Intervals, - Subjects, - Animal_lD, - Hanging_Toy_at_perch_N0Y1, - Perch_at_dividerN0Y1, - Duration, - Introduction, - diet_group)

#skim(selected_NPI_znorm_v1)

selected_NPIgroup_znorm_v1 <- NPI_znorm_v1 %>%
  select(- Result_Containers, - Statistics, - Observations, - Intervals, - Subjects, - Hanging_Toy_at_perch_N0Y1, - Perch_at_dividerN0Y1, - Duration, - Introduction)
```


The subsequent table outlines the overall effects of protein restriction whereas the latter graphs detail specific differences of either 33% or 50% PR. I did not look at individual behavioral differences (comparing individual behaviors across diet groups). **Only pro-social behavior surrounding the 50 PR showed significant differences based on diet group.** This is visualized later in the document. 

The 'v2' just capture the little changes that I thought improved the groups. 



```{r, include=F}
#View(selected_NPIgroup_znorm_v1)

#ListVariableNames <- function(DATA) {
 #   cat('"', noquote(paste(names(DATA), collapse='","')), '"\n', sep='')
#}


#anx_behavs <- c("ZStereotypy",
#               "ZScratch",
 ##              "ZYawn",
   #            "ZCrook_Tail",
    #           "ZPlay",
     #          "ZPlay_Face",
      #         "ZNo_Response",
       #        "ZLipsmack",
        #       "ZInspect",
         #      "ZMouth_peer",
          #     "ZContact_peer",
           #    "ZContact_peer_count",
            #  "ZOpen_Mouth_Threat",
             #  "ZAggress",
              # "ZFear_Grimace",
               #"ZBite",
               #"ZSteal_Resource_Attempt",
               #"ZSteal_Resource_Attempt_count",
               #"ZResource_Takeover",
               #"ZDisplace",
               #"ZAvoid",
               #"ZFollow",
               #"ZGroom",
               #"ZApproach_peer",
               #"ZLeave_peer",
               #"ZNo_Response_to_Object",
               #"ZMouth_Object",
               #"ZTouch_Object",
               #"ZManipulate_Object",
               #"ZInspect_Object")
```               

**Prosocial** = "Play",
               "Play Face",
               "Contact peer",
               "Follow",
               "Groom",
               "Approach peer",
               "Lipsmack"
               
**Prosocial_v2** = "Play",
               "Play Face",
               "Contact peer",
               "Follow",
               "Groom",
               "Approach peer",
               "Lipsmack",
               "Inspect",
               "Mouth peer"

**Nonsocial** = "Open Mouth Threat",
               "Aggress",
               "Bite",
               "Displace",
               "Avoid",
               "Crooked Tail",
               
**Nonsocial_v2** = "Open Mouth Threat",
               "Aggress",
               "Bite",
               "Resource Takeover",
               "Displace",
               "Avoid",
               "Crooked Tail",
               "Steal Resource Attempt Count"

**Object_Oriented** = "Mouth Object",
               "Touch Object",
               "Manipulate Object",
               "Inspect Object",
               "Steal Resource Attempt"

**Anxious** =  "Stereotypy",
               "Scratch",
               "Yawn",
               "Avoid",
               "Fear Grimmace"
               
**Omitted** = "No Response to Object",
              "Contact Peer Count",
              "No Response"

               
```{r, echo=F, warning=F, message=F}
#Behavior_table <- rbind(Anxious, Prosocial, Nonsocial, Object_Oriented)

#ListVariableNames(selected_NPIgroup_znorm_v1)
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}
cor_selected_NPI <- cor(selected_NPI_znorm_v1)
melted_cor_selected_NPI <- melt(cor_selected_NPI) %>%
  filter(value <= .99)
```

##Correllation Values and Matrix##

```{r, echo=F, message=F, warning=F}
datatable(formattable(melted_cor_selected_NPI))
```


```{r, echo=F, message=F, warning=F}
library(wesanderson)
corrmat_NPI <- ggplot(melted_cor_selected_NPI, aes(x=Var2, y=Var1, fill=value)) +
  geom_tile() +
  scale_fill_gradient(low = "dark blue", high = "orange") + 
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90,
                                   vjust = 1,
                                   size = 9,
                                   hjust = 1))

corrmat_NPI+labs(title="NPI Behavior Correlation Matrix", x="Diet")  +
  scale_color_manual(values = wes_palette(n=2, "GrandBudapest1")) +
  theme(plot.background = element_rect(fill = "light gray"),
        panel.background = element_rect(fill = "dark gray"),
        legend.background = element_rect(fill = "dark gray"))

```


```{r, echo=FALSE, warning=FALSE, message=FALSE}

############
selected_NPIgroup_znorm_v1$social_affiliative <- selected_NPIgroup_znorm_v1$ZContact_peer + selected_NPIgroup_znorm_v1$ZPlay + selected_NPIgroup_znorm_v1$ZPlay_Face + selected_NPIgroup_znorm_v1$ZGroom + selected_NPIgroup_znorm_v1$ZMouth_peer + selected_NPIgroup_znorm_v1$ZLipsmack + selected_NPIgroup_znorm_v1$ZApproach_peer 


#v2 are edits to my groups based on what is outlined in the above text. social_affiliative_v2 adds follow and inspect to the original prosocial group 
selected_NPIgroup_znorm_v1$social_affiliative_v2 <- selected_NPIgroup_znorm_v1$ZContact_peer + selected_NPIgroup_znorm_v1$ZPlay + selected_NPIgroup_znorm_v1$ZPlay_Face + selected_NPIgroup_znorm_v1$ZGroom + selected_NPIgroup_znorm_v1$ZMouth_peer + selected_NPIgroup_znorm_v1$ZLipsmack + selected_NPIgroup_znorm_v1$ZApproach_peer + selected_NPIgroup_znorm_v1$ZFollow + selected_NPIgroup_znorm_v1$ZInspect

############

#count not used because object orientation is how long were you attempting to steal the resource, not how much?
selected_NPIgroup_znorm_v1$object_oriented <- selected_NPIgroup_znorm_v1$ZInspect_Object + selected_NPIgroup_znorm_v1$ZMouth_Object + selected_NPIgroup_znorm_v1$ZManipulate_Object + selected_NPIgroup_znorm_v1$ZTouch_Object + selected_NPIgroup_znorm_v1$ZSteal_Resource_Attempt 

##############

selected_NPIgroup_znorm_v1$social_nonaffiliative <- selected_NPIgroup_znorm_v1$ZCrook_Tail + selected_NPIgroup_znorm_v1$ZOpen_Mouth_Threat + selected_NPIgroup_znorm_v1$ZAggress + selected_NPIgroup_znorm_v1$ZBite + selected_NPIgroup_znorm_v1$ZDisplace + selected_NPIgroup_znorm_v1$ZLeave_peer + selected_NPIgroup_znorm_v1$ZAvoid

# v2 nonaffiliative_v2 has resource takeover and resource takeover attempts (count) included. count is used becuase it is a measurement of how persistent an individual was at taking over a resource
selected_NPIgroup_znorm_v1$social_nonaffiliative_v2 <- selected_NPIgroup_znorm_v1$ZCrook_Tail + selected_NPIgroup_znorm_v1$ZOpen_Mouth_Threat + selected_NPIgroup_znorm_v1$ZAggress + selected_NPIgroup_znorm_v1$ZBite + selected_NPIgroup_znorm_v1$ZDisplace + selected_NPIgroup_znorm_v1$ZLeave_peer + selected_NPIgroup_znorm_v1$ZAvoid + selected_NPIgroup_znorm_v1$ZResource_Takeover + selected_NPIgroup_znorm_v1$ZSteal_Resource_Attempt_count


################

selected_NPIgroup_znorm_v1$anxious <- selected_NPIgroup_znorm_v1$ZStereotypy + selected_NPIgroup_znorm_v1$ZScratch + selected_NPIgroup_znorm_v1$ZYawn + selected_NPIgroup_znorm_v1$ZAvoid + selected_NPIgroup_znorm_v1$ZFear_Grimace

#View(selected_NPIgroup_znorm_v1)
```

```{r, echo=FALSE, message=FALSE, warning=FALSE}
#740 and 764 reintroduced


######## object oriented

object_oriented_v <- selected_NPIgroup_znorm_v1 %>%
  select(Animal_lD, object_oriented, diet_group) 

#..vR == the undefined diet group subjects (0, n=2) removed 
object_oriented_vR <- selected_NPIgroup_znorm_v1 %>%
  select(Animal_lD, object_oriented, diet_group)


object_oriented_vR_50 <- object_oriented_vR %>%
  filter(Animal_lD <= 35000)

object_oriented_vR_33 <- object_oriented_vR %>%
  filter(Animal_lD >= 35000)

#glimpse(object_oriented_vR_50)

####### social affiliative

social_affiliative_v <- selected_NPIgroup_znorm_v1 %>%
  select(Animal_lD, social_affiliative, diet_group)


social_affiliative_vR <- selected_NPIgroup_znorm_v1 %>%
  select(Animal_lD, social_affiliative, diet_group) 

social_affiliative_vR_v2 <- selected_NPIgroup_znorm_v1 %>%
  select(Animal_lD, social_affiliative_v2, diet_group)


social_affiliative_vR_50 <- social_affiliative_vR %>%
  filter(Animal_lD <= 35000)

social_affiliative_vR_33 <- social_affiliative_vR %>%
  filter(Animal_lD >= 35000)

social_affiliative_vR_50_v2 <- social_affiliative_vR_v2 %>%
  filter(Animal_lD <= 35000)

social_affiliative_vR_33_v2 <- social_affiliative_vR_v2 %>%
  filter(Animal_lD >= 35000)


#View(social_affiliative_vR_33)

######## social_nonaffiliative 

social_nonaffiliative_v <- selected_NPIgroup_znorm_v1 %>%
  select(Animal_lD, social_nonaffiliative, diet_group)

social_nonaffiliative_vR <- selected_NPIgroup_znorm_v1 %>%
  select(Animal_lD, social_nonaffiliative, diet_group) 

social_nonaffiliative_vR_v2 <- selected_NPIgroup_znorm_v1 %>%
  select(Animal_lD, social_nonaffiliative_v2, diet_group) 

social_nonaffiliative_vR_50 <- social_nonaffiliative_vR %>%
  filter(Animal_lD <= 35000)

social_nonaffiliative_vR_33 <- social_nonaffiliative_vR %>%
  filter(Animal_lD >= 35000)

#changes below are made to compare and include the altered versions of the composite groups 

social_nonaffiliative_vR_50_v2 <- social_nonaffiliative_vR_v2 %>%
  filter(Animal_lD <= 35000)

social_nonaffiliative_vR_33_v2 <- social_nonaffiliative_vR_v2 %>%
  filter(Animal_lD >= 35000)

#skim(social_nonaffiliative_vR_33)

######## anxious 

anxious_v <- selected_NPIgroup_znorm_v1 %>%
  select(Animal_lD, anxious, diet_group)

anxious_vR <- selected_NPIgroup_znorm_v1 %>%
  select(Animal_lD, anxious, diet_group) 

anxious_vR_50 <- anxious_vR %>%
  filter(Animal_lD <= 35000)

anxious_vR_33 <- anxious_vR %>%
  filter(Animal_lD >= 35000)

#skim(anxious_vR_33)

```





##Overall Behavioral Differences##

This table details the different effects of diet based on overall year1 and year2 PR vs Control performance.

```{r, echo=F, message=F, warning=F}
#these are totals diet group differences; no splitting between 50 and 33% PR)
ttest_anxious_vR <- wilcox.test(anxious_vR$anxious~anxious_vR$diet_group)
ttest_social_nonaffiliative_vR <- wilcox.test(social_nonaffiliative_vR$social_nonaffiliative~social_nonaffiliative_vR$diet_group)
ttest_social_affiliative_vR <- wilcox.test(social_affiliative_vR$social_affiliative~social_affiliative_vR$diet_group)
ttest_object_oriented_vR <- t.test(object_oriented_vR$object_oriented~object_oriented_vR$diet_group)
#no significant effects of subjective groupings
#ttest_social_nonaffiliative_vR$p.value

pvals_NPI_subjectivegroups <- c(ttest_anxious_vR$p.value, ttest_social_affiliative_vR$p.value, ttest_social_nonaffiliative_vR$p.value, ttest_object_oriented_vR$p.value)
pvals_NPI_titles <- c("Anxious Behavior",
                      "Pro-Social Behavior",
                      "Non-Social Behavior",
                      "Object Oriented Behavior")
table_NPI_subjectivegroups <- data.frame(pvals_NPI_titles, pvals_NPI_subjectivegroups)

datatable(formattable(table_NPI_subjectivegroups))
```



##Year 1 vs Year 2 differences##

The table below describes differences based on year one and year two. In summary, only year 1 displayed behavioral differences. Year 1 subjected to a PR diet were less likely to engage in prosocial behavior. Overall diet differences were not seen, also since the diet was changed between year 1 and year 2, I felt separating the years was more informative. Especially given the results. 

Year 1 Information:
  - n = 8
  - 4 CTR(1) / 4 PR(2)
  
Year 2 Information:
  - n = 13
  - 7 CTR(1) / 6 PR(2)
  
  
Each group was subject to a Shapiro Wilks test for normality, and the respective t-test was done to calculate the p-value. For groupings whose distribution significantly deviated from the norm, a Wilcoxon/Mann-Whitney test was used.

**50per = year 1 @ 50% PR**


**33per = year 2 @ 33% PR**

```{r, echo=FALSE, warning=FALSE, message=FALSE}
#differences dependent on PR levels 

#shapiro.test(social_affiliative_vR_33_v2$social_affiliative_v2)


##############

#anxious..50_ttest has normal distribution 
anxious_vR_33_ttest <- wilcox.test(anxious_vR_33$anxious~anxious_vR_33$diet_group, paired=FALSE)
anxious_vR_50_ttest <- t.test(anxious_vR_50$anxious~anxious_vR_50$diet_group)

##############

#social_affiliative...50 has normal distribution  
social_affiliative_vR_33_ttest <- wilcox.test(social_affiliative_vR_33$social_affiliative~social_affiliative_vR_33$diet_group, paired=FALSE)

social_affiliative_vR_33_v2_ttest <- wilcox.test(social_affiliative_vR_33_v2$social_affiliative_v2~social_affiliative_vR_33_v2$diet_group, paired=FALSE)

social_affiliative_vR_50_ttest <- t.test(social_affiliative_vR_50$social_affiliative~social_affiliative_vR_50$diet_group, paired=FALSE)

social_affiliative_vR_50_v2_ttest <- t.test(social_affiliative_vR_50_v2$social_affiliative_v2~social_affiliative_vR_50_v2$diet_group, paired=FALSE)

#social_affiliative_vR_33_ttest_unpaired$p.value

##############


#both social_non...33 and social_non..50 have normal distributions
social_nonaffiliative_vR_33_ttest <- t.test(social_nonaffiliative_vR_33$social_nonaffiliative~social_nonaffiliative_vR_33$diet_group, paired=F)

social_nonaffiliative_vR_33_v2_ttest <- wilcox.test(social_nonaffiliative_vR_33_v2$social_nonaffiliative_v2~social_nonaffiliative_vR_33_v2$diet_group, paired=F)

social_nonaffiliative_vR_50_ttest <- t.test(social_nonaffiliative_vR_50$social_nonaffiliative~social_nonaffiliative_vR_50$diet_group, paired=F)

social_nonaffiliative_vR_50_v2_ttest <- t.test(social_nonaffiliative_vR_50_v2$social_nonaffiliative_v2~social_nonaffiliative_vR_50_v2$diet_group, paired=F)

##############

#both object_ori..33 and object_ori..50 have normal distributions
object_oriented_vR_33_ttest <- t.test(object_oriented_vR_33$object_oriented~object_oriented_vR_33$diet_group, paired=F)

object_oriented_vR_50_ttest <- t.test(object_oriented_vR_50$object_oriented~object_oriented_vR_50$diet_group, paired=F)


###########
#table
###########



pvals_PRgroups_on_behaviors <- c(anxious_vR_33_ttest$p.value,
                                 anxious_vR_50_ttest$p.value,
                                 social_affiliative_vR_33_ttest$p.value,
                                 social_affiliative_vR_33_v2_ttest$p.value,
                                 social_affiliative_vR_50_ttest$p.value,
                                 social_affiliative_vR_50_v2_ttest$p.value,
                                 social_nonaffiliative_vR_33_ttest$p.value,
                                 social_nonaffiliative_vR_33_v2_ttest$p.value,
                                 social_nonaffiliative_vR_50_ttest$p.value,
                                 social_nonaffiliative_vR_50_v2_ttest$p.value,
                                 object_oriented_vR_33_ttest$p.value,
                                 object_oriented_vR_50_ttest$p.value)
pvals_PRgroups_titles <- c("Anxious_33per_PR",
                           "Anxious_50per_PR",
                           "Pro-Social_33per_PR",
                           "Pro-Social_33per_PR_v2",
                           "Pro-Social_50per_PR",
                           "Pro-Social_50per_PR_v2",
                           "Non-Affiliative_33per_PR",
                           "Non-Affiliative_33per_PR_v2",
                           "Non-Affiliative_50per_PR",
                           "Non-Affiliative_50per_PR_v2",
                           "Object_Oriented_33per_PR",
                           "Object_Oriented_50per_PR")
#this table defines differences based on 50 percent and 30 percent PR..
PR_group_differences <- data.frame(pvals_PRgroups_titles, pvals_PRgroups_on_behaviors)

datatable(formattable(PR_group_differences))
#prosocial behaviors are significan
```

Below is the boxplot overlayed with scatter plot showing both the avg difference between groups and where each individual falls per group. **This was the only behavior category that showed significant differences based on diet.**

```{r, echo=FALSE, message=FALSE, warning=FALSE}
#library(wesanderson)
library(ggthemes)

customviz_affil_vR50 <- ggplot(social_affiliative_vR_50, aes(x=diet_group, y=social_affiliative, color=diet_group)) +
  geom_boxplot() +
  geom_jitter() +
  theme_grey()
  #geom_smooth(method = "lm", se =FALSE)

customviz_affil_vR50+labs(title="Year 1: Prosocial Behaviors** (CTR=1, PR=2)", x="Diet",y="Affiliative Behavior Score")  +
  scale_color_manual(values = wes_palette(n=2, "GrandBudapest1")) +
  theme(plot.background = element_rect(fill = "light gray"),
        panel.background = element_rect(fill = "dark gray"),
        legend.background = element_rect(fill = "dark gray"))

#social_aff_color_test+labs(title="Social Affiliative Behaviors", x="Diet",y="Affiliative Composite Score") +
  #theme_economist()
```


the table below is the output when you control for housing. After speaking with Jacqui, she said that several animals in year 2 were housed differently than the others. Thus, i removed these animals to look at the **effect housing has on the beahvioral outcome**.

Take this data with a grain of salt though.. at this point, a lot of manipulation have been made and i will likely need to reevaluate this part of the analysis. Additionally, it has not been looked at as thoroughly as other measurements. 

The boxplot below is for fact checking the distribution

```{r, echo=F, warning=F, message=F}

########
#cutting out the outlier; only applies to year 2
########
#controling for housing <- select out 35375, 35455, 35476, 35712

matched_housing_object_ori <- object_oriented_vR %>%
  filter(!Animal_lD %in% 35375, 35455, 35476, 35712)

#View(matched_housing_object_ori)

matched_housing_prosocial_v1 <- social_affiliative_vR %>%
  filter(!Animal_lD %in% 35375, 35455, 35476, 35712)

matched_housing_prosocial_v2 <- social_affiliative_vR_33 %>%
  filter(!Animal_lD %in% 35375, 35455, 35476, 35712)

matched_housing_prosocial_v2_ttesst <- wilcox.test(matched_housing_prosocial_v2$social_affiliative~matched_housing_prosocial_v2$diet_group)
#matched_housing_prosocial_v2_ttesst$p.value

matched_housing_nonsocial_v1 <- social_nonaffiliative_vR %>%
  filter(!Animal_lD %in% 35375, 35455, 35476, 35712)

#matched_housing_nonsocial_v2 <- social_nonaffiliative_vR %>%
  #filter(!Animal_lD %in% 35375, 35455, 35476, 35712)

matched_housing_anxious <- anxious_vR %>%
  filter(!Animal_lD %in% 35375, 35455, 35476, 35712)

#only object has normal distrib
#shapiro.test(matched_housing_prosocial_v2$social_affiliative)

matched_housing_prosocial_v1_ttest <- wilcox.test(matched_housing_prosocial_v1$social_affiliative~matched_housing_prosocial_v1$diet_group, paired=F)
matched_housing_nonsocial_v1_ttest <- wilcox.test(matched_housing_nonsocial_v1$social_nonaffiliative~matched_housing_nonsocial_v1$diet_group, paired=F)
matched_housing_object_ori_ttest <- t.test(matched_housing_object_ori$object_oriented~matched_housing_object_ori$diet_group, paired=F)
matched_housing_anxious_ttest <- wilcox.test(matched_housing_anxious$anxious~matched_housing_anxious$diet_group, paired=F)


matched_housing_pvals <- c(matched_housing_prosocial_v1_ttest$p.value, matched_housing_nonsocial_v1_ttest$p.value, matched_housing_object_ori_ttest$p.value, matched_housing_anxious_ttest$p.value)
matched_housing_titles <- c("Prosocial Controlled Housing",
                            "Nonsocial Controlled Housing",
                            "Object Oriented Controlled Housing",
                            "Anxious Controlled Housing")

table_matched_housing <- data.frame(matched_housing_titles, matched_housing_pvals)

datatable(formattable(table_matched_housing))


customviz_matched <- ggplot(matched_housing_nonsocial_v1, aes(x=diet_group, y=social_nonaffiliative, color=diet_group)) +
  geom_boxplot() +
  geom_jitter() +
  theme_grey()
  #geom_smooth(method = "lm", se =FALSE)

customviz_matched+labs(title="Non-Social Behavior_HousMath (CTR=1, PR=2)", x="Diet",y="Non-Affiliative Behavior Score")  +
  scale_color_manual(values = wes_palette(n=2, "FantasticFox1")) +
  theme(plot.background = element_rect(fill = "light gray"),
        panel.background = element_rect(fill = "dark gray"),
        legend.background = element_rect(fill = "dark gray"))
```



##Extra Visualizations##

Here, only groups i thought were interesting are visualized; just for simplicity. We can add as many as desired. 

**Non-affiliative behaviors during year 1 (50% PR) were strongly trending but not significant.**

Refer to the list in the **Description Section** for individual behaviors that make up this group. 

```{r, echo=F, warning=F, message=F}
####nonaffiliative
#v2 includes resource takeover attempts 

customviz_nonaffil_vR50 <- ggplot(social_nonaffiliative_vR_50, aes(x=diet_group, y=social_nonaffiliative, color=diet_group)) +
  geom_boxplot() +
  geom_jitter() +
  theme_grey()
  #geom_smooth(method = "lm", se =FALSE)

customviz_nonaffil_vR50+labs(title="Year 1: Social Non-Affiliative Behaviors (CTR=1, PR=2)", x="Diet",y="Non-Affiliative Behavior Score")  +
  scale_color_manual(values = wes_palette(n=2, "Royal1")) +
  theme(plot.background = element_rect(fill = "light gray"),
        panel.background = element_rect(fill = "dark gray"),
        legend.background = element_rect(fill = "dark gray"))

customviz_nonaffil_vR50_v2 <- ggplot(social_nonaffiliative_vR_50_v2, aes(x=diet_group, y=social_nonaffiliative_v2, color=diet_group)) +
  geom_boxplot() +
  geom_jitter() +
  theme_grey()
  #geom_smooth(method = "lm", se =FALSE)

customviz_nonaffil_vR50_v2+labs(title="Year 1: Social Non-Affiliative (v2) Behaviors (CTR=1, PR=2)", x="Diet",y="Non-Affiliative Behavior (v2) Score")  +
  scale_color_manual(values = wes_palette(n=2, "Royal1")) +
  theme(plot.background = element_rect(fill = "light gray"),
        panel.background = element_rect(fill = "dark gray"),
        legend.background = element_rect(fill = "dark gray"))


#######prosocial: year 1
#v2 includes follow and 


customviz_affil_vR50 <- ggplot(social_affiliative_vR_50, aes(x=diet_group, y=social_affiliative, color=diet_group)) +
  geom_boxplot() +
  geom_jitter() +
  theme_grey()
  #geom_smooth(method = "lm", se =FALSE)

customviz_affil_vR50+labs(title="Year 1: Prosocial Behaviors (CTR=1, PR=2)", x="Diet",y="Prosocial Behavior Score")  +
  scale_color_manual(values = wes_palette(n=2, "Cavalcanti1")) +
  theme(plot.background = element_rect(fill = "light gray"),
        panel.background = element_rect(fill = "dark gray"),
        legend.background = element_rect(fill = "dark gray"))


customviz_affil_vR50_v2 <- ggplot(social_affiliative_vR_50_v2, aes(x=diet_group, y=social_affiliative_v2, color=diet_group)) +
  geom_boxplot() +
  geom_jitter() +
  theme_grey()
  #geom_smooth(method = "lm", se =FALSE)

customviz_affil_vR50+labs(title="Year 1: Prosocial (v2) Behaviors (CTR=1, PR=2)", x="Diet",y="Prosocial (v2) Behavior Score")  +
  scale_color_manual(values = wes_palette(n=2, "Cavalcanti1")) +
  theme(plot.background = element_rect(fill = "light gray"),
        panel.background = element_rect(fill = "dark gray"),
        legend.background = element_rect(fill = "dark gray"))



```


**Closing**

I'd like to talk about adjusting the groupings and how we could expand on that end. Or would these be the groups we move forward with? I thought they were useful because the interpretation of them is more simple than the community detection approach.


There are many things I need to improve on here but mainly:

-  **I want to finalize and better form my groupings** 
      - some of the composites were created using steal_resource_attempt rather than steal_resource_attempt_count.
And i am unsure of which would be the better approach for quantifying such behaviors. 
      
-  **I plan to run an SEM analysis to look at the fit of these composites** 
      - I wanted to get this to you prior to doing that so we can discuss changes/approaches to the SEMs.
      
      
-  **I am open to any changes. I am unsure if the makeup of the data allows us to use Community Detection, and I also thought grouping terms into pro-social/non-social/anxious provides a better output.** 
      - But again, we can adjust where needed.. this adjustment/documentation process has become increasingly faster for me. 
    

##new graphs and tables for NPI meeting##

```{r, echo=FALSE, warning=F, message=F}
##importing data

year1_4mo_NPI_combined <- read_excel("4mo_year1_NPI_combined.xlsx")
#str(year1_4mo_NPI_combined)

year2_4mo_NPI_combined <- read_excel("4mo_year2_NPI_combined.xlsx")
#View(year2_4mo_NPI_combined)

year2_11mo_NPI_combined <- read_excel("11mo_year2_NPI_combined.xlsx")
#str(year2_11mo_NPI_combined)

year2_4mo_11mo_NPI_combined <- read_excel("year2_4mo_11mo_combined.xlsx")
#View(year2_4mo_11mo_NPI_combined)

yr1_yr2_4mo_NPI_combined <- read_excel("yr1_yr2_4mo_combined.xlsx")
#View(yr1_yr2_4mo_NPI_combined)
```




```{r, echo=FALSE, warning=F, message=F}
##removing subjects with missing data

selected_yr1_4moNPI_v1 <- year1_4mo_NPI_combined %>%
  filter(!subj_id == 34404,
         !subj_id == 34412)

selected_yr2_4moNPI_v1 <- year2_4mo_NPI_combined %>%
  filter(!subj_id == 35408)

selected_yr2_11moNPI_v1 <- year2_11mo_NPI_combined %>%
  filter(!subj_ID == 35505,
         !subj_ID == 35408)

selected_year2_4mo_11mo_NPI_combined <- year2_4mo_11mo_NPI_combined %>%
  filter(!subj_id == 35505,
         !subj_id == 35408)

selected_yr1_yr2_4mo_NPI_combined <- yr1_yr2_4mo_NPI_combined %>%
  filter(!subj_id == 35505,
         !subj_id == 35408)

####selecting out subj_id and group

selected_yr1_4moNPI_v1_4cor <- year1_4mo_NPI_combined %>%
  filter(!subj_id == 34404,
         !subj_id == 34412) %>%
  select(- subj_id, -group)

selected_yr2_4moNPI_v14cor <- year2_4mo_NPI_combined %>%
  filter(!subj_id == 35408) %>%
  select(- subj_id, -group)

selected_yr2_11moNPI_v1_4cor <- year2_11mo_NPI_combined %>%
  filter(!subj_ID == 35505,
         !subj_ID == 35408) %>%
  select(-subj_ID, -group)

selected_year2_4mo_11mo_NPI_combined_4cor <- year2_4mo_11mo_NPI_combined %>%
  filter(!subj_id == 35505,
         !subj_id == 35408) %>%
  select(- subj_id, -group)

selected_yr1_yr2_4mo_NPI_combined_4cor <- selected_yr1_yr2_4mo_NPI_combined %>%
  filter(!subj_id == 35505,
         !subj_id == 35408, 
         !subj_id == 34412) %>%
  select(- subj_id, -group)


#View(selected_yr2_11moNPI_v1_4cor)
```



```{r, echo=FALSE, warning=F, message=F}
##calculating correlations
cor_yr1_4moNPI_v1 <- cor(selected_yr1_4moNPI_v1_4cor)
melted_cor_yr1_4moNPI_v1 <- melt(cor_yr1_4moNPI_v1)

thresh_melt_cor_yr1 <- melted_cor_yr1_4moNPI_v1 %>%
  filter(!is.na(value),
         value <= .99)

#View(thresh_melt_cor_yr1)


cor_yr2_4moNPI_v1 <- cor(selected_yr2_4moNPI_v14cor)
melted_cor_yr2_4moNPI_v1 <- melt(cor_yr2_4moNPI_v1)

cor_yr2_11moNPI_v1 <- cor(selected_yr2_11moNPI_v1_4cor)
melted_cor_yr2_11moNPI_v1 <- melt(cor_yr2_11moNPI_v1)

cor_yr2_4mo_11mo_NPI_4cor <- cor(selected_year2_4mo_11mo_NPI_combined_4cor)
metled_cor_yr2_4mo_11mo_NPI_4cor <- melt(cor_yr2_4mo_11mo_NPI_4cor)

cor_selected_yr1_yr2_4mo_NPI_combined_4cor <- cor(selected_yr1_yr2_4mo_NPI_combined_4cor)
metled_cor_selected_yr1_yr2_4mo_NPI_combined_4cor <- melt(cor_selected_yr1_yr2_4mo_NPI_combined_4cor)

#View(metled_cor_selected_yr1_yr2_4mo_NPI_combined_4cor)

```

**For the tables below, type in your measurement of interest (e.g. "FA" or "brain_volume") to see how it correlates with other behaviors. You can also type in values (e.g. ".8") to isolate variables with strong correlations. Additionally, you can extend this to ".8 FA to isolate further.** 

```{r, echo=F, warning=F, message=F}
datatable(formattable(melted_cor_yr1_4moNPI_v1))

datatable(formattable(melted_cor_yr2_4moNPI_v1))

datatable(formattable(melted_cor_yr2_11moNPI_v1))

datatable(formattable(metled_cor_yr2_4mo_11mo_NPI_4cor))

```


##new graphs for NPI meeting##

**The graphs below are visualizations for how behaviors and structural measurements relate to each other. The title describe the input data. I was unable to generate a matrix showing how the 4mo_yr1 and 4mo_yr2 data relate to behavior as well as each other, but I should be able to figure this out after the weekend.**


**Additionally, I want to investigate the percent change between yr2_4mo and yr2_11mo subjects related to structural, volume, and diffusion (MRI metrics) and relate these to expressions in behavior. The driving hypothesis being that animals subjected to a protein restricted diet will have less change (flatter slope) between these time points than subjects given a normal diet. I expect PR animals to have decreased changes relating to normal developmental trajectories and relate these differences to differences in behavior between the two diet groups.**

```{r, echo=F, warning=F, message=F}
##yr1_4mo

gplot_yr1_4mo_NPI <- ggplot(melted_cor_yr1_4moNPI_v1, aes(x=Var2, y=Var1, fill=value)) +
  geom_tile() +
  scale_fill_gradientn(colours = topo.colors(10)) + 
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90,
                                   vjust = 1,
                                   size = 9,
                                   hjust = 1))

gplot_yr1_4mo_NPI+labs(title="Year 1 Behavior plus Structure Correlation Matrix", x="Correlation Value", y="Correlation Value")

###yr2_4mo

gplot_yr2_4mo_NPI <- ggplot(melted_cor_yr2_4moNPI_v1, aes(x=Var2, y=Var1, fill=value)) +
  geom_tile() +
  scale_fill_gradientn(colours = topo.colors(10)) + 
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90,
                                   vjust = 1,
                                   size = 9,
                                   hjust = 1))

gplot_yr2_4mo_NPI + labs(title= "Year2_4mo_NPI", x="Correlation Value", y="Correlation Value")

###yr2_11mo

gplot_yr2_11mo_NPI <- ggplot(melted_cor_yr2_11moNPI_v1, aes(x=Var2, y=Var1, fill=value)) +
  geom_tile() +
  scale_fill_gradientn(colours = topo.colors(10)) + 
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90,
                                   vjust = 1,
                                   size = 9,
                                   hjust = 1))

gplot_yr2_11mo_NPI + labs(title="year2_11mo_NPI", x="Correlation Value", y="Correlation Value")

###yr_2_4mo_11mo


gplot_yr2_4mo_11mo_NPI <- ggplot(metled_cor_yr2_4mo_11mo_NPI_4cor, aes(x=Var2, y=Var1, fill=value)) +
  geom_tile() +
  scale_fill_gradientn(colours = topo.colors(10)) + 
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90,
                                   vjust = 1,
                                   size = 9,
                                   hjust = 1,
                                    ))

gplot_yr2_4mo_11mo_NPI + labs(title="Year 2: 4mo plus 11mo Correlations ", x="Correlation Value", y="Correlation Value") 
```


```{r}
##loading the new data

#structural only @ 4mo between yr1 and yr2
##did this to maximize data set - must remove significantly more animals from data when incorporating behavior
structural_yr1_yr2 <- read_excel("structural_yr1_yr2.xlsx")
str(structural_yr1_yr2)

#behavior plus structural @ 4mo between yr1 and yr2
perc_diffs_4mo_yr2_11mo_yr2 <- read_excel("perc_diffs_4mo_yr2_11mo_yr2.xlsx")
View(perc_diffs_4mo_yr2_11mo_yr2)

selected_perc_diffs_4mo_yr2_11mo_yr2 <- perc_diffs_4mo_yr2_11mo_yr2 %>%
  select(- subj_id, - group)

#differences in 4mo_yr1 and 11mo_yr2
##here i'm trying to look at the percent change between the 7 months and see if animals given a normal diet showed greater changes in areas associated with behaviors (i.e are there differences in percent change in areas related to the behavioral deficiets we see - is there a flatter slope in the PFC between diet groups and does their behavior reflect this?)
behav_plus_struc_4moyr1_4moyr2 <- read_excel("behav_plus_struc_4moyr1_4moyr2.xlsx")
str(behav_plus_struc_4moyr1_4moyr2)
```

```{r}

##testing
##select out variables (only in perc_diffs) and calculate correlations as well as pvalues and melt data
##here we are also going to start incorporating the cor.test funciton to get pvalues for Elinor
    #compare matrices to make sure they look the same - problem is i need to make my data in a different format for pvalues 

cor_structural_yr1_yr2 <- cor(structural_yr1_yr2)

cor_pval_structural_yr1_yr2 <- rcorr(as.matrix(structural_yr1_yr2))

cor_structural_yr1_yr2

cor_pval_structural_yr1_yr2$P

melted_cor_structural_yr1_yr2_cor <- melt(cor_structural_yr1_yr2)
melted_cor_structural_yr1_yr2_rcorr <- melt(cor_pval_structural_yr1_yr2$r)

thresh_melted_cor_pval_structural_yr1_yr2 <- melt(cor_pval_structural_yr1_yr2$P) %>%
  filter(value <= .05)

#gplot_thresh_melted_cor_pval_structural_yr1_yr2 <- 
  
  ggplot(melted_cor_structural_yr1_yr2_cor, aes(x=Var2, y=Var1, fill=value)) +
  geom_tile() +
  scale_fill_gradientn(colours = topo.colors(10)) + 
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90,
                                   vjust = 1,
                                   size = 9,
                                   hjust = 1,
                                    ))

gplot_thresh_melted_cor_pval_structural_yr1_yr2 + labs(title="Structural Relationships between Yr1 and Yr1: Correlations", x="P Value", y="P Value") 

```


```{r}
#starting real analysis
cor_pval_structural_yr1_yr2 <- rcorr(as.matrix(structural_yr1_yr2))
melted_cor_structural_yr1_yr2_rcorr <- melt(cor_pval_structural_yr1_yr2$r)

str(melted_cor_structural_yr1_yr2_rcorr)

thresh_melted_cor_pval_structural_yr1_yr2 <- melt(cor_pval_structural_yr1_yr2$P) %>%
  filter(value <= .05)

####

cor_perc_diffs_4mo_yr2_11mo_yr2 <- rcorr(as.matrix(selected_perc_diffs_4mo_yr2_11mo_yr2))
melted_cor_perc_diffs_4mo_yr2_11mo_yr2 <- melt(cor_perc_diffs_4mo_yr2_11mo_yr2$r)
removed_dups_perc_diffs <- melt(cor_perc_diffs_4mo_yr2_11mo_yr2$r) %>%
  filter(value <= 1,
         value >=.25)
thresh_melted_pvals_melted_cor_perc_diffs_4mo_yr2_11mo_yr2 <- melt(cor_perc_diffs_4mo_yr2_11mo_yr2$P) %>%
  filter(value <= .05)

#####

removed_dups_behave_plus_struc <- thresh_melted_pvals_cor_behav_plus_struc_4moyr1_4moyr2[!duplicated(thresh_melted_pvals_cor_behav_plus_struc_4moyr1_4moyr2$value), ]
  
cor_behav_plus_struc_4moyr1_4moyr2 <- rcorr(as.matrix(behav_plus_struc_4moyr1_4moyr2))
melted_cor_behav_plus_struc_4moyr1_4moyr2 <- melt(cor_behav_plus_struc_4moyr1_4moyr2$r)
thresh_melt_behav_struc_test4pval <- melt(cor_behav_plus_struc_4moyr1_4moyr2$r) %>%
  filter(value <= 1,
         value >= .25)
thresh_melted_pvals_cor_behav_plus_struc_4moyr1_4moyr2 <- melt(cor_behav_plus_struc_4moyr1_4moyr2$P) %>%
  filter(value <= .05)

removed_dups_behave_plus_struc <- thresh_melted_pvals_cor_behav_plus_struc_4moyr1_4moyr2[!duplicated(thresh_melted_pvals_cor_behav_plus_struc_4moyr1_4moyr2$value), ]

duplicated(thresh_melted_pvals_cor_behav_plus_struc_4moyr1_4moyr2$value) 
```

```{r}

##Structure Only

#plot of correlation values of structure only
gplot_corplot_structure_corval <- ggcorrplot(cor_structural_yr1_yr2$r, hc.order = TRUE, type = "lower", colors = topo.colors(3)) +  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90,
                                   vjust = 1,
                                   size = 9,
                                   hjust = 1
                                    ))

gplot_corplot_structure_corval + labs(title="Structural Relationships between Yr1 and Yr2: Correlations", x="Correlation Stregnth", y="Correlation Stregnth")

#plot of pvalues for structure only
gplot_corplot_structure_pval <- ggplot(thresh_melted_cor_pval_structural_yr1_yr2, aes(x=Var2, y=Var1, fill=value)) +
  geom_tile() +
  scale_fill_gradientn(colours = topo.colors(5)) + 
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90,
                                   vjust = 1,
                                   size = 9,
                                   hjust = 1
                                    ))

gplot_corplot_structure_pval + labs(title="Structural Relationships between Yr1 and Yr2: Pvalue", x="Significance Stregnth", y="Significance Stregnth")


##this is all structure
datatable(formattable(thresh_melted_cor_pval_structural_yr1_yr2))
```

```{r}

##Structure plus behavior

#plot of correlation values of structure plus behavior
##does not work bc of NaNs
gplot_corplot_behave_struct_corval <- ggcorrplot(cor_behav_plus_struc_4moyr1_4moyr2$r, hc.order = TRUE, type = "lower", colors = topo.colors(3)) +  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90,
                                   vjust = 1,
                                   size = 9,
                                   hjust = 1
                                    ))

gplot_corplot_behav_struct_corval + labs(title="Behavior and Structure Relationships between Yr1 and Yr2: Correlations", x="Correlation Stregnth", y="Correlation Stregnth")

##above does not work bc of NaNs
gplot_behav_plus_struc_4moyr1_4moyr2_v2 <- ggplot(melted_cor_behav_plus_struc_4moyr1_4moyr2, aes(x=Var2, y=Var1, fill=value)) +
  geom_tile() +
  scale_fill_gradientn(colours = topo.colors(9)) + 
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90,
                                   vjust = 1,
                                   size = 9,
                                   hjust = 1
                                    )) +
     theme(axis.text.y = element_text(margin = margin(b = .5, t=.5),
                                   vjust = 1,
                                   size = 9,
                                   hjust = 1
                                    ))

gplot_behav_plus_struc_4moyr1_4moyr2_v2 + labs(title="Year1 and Year2: Structure related to Behavior: Correlation Values", x="Correlation Stregnth", y="Correlation Stregnth")




#plot of pvalues for structure plus behavior
gplot_corplot_behave_structure_pval <- ggplot(thresh_melted_pvals_cor_behav_plus_struc_4moyr1_4moyr2, aes(x=Var2, y=Var1, fill=value)) +
  geom_tile() +
  scale_fill_gradientn(colours = topo.colors(5)) + 
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90,
                                   vjust = 1,
                                   size = 9,
                                   hjust = 1
                                    ))

gplot_corplot_behave_structure_pval + labs(title="Year1 and Year2: Structure related to Behavior: Pvalue", x="Significance Stregnth", y="Significance Stregnth")

datatable(formattable(removed_dups_behave_plus_struc))
```

```{r}
##percent differences between year 2, 4mo and 11mo 

#plot of correlation values of percent differences between year 2, 4mo and 11mo
##DOESNT WORK DUE TO NaNs
gplot_corplot_percdiff_corval <- ggcorrplot(cor_perc_diffs_4mo_yr2_11mo_yr2$r, hc.order = TRUE, type = "lower", colors = topo.colors(3)) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90,
                                   vjust = 1,
                                   size = 9,
                                   hjust = 1
                                    ))

gplot_corplot_structure_corval + labs(title="Year2: Percent Differences 4mo/11mo in Structure related to Behavior: Correlations", x="Correlation Stregnth", y="Correlation Stregnth")

#######above doesnt work bc of NaNs...

gplot_corplot_percdiff_corval_v2 <- ggplot(melted_cor_perc_diffs_4mo_yr2_11mo_yr2, aes(x=Var2, y=Var1, fill=value)) +
  geom_tile() +
  scale_fill_gradientn(colours = topo.colors(10)) + 
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90,
                                   vjust = 1,
                                   size = 9,
                                   hjust = 1
                                    )) +
     theme(axis.text.y = element_text(margin = margin(b = .5, t=.5),
                                   vjust = 1,
                                   size = 9,
                                   hjust = 1
                                    ))

gplot_corplot_percdiff_corval_v2 + labs(title="Year2: Percent Differences in 4mo/11mo Structure related to Behavior: Correlation Values", x="Correlation Stregnth", y="Correlation Stregnth")

#plot of correlation values of percent differences between year 2, 4mo and 11mo
gplot_corplot_cor_percdiff_pval <- ggplot(thresh_melted_pvals_melted_cor_perc_diffs_4mo_yr2_11mo_yr2, aes(x=Var2, y=Var1, fill=value)) +
  geom_tile() +
  scale_fill_gradientn(colours = topo.colors(5)) + 
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90,
                                   vjust = 1,
                                   size = 9,
                                   hjust = 1
                                    )) +
   theme(axis.text.y = element_text(
                                   vjust = 1,
                                   size = 9,
                                   hjust = 1
                                    ))

gplot_corplot_cor_percdiff_pval + labs(title="Year2: Percent Differences in 4mo/11mo Structure related to Behavior: Pvalue", x="Significance Stregnth", y="Significance Stregnth")


datatable(formattable(melted_cor_perc_diffs_4mo_yr2_11mo_yr2))
```




```{r}
##need to go back and look at this with the composite variable included
#additionally need to make tables.. 


##adding the composite


```







![I could not figure out how to get the proportions right in the output of the figure so i took an extended screenshot to represent the last figure; detailing the year2 4 and 11 mo comparison with behavior.](spread_last_fig.PNG)