# 1- EXAMPLE CODE FOR COMPUTING COLOUR SATURATION USING AN AVIAN VISUAL MODEL.

# Before running this, load the package "pavo" (used version 0.5-2), and use any standard 
# method (e.g., the function "getspec") to import reflectance spectra (in our case, from 300 
# to 700 nm, with reflectance <1% flattened to 1%) into a rspec object named "refectance_data".


visual_model <- vismodel(refectance_data, visual = "avg.uv", illum='D65', relative=TRUE)
# This visual model is in relative mode to computed colour metrics from a visual tetracolour space 
# (as opposed to absolute mode to calculate colour distances in Just Noticeaable Differences)

# Asks for pavo's standard set of color metrics in the tetrachomatic space, and then exports them 
# to a text file. Colour saturation is labelled "r.vec"
tetracolour_space_metrics <- tcs(visual_model)
write.table(tetracolour_space_metrics, file="tetracolour_space_metrics.txt", quote=F, sep="\t")


 

# -----------------------------------------------------------------------------------------------


# 2- EXAMPLE CODE FOR QUANTILE REGRESSIONS.

# Before running this, load the package "quantreg" (used version 5.35), and import the data
# in the Supplementary Material file "1_Mean_values_per_species", excluding red, pied or 
# black-rumped species (marked **),into an object named "spp_data". 
# This code exemplifies quantile regression of female on male breast colour saturation.


# deleting lines with missing data
data_breast_sat <- spp_data[complete.cases(spp_data[,c("MB_sat","FB_sat")]),]
# 90% quantile regression
fit_90FMb_sat <- rq ( FB_sat ~ MB_sat, tau = .9, data = data_breast_sat )
# 10% quantile regression
fit_10FMb_sat <- rq ( FB_sat ~ MB_sat, tau = .1, data = data_breast_sat )

# Show on screen stats from these models; t and P values use quantreg's default Hall-Sheather 
# bandwidth rule. Stats positions @ coefficients(stats_...) are: [1,] intercept, [2,] slope, 
# [,1] value, [,2] st. error, [,3] t, [,4] p. Also shows 95% confidence intervals from these models.
stats_90FMb_sat <- summary(fit_90FMb_sat, se = "nid")
CI_90FMb_sat <- summary(fit_90FMb_sat, alpha = .05)
stats_90FMb_sat
CI_90FMb_sat
stats_10FMb_sat <- summary(fit_10FMb_sat, se = "nid")
CI_10FMb_sat <- summary(fit_10FM_satb, alpha = .05)
stats_10FMb_sat
CI_10FMb_sat
# Exports these stats to text files
write.table(stats_90FMb_sat[3], file="stats_90FMb_sat.txt", sep="\t")
write.table(CI_90FMb_sat[3], file="CI_90FMb_sat.txt", sep="\t")
write.table(stats_10FMb_sat[3], file="stats_10FMb_sat.txt", sep="\t")
write.table(CI_10FMb_sat[3], file="CI_10FMb_sat.txt", sep="\t")




# -----------------------------------------------------------------------------------------------


# 3- EXAMPLE CODE FOR PGLS REGRESSION

# Before running this, load the package "caper" (used version 1.0.1), and use any standard method 
# (e.g., function "read.nexus") to import a set of phylogenetic trees (in our case, 1000 chronograms 
# taken from the Supplementary Material of Cardoso et al. 2020 Evolution doi:10.1111/evo.13987) into 
# a MultiPhylo object named "trees".
# Also, as in the previous example code (2), import the data in the Supplementary Material file 
# "1_Mean_values_per_species" into an object named "spp_data". In our analyses, we did not included 
# red, pied or black-rumped species (marked **) and species lacking molecular information (marked **). 
# This code exemplifies PGLS regression of female on male breast colour saturation.


# Subseting data
data_breast_sat<-spp_data[,c("Species_name","FB_sat","MB_sat")]
data_breast_sat<-data_breast_sat[complete.cases(data_breast_sat), ]

# compute SD for each variable, to later standardize PGLS beta values
sd.FB_sat<-sd(data_breast_sat$FB_sat,na.rm=T)
sd.MB_sat<-sd(data_breast_sat$MB_sat,na.rm=T)

# Empty table with headers, to later collate results
table_results<-matrix(nrow=1,ncol=9,dimnames=list(c("FB_sat_on_MB_sat"),c("beta_unstandardized","beta_standardized","SE_unstandardized","SE_standardized","new_p-value","CI-lower","CI-upper","lambda","n_spp")))

# Run one PGLS regression model for each phylogenetic tree
FB_sat_on_MB_sat=c() # empty object, to receive the results of each PGLS model in the following loop
for(i in 1:length(trees)){
  tryCatch({
    temp_tree<-trees[[i]]
    class(temp_tree)<-"phylo"
    tree.data=comparative.data(phy=temp_tree,data=data_breast_sat,names.col=Species_name,vcv=T,na.omit=T,warn.dropped=T)
    mod=pgls(FB_sat~MB_sat,data=tree.data,lambda='ML')
    FB_sat_on_MB_sat=rbind(FB_sat_on_MB_sat, data.frame(row.names = NULL, tree = i, logLik = mod$model$log.lik,AIC=mod$aicc,lambda=mod$param[2],int=mod$model$coef[1],int_err=mod$sterr[1],beta=mod$model$coef[2],beta_err=mod$sterr[2]))
  },error=function(e){cat("ERROR:",conditionMessage(e),"\n")})
}
# If the phylogeny has more species than those in the dataset (which is fine), there will be warnings like this:
# Data dropped in compiling comparative data object
# To check which species were dropped, we can type mod$data$dropped$tips

# Compute AIC delta and w for each model
FB_sat_on_MB_sat$delta<-FB_sat_on_MB_sat$AIC-min(FB_sat_on_MB_sat$AIC)
FB_sat_on_MB_sat$w<-exp(-0.5*FB_sat_on_MB_sat$delta)/sum(exp(-0.5*FB_sat_on_MB_sat$delta))

# Exports a text file with results of each individual PGLS model
write.table(FB_sat_on_MB_sat,file = "FB_sat on MB_sat models.txt")

# AIC-weighted model averaging, standardizing PGLS beta values and standard errors (SE)
FB_sat_on_MB_sat.beta<-weighted.mean(FB_sat_on_MB_sat$beta,FB_sat_on_MB_sat$w) # AIC-weighted average beta
FB_sat_on_MB_sat.beta_stand<-FB_sat_on_MB_sat.beta*sd.MB_sat/sd.FB_sat # standardized
w_err=c()
for(i in 1:nrow(FB_sat_on_MB_sat)) {
  w_err[i]<-sqrt(FB_sat_on_MB_sat[i, ]$beta_err^2+(FB_sat_on_MB_sat[i, ]$beta-weighted.mean(FB_sat_on_MB_sat$beta,FB_sat_on_MB_sat$w))^2)
}
FB_sat_on_MB_sat.SE<-weighted.mean(w_err,FB_sat_on_MB_sat$w) # AIC-weighted average SE
FB_sat_on_MB_sat.SE_stand<-FB_sat_on_MB_sat.SE*sd.MB_sat/sd.FB_sat # standardized

# Computing P value and 95% confidence interval of beta
FB_sat_on_MB_sat.t<-abs(FB_sat_on_MB_sat.beta/FB_sat_on_MB_sat.SE) # Value of t
FB_sat_on_MB_sat.tdist<-pt(FB_sat_on_MB_sat.t,34) # This formula used 34 degrees of freedom (36 spp - 1 - 1)
FB_sat_on_MB_sat.p<-2-2*FB_sat_on_MB_sat.tdist # P value
FB_sat_on_MB_sat.lowci<-weighted.mean(FB_sat_on_MB_sat$beta, FB_sat_on_MB_sat$w) - 1.96 * weighted.mean(w_err, FB_sat_on_MB_sat$w) # upper limit of CI
FB_sat_on_MB_sat.lowci<-FB_sat_on_MB_sat.lowci*sd.MB_sat/sd.FB_sat # standardized
FB_sat_on_MB_sat.uppci<-weighted.mean(FB_sat_on_MB_sat$beta, FB_sat_on_MB_sat$w) + 1.96 * weighted.mean(w_err, FB_sat_on_MB_sat$w) # lower limit of CI
FB_sat_on_MB_sat.uppci<-FB_sat_on_MB_sat.uppci*sd.MB_sat/sd.FB_sat # standardized

# AIC-weighted model lambda
FB_sat_on_MB_sat.lambda<-weighted.mean(FB_sat_on_MB_sat$lambda,FB_sat_on_MB_sat$w) # average lambda

# Collate values in the results table, and exports table
table_results[1,]<-c(FB_sat_on_MB_sat.beta,FB_sat_on_MB_sat.beta_stand,FB_sat_on_MB_sat.SE,FB_sat_on_MB_sat.SE_stand,FB_sat_on_MB_sat.p,FB_sat_on_MB_sat.lowci,FB_sat_on_MB_sat.uppci,FB_sat_on_MB_sat.lambda,mod$n)
write.table(format(table_results, scientific=FALSE), "AIC-weighted result, FB_sat on MB_sat.txt")




# -----------------------------------------------------------------------------------------------


# 4- EXAMPLE CODE FOR LAMBDA OF SINGLE TRAITS

# Same notes on R package and data as in the previous example code (3).
# This code exemplifies estimating lambda for male breast colour saturation.


# subseting data
data_lambda_MB_sat<-spp_data[,c("Species_name","MB_sat")]
data_lambda_MB_sat<-data_lambda_MB_sat[complete.cases(data_lambda_MB_sat), ]

# Run one single-trait PGLS for each phylogenetic tree
lambda_MB_sat=c() # empty object, to receive lambda values in the following loop
for(i in 1:length(trees)){
  tryCatch({
    temp_tree<-trees[[i]]
    class(temp_tree)<-"phylo"
    tree.data=comparative.data(phy=temp_tree,data=data_lambda_MB_sat,names.col=Species_name,vcv=T,na.omit=T,warn.dropped=T)
    mod=pgls(MB_sat~1,data=tree.data,lambda='ML')
    lambda_MB_sat=rbind(lambda_MB_sat, data.frame(row.names = NULL, tree = i, logLik = mod$model$log.lik,AIC=mod$aicc,lambda=mod$param[2],int=mod$model$coef[1],int_err=mod$sterr[1],beta=mod$model$coef[2],beta_err=mod$sterr[2]))
  },error=function(e){cat("ERROR:",conditionMessage(e),"\n")})
}

# Compute AIC delta and w for each model
lambda_MB_sat$delta<-lambda_MB_sat$AIC-min(lambda_MB_sat$AIC)
lambda_MB_sat$w<-exp(-0.5*lambda_MB_sat$delta)/sum(exp(-0.5*lambda_MB_sat$delta))

# Exports a text file with the lambda value for each individual PGLS model
write.table(lambda_MB_sat,file = "lambda_MB_sat.txt")

# AIC-weighted lambda
lambda_MB_sat.lambda<-weighted.mean(lambda_MB_sat$lambda,lambda_MB_sat$w)