#Female contraception
library(tidyverse)
xdata <- read_delim(file="chimpanzeeDHEASdata.txt",col_names=T,delim="\t") %>%
  mutate(across(where(is.character), factor))

summary(xdata)
xdata<-droplevels(xdata)

xdata<-subset(xdata, sex== "F")

library(car)
model.test<-lm(dheas~age.serum+contraception, data=xdata)
vif(model.test)
par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(model.test)
trans_cort<-xdata$dheas
trans_model<-model.test

#TESTING ASSUMPTIONS
#Generate residual and predicted values
trans_resids <- residuals(trans_model)
trans_preds <- predict(trans_model)
sq_preds <- trans_preds^2
#Look at a plot of residual vs. predicted values
plot(trans_resids ~ trans_preds, data = xdata,
     xlab = "Predicted Values",
     ylab = "Residuals")
#Perform a Shapiro-Wilk test for normality of residuals
shapiro.test(trans_resids) #significant. Transform
#Perform Levene's Test for homogeneity of variances
leveneTest(trans_cort~contraception,data=xdata, center = mean) 
leveneTest(trans_cort~contraception,data=xdata, center = median)

library(MASS)
b=boxcox(dheas~age.serum+contraception, data=xdata)
lamda=b$x
lik=b$y
bc=cbind(lamda, lik)
bc [order(lik), ]
trans_model<-lm(dheas^0.42~contraception+age.serum, data=xdata)
par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(trans_model)
trans_cort<-xdata$dheas^0.42

#TESTING ASSUMPTIONS
#Generate residual and predicted values
trans_resids <- residuals(trans_model)
trans_preds <- predict(trans_model)
sq_preds <- trans_preds^2
#Look at a plot of residual vs. predicted values
plot(trans_resids ~ trans_preds, data = xdata,
     xlab = "Predicted Values",
     ylab = "Residuals")
#Perform a Shapiro-Wilk test for normality of residuals
shapiro.test(trans_resids)
#Perform Levene's Test for homogeneity of variances
leveneTest(trans_cort~contraception,data=xdata, center = mean) #Significative. Non-homogeneous
leveneTest(trans_cort~contraception,data=xdata, center = median)

library(AICcmodavg)
library(lme4)
m0 <- glm(trans_cort~1, data=xdata, family = gaussian)
m1 <- glm(trans_cort~age.serum*contraception, family=gaussian, data=xdata)
AIC(m0,m1)
m2 <- glm(trans_cort~age.serum+contraception, family=gaussian, data=xdata)
AIC(m2,m0)
summary(m2)
m3 <- glm(trans_cort~age.serum, family=gaussian, data=xdata)
AIC(m3,m0)
m4 <- glm(trans_cort~contraception, family=gaussian, data=xdata)
AIC(m4,m0)

AIC(m0,m1,m2,m3,m4)

library(MuMIn)
model.sel(m0,m1,m2,m3,m4)
summary(m3)

#########CORTISOL
library(car)
model.test<-lm(cortisol~age.serum+contraception, data=xdata)
vif(model.test)
par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(model.test)
trans_cort<-xdata$cortisol
trans_model<-model.test

#TESTING ASSUMPTIONS
#Generate residual and predicted values
trans_resids <- residuals(trans_model)
trans_preds <- predict(trans_model)
sq_preds <- trans_preds^2
#Look at a plot of residual vs. predicted values
plot(trans_resids ~ trans_preds, data = xdata,
     xlab = "Predicted Values",
     ylab = "Residuals")
#Perform a Shapiro-Wilk test for normality of residuals
shapiro.test(trans_resids) #significant. Transform
#Perform Levene's Test for homogeneity of variances
leveneTest(trans_cort~contraception,data=xdata, center = mean) 
leveneTest(trans_cort~contraception,data=xdata, center = median)

library(MASS)
b=boxcox(cortisol~age.serum+contraception, data=xdata)
lamda=b$x
lik=b$y
bc=cbind(lamda, lik)
bc [order(lik), ]
trans_model<-lm(cortisol^0.5~contraception+age.serum, data=xdata)
par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(trans_model)
trans_cort<-xdata$cortisol^0.5

#TESTING ASSUMPTIONS
#Generate residual and predicted values
trans_resids <- residuals(trans_model)
trans_preds <- predict(trans_model)
sq_preds <- trans_preds^2
#Look at a plot of residual vs. predicted values
plot(trans_resids ~ trans_preds, data = xdata,
     xlab = "Predicted Values",
     ylab = "Residuals")
#Perform a Shapiro-Wilk test for normality of residuals
shapiro.test(trans_resids)
#Perform Levene's Test for homogeneity of variances
leveneTest(trans_cort~contraception,data=xdata, center = mean) #Significative. Non-homogeneous
leveneTest(trans_cort~contraception,data=xdata, center = median)

library(AICcmodavg)
library(lme4)
m0 <- glm(trans_cort~1, data=xdata, family = gaussian)
m1 <- glm(trans_cort~age.serum*contraception, family=gaussian, data=xdata)
AIC(m0,m1)
m2 <- glm(trans_cort~age.serum+contraception, family=gaussian, data=xdata)
AIC(m2,m0)
summary(m2)
m3 <- glm(trans_cort~age.serum, family=gaussian, data=xdata)
AIC(m3,m0)
summary(m3)
m4 <- glm(trans_cort~contraception, family=gaussian, data=xdata)
AIC(m4,m0)
summary(m4)

AIC(m0,m1,m2,m3,m4)

library(MuMIn)
model.sel(m0,m1,m2,m3,m4)

#########DHEAS/CORTISOL ratio
library(car)
model.test<-lm(ratio~age.serum+contraception, data=xdata)
vif(model.test)
par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(model.test)
trans_cort<-xdata$ratio
trans_model<-model.test

#TESTING ASSUMPTIONS
#Generate residual and predicted values
trans_resids <- residuals(trans_model)
trans_preds <- predict(trans_model)
sq_preds <- trans_preds^2
#Look at a plot of residual vs. predicted values
plot(trans_resids ~ trans_preds, data = xdata,
     xlab = "Predicted Values",
     ylab = "Residuals")
#Perform a Shapiro-Wilk test for normality of residuals
shapiro.test(trans_resids) #significant. Transform
#Perform Levene's Test for homogeneity of variances
leveneTest(trans_cort~contraception,data=xdata, center = mean) 
leveneTest(trans_cort~contraception,data=xdata, center = median)

library(MASS)
b=boxcox(ratio~age.serum+contraception, data=xdata)
lamda=b$x
lik=b$y
bc=cbind(lamda, lik)
bc [order(lik), ]
trans_model<-lm(ratio^0.18~contraception+age.serum, data=xdata)
par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(trans_model)
trans_cort<-xdata$ratio^0.18

#TESTING ASSUMPTIONS
#Generate residual and predicted values
trans_resids <- residuals(trans_model)
trans_preds <- predict(trans_model)
sq_preds <- trans_preds^2
#Look at a plot of residual vs. predicted values
plot(trans_resids ~ trans_preds, data = xdata,
     xlab = "Predicted Values",
     ylab = "Residuals")
#Perform a Shapiro-Wilk test for normality of residuals
shapiro.test(trans_resids)
#Perform Levene's Test for homogeneity of variances
leveneTest(trans_cort~contraception,data=xdata, center = mean) #Significative. Non-homogeneous
leveneTest(trans_cort~contraception,data=xdata, center = median)

library(AICcmodavg)
library(lme4)
m0 <- glm(trans_cort~1, data=xdata, family = gaussian)
m1 <- glm(trans_cort~age.serum*contraception, family=gaussian, data=xdata)
AIC(m0,m1)
m2 <- glm(trans_cort~age.serum+contraception, family=gaussian, data=xdata)
AIC(m2,m0)
summary(m2)
m3 <- glm(trans_cort~age.serum, family=gaussian, data=xdata)
AIC(m3,m0)
summary(m3)
m4 <- glm(trans_cort~contraception, family=gaussian, data=xdata)
AIC(m4,m0)
summary(m4)

AIC(m0,m1,m2,m3,m4)

library(MuMIn)
model.sel(m0,m1,m2,m3,m4)

#FULL DATA set
#Effects of sex, age, and colony
library(tidyverse)
xdata <- read_delim(file="chimpanzeeDHEASdata.txt",col_names=T,delim="\t") %>%
  mutate(across(where(is.character), factor))

summary(xdata)

library(car)
model.test<-lm(dheas~age.serum+sex+colony, data=xdata)
vif(model.test)

par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(model.test)
trans_cort<-xdata$dheas
trans_model<-model.test

#TESTING ASSUMPTIONS
#Generate residual and predicted values
trans_resids <- residuals(trans_model)
trans_preds <- predict(trans_model)
sq_preds <- trans_preds^2
#Look at a plot of residual vs. predicted values
plot(trans_resids ~ trans_preds, data = xdata,
     xlab = "Predicted Values",
     ylab = "Residuals")
#Perform a Shapiro-Wilk test for normality of residuals
shapiro.test(trans_resids)  #Shapiro test significant. Transform
#Perform Levene's Test for homogeneity of variances
leveneTest(trans_cort~sex,data=xdata, center = mean)
leveneTest(trans_cort~sex,data=xdata, center = median)
leveneTest(trans_cort~colony,data=xdata, center = mean)
leveneTest(trans_cort~colony,data=xdata, center = median)

library(MASS)
b=boxcox(dheas~age.serum+sex+colony, data=xdata)
lamda=b$x
lik=b$y
bc=cbind(lamda, lik)
bc [order(lik), ]# choose the value o lamda that has the higher likehood
trans_model<-lm(dheas^0.46~sex+age.serum+colony, data=xdata)
par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(trans_model)
trans_cort<-xdata$dheas^0.46

#TESTING ASSUMPTIONS
#Generate residual and predicted values
trans_resids <- residuals(trans_model)
trans_preds <- predict(trans_model)
sq_preds <- trans_preds^2
#Look at a plot of residual vs. predicted values
plot(trans_resids ~ trans_preds, data = xdata,
     xlab = "Predicted Values",
     ylab = "Residuals")
#Perform a Shapiro-Wilk test for normality of residuals
shapiro.test(trans_resids)
#Perform Levene's Test for homogeneity of variances
leveneTest(trans_cort~sex,data=xdata, center = mean)
leveneTest(trans_cort~sex,data=xdata, center = median)
leveneTest(trans_cort~colony,data=xdata, center = mean)
leveneTest(trans_cort~colony,data=xdata, center = median)

library(AICcmodavg)
library(lme4)
m0 <- glm(trans_cort~1, data=xdata)
m1 <- glm(trans_cort~age.serum+sex+colony+age.serum*sex+age.serum*colony+sex*colony, data=xdata,family=gaussian)
AIC(m0,m1)
summary(m1)
m2 <- glm(trans_cort~age.serum+sex+colony+age.serum*colony+sex*colony, data=xdata,family=gaussian)
AIC(m2,m0)
summary(m2)
m3 <- glm(trans_cort~age.serum+sex+colony+age.serum*colony, data=xdata,family=gaussian)
AIC(m3,m0)
summary(m3)
m4 <- glm(trans_cort~age.serum+sex+colony, data=xdata,family=gaussian)
AIC(m4,m0)
summary(m4)
m5 <- glm(trans_cort~age.serum+sex, data=xdata,family=gaussian)
AIC(m5,m0)
summary(m5)
m6 <- lm(trans_cort~age.serum, data=xdata)
AIC(m6,m0)
summary(m6)

#Polynomial
xdata$age.serum2<-xdata$age.serum^2
xdata$age.serum3<-xdata$age.serum^3

m8<-lm(trans_cort~age.serum+age.serum2+age.serum3+sex+colony,data=xdata)
summary(m8)
m9<-lm(trans_cort~age.serum+age.serum2+age.serum3,data=xdata)
summary(m9)


library(MuMIn)
model.sel(m0,m1,m2,m3,m4,m5,m6,m8,m9)

summary(m6)

m6 <- lm(trans_cort~age.serum, data=xdata)
parameters::model_parameters(anova(m6))
library(heplots)
etasq(m6, anova = 2)
library(rsq)
r.squaredGLMM(m6) 

#Figure 1
p2<-ggplot(xdata, aes(x=age.serum, y=trans_cort)) +
  geom_point(size=3) +
  geom_smooth(method=lm, col="black")+
  labs(y= "[Power-transformed] DHEAS levels", x="Age")+ #, title=c("A")
  theme_bw()+
  theme(
    plot.background = element_blank()
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()  ) +
  theme(plot.title = element_text(hjust = 0, size=rel(1.5), face = "bold"))+
  theme(plot.margin=unit(c(0.2,0.2,0.2,0.2),units="cm"))+
  theme(axis.title.y = element_text(size = rel(1.5), angle = 90, vjust=0.1, face="bold",margin = margin(r = 15)))+
  theme(axis.title.x = element_text(size = rel(1.5), angle = 00, vjust=1.5,face="bold",margin = margin(10, 15,0,0)))+
  theme(axis.text.y = element_text(colour="black", size = rel(1.25), angle = 00))+
  theme(axis.text.x = element_text(colour="black", size = rel(1.25), angle = 00))
p2


#Cortisol ~ age, sex and colony

model.test<-lm(cortisol~age.serum+sex+colony, data=xdata)
vif(model.test)

par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(model.test)
trans_cort<-xdata$cortisol
trans_model<-model.test

library(MASS)
b=boxcox(cortisol~age.serum+sex+colony, data=xdata)
lamda=b$x
lik=b$y
bc=cbind(lamda, lik)
bc [order(lik), ]
trans_model<-lm(cortisol^0.46~sex+age.serum+colony, data=xdata)
par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(trans_model)
trans_cort<-xdata$cortisol^0.46

#TESTING ASSUMPTIONS
#Generate residual and predicted values
trans_resids <- residuals(trans_model)
trans_preds <- predict(trans_model)
sq_preds <- trans_preds^2
#Look at a plot of residual vs. predicted values
plot(trans_resids ~ trans_preds, data = xdata,
     xlab = "Predicted Values",
     ylab = "Residuals")
#Perform a Shapiro-Wilk test for normality of residuals
shapiro.test(trans_resids)
#Perform Levene's Test for homogeneity of variances
leveneTest(trans_cort~sex,data=xdata, center = mean) 
leveneTest(trans_cort~sex,data=xdata, center = median)
leveneTest(trans_cort~colony,data=xdata, center = mean) #Significative. Non-homogeneous
leveneTest(trans_cort~colony,data=xdata, center = median) #Significative. Non-homogeneous

#Remove outlier
xdata<-subset(xdata, cortisol < 170)

library(MASS)
b=boxcox(cortisol~age.serum+sex+colony, data=xdata)
lamda=b$x
lik=b$y
bc=cbind(lamda, lik)
bc [order(lik), ]
trans_model<-lm(cortisol^0.5~sex+age.serum+colony, data=xdata)
par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(trans_model)
trans_cort<-xdata$cortisol^0.5

library(AICcmodavg)
library(lme4)
m0 <- glm(trans_cort~1, data=xdata)
m1 <- glm(trans_cort~age.serum+sex+colony+age.serum*sex+age.serum*colony+sex*colony, data=xdata,family=gaussian)
AIC(m0,m1)
summary(m1) #remove age*colony
m2 <- glm(trans_cort~age.serum+sex+colony+age.serum*sex+sex*colony, data=xdata,family=gaussian)
AIC(m2,m0)
summary(m2) #remove age*sex
m3 <- glm(trans_cort~age.serum+sex+colony+sex*colony, data=xdata,family=gaussian)
AIC(m3,m0)
summary(m3) #remove sex*colony
m4 <- glm(trans_cort~age.serum+sex+colony, data=xdata,family=gaussian)
AIC(m4,m0)
summary(m4) #remove colony
m5 <- glm(trans_cort~age.serum+sex, data=xdata,family=gaussian)
AIC(m5,m0)
summary(m5)
m6 <- glm(trans_cort~age.serum, data=xdata,family=gaussian)
AIC(m6,m0)
summary(m6)

library(MuMIn)
AIC(m0,m1,m2,m3,m4,m5,m6)
model.sel(m0,m1,m2,m3,m4,m5,m6)


#PCTB models
library(tidyverse)
xdata <- read_delim(file="chimpanzeeDHEASdata.txt",col_names=T,delim="\t") %>%
  mutate(across(where(is.character), factor))

summary(xdata)
#Remove age discrepancy > 1
xdata<-subset(xdata,  age.difference  < 1.1 &  age.difference > -1.1)
summary(xdata)


###GLMM
#PC1
library(car)
model.test<-lm(pc1~age.pctb+sex+dheas+ratio, data=xdata)
vif(model.test) #cort and dheas should be analyzes separately

data <- xdata[ , c("dheas", "ratio")]
data<-na.omit(data)
cor(data)

par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(model.test)
trans_cort<-xdata$pc1

#TESTING ASSUMPTIONS
#Generate residual and predicted values
trans_resids <- residuals(model.test)
trans_preds <- predict(model.test)
sq_preds <- trans_preds^2
#Look at a plot of residual vs. predicted values
plot(trans_resids ~ trans_preds, data = xdata,
     xlab = "Predicted Values",
     ylab = "Residuals")
#Perform a Shapiro-Wilk test for normality of residuals
shapiro.test(trans_resids)
#Perform Levene's Test for homogeneity of variances
leveneTest(trans_cort~sex,data=xdata, center = mean) #Significative. Non-homogeneous
leveneTest(trans_cort~sex,data=xdata, center = median)

library(AICcmodavg)
library(lme4)
m0 <- lm(trans_cort~1, data=xdata)
m1 <- lm(trans_cort~age.pctb+sex+ratio+dheas+age.pctb*sex+age.pctb*ratio+age.pctb*dheas+sex*ratio+dheas*sex, data=xdata)
summary(m1) #remove sex*ratio
m2 <- lm(trans_cort~age.pctb+sex+ratio+dheas+age.pctb*sex+age.pctb*ratio+age.pctb*dheas+dheas*sex, data=xdata)
summary(m2) #remove age*dheas
m3 <- lm(trans_cort~age.pctb+sex+ratio+dheas+age.pctb*sex+age.pctb*ratio+dheas*sex, data=xdata)
summary(m3) #remove dheas*sex
m4 <- lm(trans_cort~age.pctb+sex+ratio+dheas+age.pctb*sex+age.pctb*ratio, data=xdata)
summary(m4) #remove dheas
m5 <- lm(trans_cort~age.pctb+sex+ratio+age.pctb*sex+age.pctb*ratio, data=xdata)
summary(m5) #remove age*ratio
m6 <- lm(trans_cort~age.pctb+sex+ratio+age.pctb*sex, data=xdata)
summary(m6) #remove sex*age
m7 <- lm(trans_cort~age.pctb+sex+ratio, data=xdata)
summary(m7) #remove sex
m8 <- lm(trans_cort~age.pctb+ratio, data=xdata)
summary(m8) #remove age
m9 <- lm(trans_cort~ratio, data=xdata)
summary(m9)

AIC(m0,m1,m2,m3,m4,m5,m6,m7,m8,m9)
library(MuMIn)
model.sel(m0,m1,m2,m3,m4,m5,m6,m7,m8,m9)

summary(m9)

parameters::model_parameters(anova(m9))
library(heplots)
etasq(m9)

#GGPlot
library(ggplot2)
p1<-ggplot(xdata, aes(x=ratio, y=pc1)) +
  geom_point(size=3) +
  stat_smooth(method = "lm", size = 1) +
  labs(y= "Spatial Relationships (PC1)", x="DHEAS/Cortisol ratio")+ #, title=c("A")
  theme_bw()+
  theme(
    plot.background = element_blank()
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()  ) +
  #theme(legend.title = element_text(size=15, face = "bold"))+ #title appearance
  #theme(legend.text = element_text(colour="black", size = 15))+ #legend text
  #theme(legend.title = element_text(size=15, face = "bold"))+ #title appearance
  #theme(legend.position="bottom")+
  theme(plot.margin=unit(c(0.2,0.2,0.2,0.2),units="cm"))+
  theme(axis.title.y = element_text(size = rel(1.5), angle = 90, vjust=0.1, face="bold",margin = margin(r = 15)))+
  theme(axis.title.x = element_text(size = rel(1.5), angle = 00, vjust=1.5,face="bold",margin = margin(10, 15,0,0)))+
  theme(axis.text.y = element_text(colour="black", size = rel(1.25), angle = 00))+
  theme(axis.text.x = element_text(colour="black", size = rel(1.25), angle = 00))

p1

#################PC2
library(car)
model.test<-lm(pc2~age.pctb+sex+dheas+ratio, data=xdata)
vif(model.test) #cort and dheas should be analyzes separately

data <- xdata[ , c("dheas", "ratio")]
data<-na.omit(data)
cor(data)

par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(model.test)
trans_cort<-xdata$pc2

#TESTING ASSUMPTIONS
#Generate residual and predicted values
trans_resids <- residuals(model.test)
trans_preds <- predict(model.test)
sq_preds <- trans_preds^2
#Look at a plot of residual vs. predicted values
plot(trans_resids ~ trans_preds, data = xdata,
     xlab = "Predicted Values",
     ylab = "Residuals")
#Perform a Shapiro-Wilk test for normality of residuals
shapiro.test(trans_resids)
#Perform Levene's Test for homogeneity of variances
leveneTest(trans_cort~sex,data=xdata, center = mean) #Significative. Non-homogeneous
leveneTest(trans_cort~sex,data=xdata, center = median)

library(AICcmodavg)
library(lme4)
m0 <- lm(trans_cort~1, data=xdata)
m1 <- lm(trans_cort~age.pctb+sex+ratio+dheas+age.pctb*sex+age.pctb*ratio+age.pctb*dheas+sex*ratio+dheas*sex, data=xdata)
summary(m1) #remove ratio*sex
m2 <- lm(trans_cort~age.pctb+sex+ratio+dheas+age.pctb*sex+age.pctb*ratio+age.pctb*dheas+dheas*sex, data=xdata)
summary(m2) #remove age*dheas
m3 <- lm(trans_cort~age.pctb+sex+ratio+dheas+age.pctb*sex+age.pctb*ratio+dheas*sex, data=xdata)
summary(m3) #remove age*ratio
m4 <- lm(trans_cort~age.pctb+sex+ratio+dheas+age.pctb*sex+dheas*sex, data=xdata)
summary(m4) #remove age*sex
m5 <- lm(trans_cort~age.pctb+sex+ratio+dheas+dheas*sex, data=xdata)
summary(m5) #remove sex*dheas
m6 <- lm(trans_cort~age.pctb+sex+ratio+dheas, data=xdata)
summary(m6) #remove sex
m7 <- lm(trans_cort~age.pctb+ratio+dheas, data=xdata)
summary(m7) #remove dheas
m8 <- lm(trans_cort~age.pctb+ratio, data=xdata)
summary(m8) #remove ratio
m9 <- lm(trans_cort~age.pctb, data=xdata)
summary(m9)

library(MuMIn)
model.sel(m0,m1,m2,m3,m4,m5,m6,m7,m8,m9)

summary(m9)
anova(m9)
eta.m9<-eta_squared(anova(m9)) # partial = TRUE by default
data.frame(eta.m9)
sum(eta.m9)

#GGPlot
library(ggplot2)
p1<-ggplot(xdata, aes(x=age.pctb, y=pc2)) +
  geom_point(size=3) +
  stat_smooth(method = "lm", size = 1) +
  labs(y= "Tool use and social communication (PC2)", x="Age (years)")+ #, title=c("A")
  theme_bw()+
  theme(
    plot.background = element_blank()
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()  ) +
  #theme(legend.title = element_text(size=15, face = "bold"))+ #title appearance
  #theme(legend.text = element_text(colour="black", size = 15))+ #legend text
  #theme(legend.title = element_text(size=15, face = "bold"))+ #title appearance
  #theme(legend.position="bottom")+
  theme(plot.margin=unit(c(0.2,0.2,0.2,0.2),units="cm"))+
  theme(axis.title.y = element_text(size = rel(1.5), angle = 90, vjust=0.1, face="bold",margin = margin(r = 15)))+
  theme(axis.title.x = element_text(size = rel(1.5), angle = 00, vjust=1.5,face="bold",margin = margin(10, 15,0,0)))+
  theme(axis.text.y = element_text(colour="black", size = rel(1.25), angle = 00))+
  theme(axis.text.x = element_text(colour="black", size = rel(1.25), angle = 00))

p1

#################PC3
#Assumptions
library(moments)
skewness(xdata$pc3, na.rm = TRUE) #between -0.5 and 0.5 is faily symmetrical. > or < 1 is highly skewed
hist(xdata$pc3)
qqnorm(xdata$pc3)
qqline(xdata$pc3)
shapiro.test(xdata$pc3) #very sensitive at high sample size

library(ggpubr)
bxp <- ggboxplot(xdata, x = "sex", y = "pc3")
bxp

library(rstatix)
xdata %>%
  identify_outliers(pc3)

ggqqplot(xdata, "pc3", facet.by = "sex")

library(car)
model.test<-lm(pc3~age.pctb+sex+dheas+ratio, data=xdata)
vif(model.test) #cort and dheas should be analyzes separately

data <- xdata[ , c("dheas", "ratio")]
data<-na.omit(data)
cor(data)

par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(model.test)
trans_cort<-xdata$pc3

#TESTING ASSUMPTIONS
#Generate residual and predicted values
trans_resids <- residuals(model.test)
trans_preds <- predict(model.test)
sq_preds <- trans_preds^2
#Look at a plot of residual vs. predicted values
plot(trans_resids ~ trans_preds, data = xdata,
     xlab = "Predicted Values",
     ylab = "Residuals")
#Perform a Shapiro-Wilk test for normality of residuals
shapiro.test(trans_resids)
#Perform Levene's Test for homogeneity of variances
leveneTest(trans_cort~sex,data=xdata, center = mean) #Significative. Non-homogeneous
leveneTest(trans_cort~sex,data=xdata, center = median)

library(AICcmodavg)
library(lme4)
m0 <- lm(trans_cort~1, data=xdata)
m1 <- lm(trans_cort~age.pctb+sex+ratio+dheas+age.pctb*sex+age.pctb*ratio+age.pctb*dheas+sex*ratio+dheas*sex, data=xdata)
summary(m1) #remove age*dheas
m2 <- lm(trans_cort~age.pctb+sex+ratio+dheas+age.pctb*sex+age.pctb*ratio+sex*ratio+dheas*sex, data=xdata)
summary(m2) #remove sex*dheas

library(MuMIn)
model.sel(m0,m1,m2)

m3 <- lm(trans_cort~age.pctb+sex+ratio+dheas+age.pctb*sex+age.pctb*ratio+sex*ratio+dheas*sex, data=xdata)
summary(m3) #remove age*ratio
m4 <- lm(trans_cort~age.pctb+sex+ratio+dheas+age.pctb*sex+sex*ratio+dheas*sex, data=xdata)
summary(m4) #remove ratio*sex
m5 <- lm(trans_cort~age.pctb+sex+ratio+dheas+age.pctb*sex+dheas*sex, data=xdata)
summary(m5) #remove sex*dheas
m6 <- lm(trans_cort~age.pctb+sex+ratio+dheas+age.pctb*sex, data=xdata)
summary(m6) #remove sex*age
m7 <- lm(trans_cort~age.pctb+sex+ratio+dheas, data=xdata)
summary(m7) #remove ratio
m8 <- lm(trans_cort~age.pctb+sex+dheas, data=xdata)
summary(m8) #remove age
m9 <- lm(trans_cort~sex+dheas, data=xdata)
summary(m9) #remove dheas
m10 <- lm(trans_cort~sex, data=xdata)
summary(m10)


library(MuMIn)
model.sel(m0,m1,m2,m3,m4,m5,m6,m7,m8,m9,m10)