Sexual signaling and sociosexual behaviours in relation to rank, parasites, hormones, and age in male vervet monkeys (Chlorocebus pygerythrus) in Uganda

Authors: Karin P. Snyder, Dina Greenberg, Taylor Fane, Alessandro Filazzola, Gabriela F. Mastromonaco and Valérie A.M. Schoof

Load Libraries and Data

```{r}
## Libraries
library(vegan)
library(tidyverse)
library(MuMIn)
library(car) 
library(lme4)
library(lmerTest)  
library(performance) 
library(glmmTMB) 
library(fitdistrplus) 
library(logspline)
library(pwr)
library(Hmisc)

## Load data
data <- read.csv("NewData8.csv", stringsAsFactors = FALSE)
data[data == "#VALUE!"] <- NA ## convert Excel errors to NA
data[,15:22] <- apply(data[,15:22], 2, as.numeric)
data$AgeClass <- as.factor(data$AgeClass) ## convert ageclass to factor

## Select response variables
responseVars <- c("BlueLUM","BlueBY","BlueRG","RedLUM","RedBY", "RedRG")

## Select predictor variables
predictorVars<- c( "Age", "Prevalence", "Richness", "GCMEAN", "ARMEAN", "Ordinal")

#Drop AgeClass
data2 <- data[ -c(4) ]

## Drop NA

data6 <- na.omit(data2)

##Load Sociosexual Data

data_socio <- read.csv("Data_counts_9Aug.csv", stringsAsFactors = FALSE)

```

Descriptive Statistics- Outcome Variables

```{r}
# Load necessary libraries
library(dplyr)

# Select relevant columns and group by Session
descriptive_stats <- data6 %>% group_by(Session) %>% summarise(
    BlueLUM_range = range(BlueLUM),
    BlueLUM_mean = mean(BlueLUM, na.rm = TRUE),
    BlueLUM_sd = sd(BlueLUM, na.rm = TRUE),
    
    BlueBY_range = range(BlueBY),
    BlueBY_mean = mean(BlueBY, na.rm = TRUE),
    BlueBY_sd = sd(BlueBY, na.rm = TRUE),
    
    BlueRG_range = range(BlueRG),
    BlueRG_mean = mean(BlueRG, na.rm = TRUE),
    BlueRG_sd = sd(BlueRG, na.rm = TRUE),
    
    RedLUM_range = range(RedLUM),
    RedLUM_mean = mean(RedLUM, na.rm = TRUE),
    RedLUM_sd = sd(RedLUM, na.rm = TRUE),
    
    RedBY_range = range(RedBY),
    RedBY_mean = mean(RedBY, na.rm = TRUE),
    RedBY_sd = sd(RedBY, na.rm = TRUE),
    
    RedRG_range = range(RedRG),
    RedRG_mean = mean(RedRG, na.rm = TRUE),
    RedRG_sd = sd(RedRG, na.rm = TRUE)
  )

# Save the results as a CSV file to the specified path
write.csv(descriptive_stats, "C:/Users/karin/OneDrive/Desktop/FINAL_outcome_descriptive.csv", row.names = FALSE)

Power analyses Model 1

```{r}


predictors <- c("BlueLUM", "BlueBY", "BlueRG", "RedLUM", "RedBY", "RedRG")
outcomes <- c("Ordinal", "Prevalence", "Richness", "ARMEAN", "GCMEAN")

power_results <- data.frame(Outcome = character(),
                            Predictor = character(),
                            R2 = numeric(),
                            f2 = numeric(),
                            v = numeric(),  # Degrees of freedom for error term
                            u = numeric(),  # Number of predictors
                            Power = numeric(),
                            stringsAsFactors = FALSE)


alpha <- 0.05  

# Loop over each combination of predictor and outcome
for (outcome in outcomes) {
  for (predictor in predictors) {
    # Fit the linear model for the current combination
    model <- lm(as.formula(paste(outcome, "~", predictor)), data = data6)
    
R2 <- summary(model)$r.squared
    
f2 <- R2 / (1 - R2)
    
df_error <- summary(model)$df[2]
    
numPredictors <- 1  

power_result <- pwr.f2.test(u = numPredictors, v = df_error, f2 = f2, sig.level = alpha)
    
power_results <- rbind(power_results, data.frame(Outcome = outcome, Predictor = predictor, R2 = R2, f2 = f2, v = df_error, u = numPredictors,Power = power_result$power))}}

desktop_path <- "C:/Users/karin/OneDrive/Desktop/FINAL_power_analysis_results.csv"
desktop_path <- gsub("YourUsername", Sys.getenv("USERNAME"), desktop_path)
write.csv(power_results, file = desktop_path, row.names = FALSE)

Power analyses Model 2

```{r}


#Rename columns
data_socio <- data_socio %>%
  rename(BlueBY = BlueB.Y, BlueRG = BlueR.G, RedBY = RedB.Y, RedRG = RedR.G)

# Define outcomes and predictors
outcomes <- c("BlueLUM", "BlueBY", "BlueRG", "RedLUM", "RedBY", "RedRG")
predictors <- c("COP", "PR", "MR","Ordinal", "Session")

power_results <- data.frame(Outcome = character(),
                            Predictor = character(),
                            R2 = numeric(),
                            f2 = numeric(),
                            v = numeric(),  # Degrees of freedom for error term
                            u = numeric(),  # Number of predictors
                            Power = numeric(),
                            stringsAsFactors = FALSE)

alpha <- 0.05  # Significance level

for (outcome in outcomes) {
  for (predictor in predictors) {
    # Fit the linear model for the current combination
    model <- lm(as.formula(paste(outcome, "~", predictor)), data = data_socio)
    

R2 <- summary(model)$r.squared

f2 <- R2 / (1 - R2)
    
df_error <- summary(model)$df[2]
    
numPredictors <- 1  # Only one predictor in each model
    
  
power_result <- pwr.f2.test(u = numPredictors, v = df_error, f2 = f2, sig.level = alpha)

power_results <- rbind(power_results, data.frame(Outcome = outcome,
                                                     Predictor = predictor,
                                                     R2 = R2,
                                                     f2 = f2,
                                                     v = df_error,
                                                     u = numPredictors,
                                                     Power = power_result$power))}}


desktop_path <- "C:/Users/karin/OneDrive/Desktop/FINAL2_power_analysis_results_sociosexual.csv"
write.csv(power_results, file = desktop_path, row.names = FALSE)

Predictor/Outcome correlations

```{r}

predictors <- c("BlueLUM", "BlueBY", "BlueRG", "RedLUM", "RedBY", "RedRG")
outcomes <- c("Ordinal", "Prevalence", "Richness", "ARMEAN", "GCMEAN")

results <- data.frame(Predictor = character(),
                      Outcome = character(),
                      Correlation = numeric(),
                      Uncorrected_p = numeric(),
                      BH_p = numeric(),
                      stringsAsFactors = FALSE)

for (pred in predictors) {
  for (out in outcomes) {
    # Calculate Spearman correlation and p-value
    test_result <- cor.test(data6[[pred]], data6[[out]], method = "spearman", use = "complete.obs")
    
  
results <- rbind(results, data.frame(Predictor = pred,
                                         Outcome = out,
                                         Correlation = test_result$estimate,
                                         Uncorrected_p = test_result$p.value,
                                         BH_p = NA)) # Initialize BH_p as NA
  }
}

# Additional correlations to calculate
additional_pairs <- list(
  c("ARMEAN", "Prevalence"),
  c("ARMEAN", "Richness"),
  c("GCMEAN", "Prevalence"),
  c("GCMEAN", "Richness"),
  c("GCMEAN", "ARMEAN")
)

# Loop over additional pairs
for (pair in additional_pairs) {
  pred <- pair[1]
  out <- pair[2]
  
  # Calculate Spearman correlation and p-value
  test_result <- cor.test(data6[[pred]], data6[[out]], method = "spearman", use = "complete.obs")
  
  # Add the results to the results data frame
  results <- rbind(results, data.frame(Predictor = pred,
                                       Outcome = out,
                                       Correlation = test_result$estimate,
                                       Uncorrected_p = test_result$p.value,
                                       BH_p = NA)) # Initialize BH_p as NA
}

# Apply Benjamini-Hochberg correction to p-values
results$BH_p <- p.adjust(results$Uncorrected_p, method = "BH")


desktop_path <- "C:/Users/karin/OneDrive/Desktop/FINAL_pred_outcome_corrs.csv"
write.csv(results, file = desktop_path, row.names = FALSE)

Intermale Analyses

Scrotal Luminance

```{r}
#model1A_group_maleID <- lmer(BlueLUM ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season + (1|Group) + (1|MaleID), data=data6, na.action = "na.fail")
#model convergence issues 

#model1A_group <- lmer(BlueLUM ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season + (1|Group), data=data6, na.action = "na.fail")
#model convergence issues 

#model1A_maleID <- lmer(BlueLUM ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season + (1|MaleID), data=data6, na.action = "na.fail")
#check_collinearity(model1A_maleID) # low VIF <4
#options(na.action = na.fail)
#dd1A_maleID <- dredge(model1A_maleID) # Warning: comparing models fitted by REML
#options(na.action = na.omit)
## model convergence issues when dredging, cannot use maleID as random variable

#Change column headings

data6 <- data6 %>%
  rename(Diversity = Richness,PropInfected = Prevalence, Season = Session)


## Final model with no random variables: 
model1A <- lm(BlueLUM ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season, data=data6, na.action = "na.fail")
check_collinearity(model1A) # low VIF <4

#dredge model selection
options(na.action = na.fail) 
dredge1A <- dredge(model1A)
options(na.action = na.omit)

#model averaging
dredge1A_modelaverage <- model.avg(dredge1A, subset = delta <= 7)
summary(dredge1A_modelaverage) 
#Ordinal Rank is the only significant variable 

#Getting 95% confidence intervals 
confint(dredge1A_modelaverage, full=TRUE) 
#confirms ordinal rank is the only significant variable 
```

Scrotal B/Y Channel

```{r}
#model1B_group_maleID <- lmer(BlueBY ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season + (1|Group) + (1|MaleID), data=data6, na.action = "na.fail")
#check_collinearity(model1B_group_maleID) # low VIF < 3
#options(na.action = na.fail)
#dd1B_group_maleID <- dredge(model1B_group_maleID) ## boundary (singular) fit: see help('isSingular')
#options(na.action = na.omit)
## model convergence issues when dredging, cannot use group and maleID as random variable 

#model1B_group <- lmer(BlueBY ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season + (1|Group), data=data6, na.action = "na.fail")
## model convergence issues

#model1B_maleID <- lmer(BlueBY ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season + (1|MaleID), data=data6, na.action = "na.fail")
#check_collinearity(model1B_maleID) # low VIF < 4
#options(na.action = na.fail)
#dd1B_maleID <- dredge(model1B_maleID) ## boundary (singular) fit: see help('isSingular')
#options(na.action = na.omit)
## model convergence issues when dredging, cannot use maleID as random variable


## Final model with no random variables: 
model1B <- lm(BlueBY ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season, data=data6, na.action = "na.fail")
check_collinearity(model1B) # low VIF < 4

#dredge model selection
options(na.action = na.fail) 
dredge1B <- dredge(model1B)
options(na.action = na.omit)

#model averaging
dredge1B_modelaverage <- model.avg(dredge1B, subset = delta <= 7)
summary(dredge1B_modelaverage) 
#No significant variables 

#Getting 95% confidence intervals 
confint(dredge1B_modelaverage, full=TRUE) 
#confirms NULL MODEL is the best model 

```

Scrotal R/G Channel

```{r}
#model1C_group_maleID <- lmer(BlueRG ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season + (1|Group) + (1|MaleID), data=data6, na.action = "na.fail")
#check_collinearity(model1C_group_maleID) # low VIF >3
#options(na.action = na.fail)
#dd1C_group_maleID <- dredge(model1C_group_maleID) ## boundary (singular) fit: see help('isSingular')
#options(na.action = na.omit)
## model convergence issues when dredging, cannot use group and maleID as random variable 

#model1C_group <- lmer(BlueRG ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season + (1|Group), data=data6, na.action = "na.fail")
## model convergence issues 

#model1C_maleID <- lmer(BlueRG ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season + (1|MaleID), data=data6, na.action = "na.fail")
#check_collinearity(model1C_maleID) # low VIF >3
#options(na.action = na.fail)
#dd1C_maleID <- dredge(model1C_maleID) ## boundary (singular) fit: see help('isSingular')
#options(na.action = na.omit)
## model convergence issues when dredging, cannot use maleID as random variable

## Final model with no random variables: 
model1C <- lm(BlueRG ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season, data=data6, na.action = "na.fail")
check_collinearity(model1C) # low VIF < 4

#dredge model selection
options(na.action = na.fail) 
dredge1C <- dredge(model1C)
options(na.action = na.omit)

#model averaging
dredge1C_modelaverage <- model.avg(dredge1C, subset = delta <= 7)
summary(dredge1C_modelaverage) 
#No significant variables 

#Getting 95% confidence intervals 
confint(dredge1C_modelaverage, full=TRUE) 
#confirms NULL MODEL is the best model 
```

Penile Luminance

```{r}
#model1D_group_maleID <- lmer(RedLUM ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season + (1|Group) + (1|MaleID), data=data6, na.action = "na.fail")
## model convergence issues

#model1D_group <- lmer(RedLUM ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season + (1|Group), data=data6, na.action = "na.fail")
## model convergence issues

#model1D_maleID <- lmer(RedLUM ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season + (1|MaleID), data=data6, na.action = "na.fail")
#check_collinearity(model1D_maleID) # low VIF < 4
#options(na.action = na.fail)
#dd1D_maleID <- dredge(model1D_maleID) ## Warning: comparing models fitted by REML
#options(na.action = na.omit)
## model convergence issues when dredging, cannot use maleID as random variable

## Final model with no random variables: 
model1D <- lm(RedLUM ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season, data=data6, na.action = "na.fail")
check_collinearity(model1D) # low VIF < 4

#dredge model selection
options(na.action = na.fail) 
dredge1D <- dredge(model1D)
options(na.action = na.omit)

#model averaging
dredge1D_modelaverage <- model.avg(dredge1D, subset = delta <= 7)
summary(dredge1D_modelaverage) 
#No significant variables 

#Getting 95% confidence intervals 
confint(dredge1D_modelaverage, full=TRUE) 
#confirms NULL MODEL is the best model 
```

Penile B/Y Channel

```{r}
#model1E_group_maleID <- lmer(RedBY ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season + (1|Group) + (1|MaleID), data=data6, na.action = "na.fail")
## model convergence issues 

#model1E_group <- lmer(RedBY ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season + (1|Group), data=data6, na.action = "na.fail")
## model convergence issues 

#model1E_maleID <- lmer(RedBY ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season + (1|MaleID), data=data6, na.action = "na.fail")
## model convergence issues 

## Final model with no random variables: 
model1E <- lm(RedBY ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season, data=data6, na.action = "na.fail")
check_collinearity(model1E) # low VIF < 4

#dredge model selection
options(na.action = na.fail) 
dredge1E <- dredge(model1E)
options(na.action = na.omit)

#model averaging
dredge1E_modelaverage <- model.avg(dredge1E, subset = delta <= 7)
summary(dredge1E_modelaverage) 
#NULL model is best model 

#Getting 95% confidence intervals 
confint(dredge1E_modelaverage, full=TRUE) 
#confirms NULL model 

```

Penile R/G Channel

```{r}
#model1F_group_maleID <- lmer(RedRG ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season + (1|Group) + (1|MaleID), data=data6, na.action 
## model convergence issues 

#model1F_group <- lmer(RedRG ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season + (1|Group), data=data6, na.action = "na.fail")
#check_collinearity(model1F_group) ## low VIF > 4
#options(na.action = na.fail)
#dd1F_group <- dredge(model1F_group) ## boundary (singular) fit: see help('isSingular')
##options(na.action = na.omit)
## model convergence issues when dredging, cannot use group as random variable 

#model1F_maleID <- lmer(RedRG ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season + (1|MaleID), data=data6, na.action = "na.fail")
## model convergence issues 

## Final model with no random variables: 
model1F <- lm(RedRG ~ Ordinal + Diversity + GCMEAN + log(ARMEAN) + PropInfected + Season, data=data6, na.action = "na.fail")
check_collinearity(model1F) # low VIF > 4

#dredge model selection
options(na.action = na.fail) 
dredge1F <- dredge(model1F)
options(na.action = na.omit)

#model averaging
dredge1F_modelaverage <- model.avg(dredge1F, subset = delta <= 7)
summary(dredge1F_modelaverage) 
#Ordinal rank is significant 

#Getting 95% confidence intervals 
confint(dredge1F_modelaverage, full=TRUE) 
#confirms ordinal rank is significant 

```

Sociosexual Analysis (Model 2)

Copulations

```{r}
## MaleID and group as random variables 
#model2A_group_maleID <- glmer(COPcount ~ offset(log(FollowTime)) + BlueR.G + BlueLUM + RedLUM + BlueB.Y + RedB.Y + RedR.G + Ordinal + Season + (1|Group) + (1|MaleID), family = poisson, data = data7, na.action = "na.fail") 
#model convergence issues 

## Male ID as random variable
#model2A_maleID <- glmer(COPcount ~ offset(log(FollowTime)) + BlueR.G + BlueLUM + RedLUM + BlueB.Y + RedB.Y + RedR.G + Ordinal + Season + (1|MaleID), family = poisson, data = data7, na.action = "na.fail") 
#check_collinearity(model2A_maleID) # low VIF (<4)
#check_overdispersion(model2A_maleID) 
#OVER DISPERSION DETECTED (dispersion ratio: 1.649)

# re-run as negative binomial to fix over-dispersion 
#model2A_maleID_NB <- glmer.nb(COPcount ~ offset(log(FollowTime)) + BlueR.G + BlueLUM + RedLUM + BlueB.Y + RedB.Y + RedR.G + Ordinal + Season + (1|MaleID), data = data7, na.action = "na.fail") 
#model convergence issues 

## group as random variable 
#model2A_group <- glmer(COPcount ~ offset(log(FollowTime)) + BlueR.G + BlueLUM + RedLUM + BlueB.Y + RedB.Y + RedR.G + Ordinal + Season + (1|Group), family = poisson, data = data7, na.action = "na.fail") 
#check_collinearity(model2A_group) # low VIF (<4)
#check_overdispersion(model2A_group) 
#OVER DISPERSION DETECTED (dispersion ratio: 1.785)

# re-run as negative binomial to fix over-dispersion 
#model2A_group_NB <- glmer.nb(COPcount ~ offset(log(FollowTime)) + BlueR.G + BlueLUM + RedLUM + BlueB.Y + RedB.Y + RedR.G + Ordinal + Season + (1|Group), data = data7, na.action = "na.fail") 
#model convergence issues (will try without randoms)

## No random variables 
#model2A <- glm(COPcount ~ offset(log(FollowTime)) + BlueR.G + BlueLUM + RedLUM + BlueB.Y + RedB.Y + RedR.G + Ordinal + Season, family = poisson, data = data7, na.action = "na.fail") 
#check_collinearity(model2A) # low VIF (<4)
#check_overdispersion(model2A) 
#OVER DISPERSION DETECTED (dispersion ratio: 1.717)

#Rename dataframe
data7 <- data_socio

#Change column headings
data7 <- data7 %>%
  rename(Diversity = Richness,PropInfected = Prevalence, Season = Session)

# re-run as negative binomial to fix over-dispersion 
model2A_NB <- glm.nb(COPcount ~ offset(log(FollowTime)) + BlueR.G + BlueLUM + RedLUM + BlueB.Y + RedB.Y + RedR.G + Ordinal + Season, data = data7, na.action = "na.fail", control=glm.control(maxit=50)) 
summary(model2A_NB)
check_collinearity(model2A_NB) # low VIF (<4)
check_overdispersion(model2A_NB) # no overdispersion 

#dredge model selection 
options(na.action = na.fail)
dredge_model2A_NB <- dredge(model2A_NB)
options(na.action = "na.omit")

#model averaging 
model2A_NB_modelaverage <- model.avg(dredge_model2A_NB, subset = delta <= 7)
summary(model2A_NB_modelaverage) #No significant variables 

#Getting 95% confidence intervals
confint(model2A_NB_modelaverage, full=TRUE) #confirms NULL model is best model

```

Mate Presentations

```{r}
## MaleID and group as random variables 
#model2B_group_maleID <- glmer(PRcount ~ offset(log(FollowTime)) + BlueR.G + BlueLUM + RedLUM + BlueB.Y + RedB.Y + RedR.G + Ordinal + Season + (1|Group) + (1|MaleID), family = poisson, data = data7, na.action = "na.fail") 
#model convergence issues 

## group as random variable 
#model2B_group <- glmer(PRcount ~ offset(log(FollowTime)) + BlueR.G + BlueLUM + RedLUM + BlueB.Y + RedB.Y + RedR.G + Ordinal + Season + (1|Group), family = poisson, data = data7, na.action = "na.fail") 
#model convergence issues 

## Male ID as random variable
#model2B_maleID <- glmer(PRcount ~ offset(log(FollowTime)) + BlueR.G + BlueLUM + RedLUM + BlueB.Y + RedB.Y + RedR.G + Ordinal + Season + (1|MaleID), family = poisson, data = data7, na.action = "na.fail") 
#model convergence issues 

## No random variables
model2B <- glm(PRcount ~ offset(log(FollowTime)) + BlueR.G + BlueLUM + RedLUM + BlueB.Y + RedB.Y + RedR.G + Ordinal + Season, family = poisson, data = data7, na.action = "na.fail") 
check_collinearity(model2B) # low VIF (<6)
check_overdispersion(model2B) # no over dispersion 

#dredge model selection
options(na.action = na.fail)
dredge_model2B <- dredge(model2B)
options(na.action = "na.omit")

#model averaging No random variables
model2B_modelaverage <- model.avg(dredge_model2B, subset = delta <= 7)
summary(model2B_modelaverage) #RedLUM and Season are significant  

#Getting 95% confidence intervals No random variables
confint(model2B_modelaverage, full=TRUE) #confirms RedLUM and Season are significant  

```

Mate Refusals

```{r}
## MaleID and group as random variables
#model2C_group_maleID <- glmer(MRcount ~ offset(log(FollowTime)) + BlueR.G + BlueLUM + RedLUM + BlueB.Y + RedB.Y + RedR.G + Ordinal + Season + (1|Group) + (1|MaleID), family = poisson, data = data7, na.action = "na.fail") 
#model convergence issues 

## group as random variable 
#model2C_group <- glmer(MRcount ~ offset(log(FollowTime)) + BlueR.G + BlueLUM + RedLUM + BlueB.Y + RedB.Y + RedR.G + Ordinal + Season + (1|Group), family = poisson, data = data7, na.action = "na.fail") 
#model convergence issues 

## Male ID as random variable 
#model2C_maleID <- glmer(MRcount ~ offset(log(FollowTime)) + BlueR.G + BlueLUM + RedLUM + BlueB.Y + RedB.Y + RedR.G + Ordinal + Season + (1|MaleID), family = poisson, data = data7, na.action = "na.fail") 
#model convergence issues 

## No random variables 
#model2C <- glm(MRcount ~ offset(log(FollowTime)) + BlueR.G + BlueLUM + RedLUM + BlueB.Y + RedB.Y + RedR.G + Ordinal + Season, family = poisson, data = data7, na.action = "na.fail") 
#check_collinearity(model2C) ## VIF > 8 
#check_overdispersion(model2C) 
# OVER DISPERSION DETECTED (RATIO = 2.025) 

# re-run as negative binomial to fix over-dispersion 
#model2C_NB <- glmer.nb(MRcount ~ offset(log(FollowTime)) + BlueR.G + BlueLUM + RedLUM + BlueB.Y + RedB.Y + RedR.G + Ordinal + Season + (1|Group), data = data7, control=glmerControl(optCtrl=list(maxfun=30000)), na.action= "na.fail") # Warning: NaNs produced
#check_collinearity(model2C_NB) ## VIF < 8 
#check_overdispersion(model2C_NB) #OVERDISPERSION DETECTED (RATIO = 2.025)

## re-run as quasi-poisson 
#model2C_QP <- glm(MRcount ~ offset(log(FollowTime)) + BlueR.G + BlueLUM + RedLUM + BlueB.Y + RedB.Y + RedR.G + Ordinal + Season, family = quasipoisson, data = data7, na.action = "na.fail") 
#check_collinearity(model2C_QP) ## VIF > 8 
#check_overdispersion(model2C_QP) #OVERDISPERSION DETECTED (RATIO = 2.025)

## zero inflated poisson
#model2C_group_maleID <-glmmTMB(MRcount ~ offset(log(FollowTime)) + BlueR.G + BlueLUM + RedLUM + BlueB.Y + RedB.Y + RedR.G + Ordinal + Season, data = data7, family = poisson() , ziformula=~1)
#check_collinearity(model2C_group_maleID) # VIF < 8
#check_zeroinflation(model2C_group_maleID) # model ok (no zero inflation)
#check_overdispersion(model2C_group_maleID) #OVERDISPERSION DETECTED (RATIO = 2.121)

## zero-inflated negative binomial 2 
#model2C_group_maleID_NB <- update(model2C_group_maleID,family=nbinom2)
#model failed to converge

## zero-inflated negative binomial 1 
#model2C_group_maleID_NB <- update(model2C_group_maleID,family=nbinom1)
#model failed to converge

## Mating refusals = extreme zero inflation, cannot analyze the data

```

Intramale Analyses

Short Term with Benjamini-Hochberg correction

```{r}

STDATA <- read.csv("Intramale_ST.csv")
# Select the relevant columns
variables <- c('BLUMST', 'BBYST', 'BRGST', 'RLUMST', 'RBYST', 'RRGST', 'ORDST', 'PREVST', 'RCHST', 'GCST', 'ARST')
STDATA_selected <- STDATA[, variables]

# Calculate the Spearman correlation matrix
correlation_matrix <- rcorr(as.matrix(STDATA_selected), type = "spearman")

# Extract the p-values
p_values <- correlation_matrix$P

# Apply Benjamini-Hochberg correction
p_values_corrected_bh <- p.adjust(p_values, method = "BH")

# Reshape the corrected p-values back into a matrix
p_values_corrected_matrix_bh <- matrix(p_values_corrected_bh, nrow = length(variables), ncol = length(variables))
dimnames(p_values_corrected_matrix_bh) <- list(variables, variables)

print("Spearman Correlation Matrix:")
print(correlation_matrix$r)
print("Corrected P-Values Matrix (BH):")
print(p_values_corrected_matrix_bh)

##CSV to Desktop:

combined_matrix <- cbind(as.data.frame(correlation_matrix$r), as.data.frame(p_values_corrected_matrix_bh))
colnames(combined_matrix) <- paste0(rep(variables, 2), "_", rep(c("Corr", "PValue"), each = length(variables)))

file_path <- "C:/Users/karin/OneDrive/Desktop/FINAL_RII_Corrs_ST.csv"

#write.csv(combined_matrix, file = file_path, row.names = TRUE)

```

Long term with Benjamini-Hochberg correction

```{r}

LTDATA <- read.csv("Intramale_LT(1).csv")
# List of expected columns for LTDATA
variables_lt <- c('BLUMLT', 'BBYLT', 'BRGLT', 'RLUMLT', 'RBYRT', 'RRGLT', 'ORDLT', 'GCLT', 'ARLT')
LTDATA_selected <- LTDATA[, variables_lt]

# Calculate the Spearman correlation matrix
correlation_matrix_lt <- rcorr(as.matrix(LTDATA_selected), type = "spearman")

# Extract the p-values
p_values_lt <- correlation_matrix_lt$P

# Apply Benjamini-Hochberg correction
p_values_corrected_bh_lt <- p.adjust(p_values_lt, method = "BH")

# Reshape the corrected p-values back into a matrix
p_values_corrected_matrix_bh_lt <- matrix(p_values_corrected_bh_lt, nrow = length(variables_lt), ncol = length(variables_lt))
dimnames(p_values_corrected_matrix_bh_lt) <- list(variables_lt, variables_lt)

# Print the correlation matrix and the corrected p-values matrix
print("Spearman Correlation Matrix:")
print(correlation_matrix_lt$r, digits=2)
print("Corrected P-Values Matrix (BH):")
print(p_values_corrected_matrix_bh_lt)

##CSV to Desktop:

combined_matrix_lt <- cbind(as.data.frame(correlation_matrix_lt$r), as.data.frame(p_values_corrected_matrix_bh_lt))
colnames(combined_matrix_lt) <- paste0(rep(variables_lt, 2), "_", rep(c("Corr", "PValue"), each = length(variables_lt)))
file_path_lt <- "C:/Users/karin/OneDrive/Desktop/RII_Corrs_LT.csv"
#write.csv(combined_matrix_lt, file = file_path_lt, row.names = TRUE)

```