#                            COMPARATIVE ANALYSIS                             #

#S Cambreling
#23/06/25

#Comparative analyses of male reproductive senescence with pglmm

rm(list=ls())

{
  #Libraries --------------------------------------------------------------------
  library(ggplot2)
  library(lme4)
  library(MuMIn)
  library(Matrix)
  library(carData)
  library(Rmisc)
  library(lmerTest)
  library(lmtest)
  library(ape)
  library(caper)
  library(lmodel2)
  library(plotrix)
  library(ggeffects)
  library(cowplot)
  library(dplyr)
  library(optimx)
  library(nlme)
  library(rptR)
  library(segmented)
  library(ghibli)
  library(report)
  library(ggpubr)
  library(mgcv)
  library(stringr)
  library(phylobase)
  library(adephylo)
  library(MCMCglmm)
  library(orthopolynom)
  library(MuMIn)
  library(phytools)
  library(coda)
  library(phyr)
  library(rr2)
}

# Data -------------------------------------------------------------------------

# Function dredge_pglmm
# install.packages("pak")
pak::pak("LucasDLalande/selectools")

{
  dataset <- read.table("Data_meta_analysis.csv",sep=";",dec=".",header=TRUE) 

  # Adjusting Species for phylogeny Time tree
  dataset$Species <- dataset$Species_valid
  dataset[which(dataset$Species == "Bison_bison_athabascae"), "Species"] <- "Bison_bison"
  dataset[which(dataset$Species == "Diceros_bicornis_michaeli"), "Species"] <- "Diceros_bicornis"
  dataset[which(dataset$Species == "Rupicapra_rupicapra"), "Species"] <- "Rupicapra_rupicapra_rupicapra"
  
  # Phylogeny ------------------------------------------------------------------
  {
    tree_sen2=read.tree("phylogeny_timetree.nwk") 
    tree_sen2=as(tree_sen2, "phylo4")
    art <- as.data.frame(unique(dataset$Species))
    names(art)<-"Species"
    dim(art)
    row.names(art)=art[,1]
    a2=which(!(tree_sen2@label %in% row.names (art)))
    toto2=prune(tree_sen2,tip = a2)
    phylo_data2=phylo4d(toto2,tip.data = art)
    tree_sen2 <-as(phylo_data2,"phylo")
  }
  
  # Making the tree ultrametric for pglmms
  is.ultrametric(tree_sen2)
  tree_ultra <- force.ultrametric(tree_sen2, method=c("nnls","extend"))
  is.ultrametric(tree_ultra)
  
  # Variables --------------------------------------------------------------------
  {
    head(dataset)
    {
      dataset$afr <- as.numeric(dataset$afr)
      dataset$log_afr <- log(dataset$afr)
      dataset$log_mbm_testes <- log(dataset$mbm_testes)
      dataset$log_testes_mass <- log(dataset$testes_mass)
      dataset$rate <- as.numeric(dataset$rate)
      dataset$abs_rate <- abs(dataset$rate)
      dataset$log_abs_rate <- log(dataset$abs_rate)
      dataset$onset <- as.numeric(dataset$onset)
      dataset$log_onset <- log(dataset$onset)
      dataset$Species<- as.factor(dataset$Species)
      dataset$onset_sd <- as.numeric(dataset$onset_sd)
      dataset$rate_sd <- as.numeric(dataset$rate_sd)
      dataset$trait_std<- as.factor(dataset$trait_std)
      dataset$Species_valid<- as.factor(dataset$Species_valid)
      dataset$population<- as.factor(dataset$population)
      dataset$timing <- as.factor(dataset$timing)
      dataset$data_type <- as.factor(dataset$data_type)
      dataset$hunting_status <- as.factor(dataset$hunting_status)
      dataset$social_monogamy <- as.factor(dataset$social_monogamy)
      dataset$mass_article <- as.numeric(dataset$mass_article)
      dataset$mbm_ssd <- as.numeric(dataset$mbm_ssd)
      dataset$fbm_ssd <- as.numeric(dataset$fbm_ssd)
    }
    
    #Computing SSD
    dataset$ssd <- log(dataset$mbm_ssd/dataset$fbm_ssd)
    
    #Change longevity values for Marmota flaviventris, Urocitellus columbianus and Octodon_degus, as longevity observed in article is longer than the one extracted from bibliography
    dataset$longevity <- ifelse(dataset$Species=="Marmota_flaviventris", dataset$longevity_observed, dataset$longevity)
    dataset$longevity <- ifelse(dataset$Species=="Urocitellus_columbianus", dataset$longevity_observed, dataset$longevity)
    dataset$longevity <- ifelse(dataset$Species=="Octodon_degus", dataset$longevity_observed, dataset$longevity)
    
    
    #Variable longevity
    dataset$age_range <- (dataset$longevity_observed/dataset$longevity)
    dataset$age_range <- ifelse(dataset$age_range>1, 1, dataset$age_range) #if the lifespan coverage is above 1, reput to 1
    
    #Male body mass
    dataset$mass <- ifelse(is.na(dataset$mass_article), dataset$mbm_ssd, dataset$mass_article) #Take the mass depicted in the article, otherwise the male body mass used to compute SSD
    dataset$mass <- as.numeric(dataset$mass)
    dataset$log_mass <- log(dataset$mass)

    #Litter size included as number of offspring
    levels(dataset$trait_std)[levels(dataset$trait_std) == "litter_size"] <- "nb_off_breeders"
    dataset$trait_std <- droplevels(dataset$trait_std)

    #Scale dataset
    {
      dataset$log_afr_sc <- scale(dataset$log_afr)
      dataset$log_abs_rate_sc <- scale(dataset$log_abs_rate)
      dataset$log_onset_sc <- scale(dataset$log_onset)
      dataset$log_mass_sc <- scale(dataset$log_mass)
      dataset$ssd_sc <- scale(dataset$ssd)
      dataset$age_range_sc <- scale(dataset$age_range)
      dataset$log_testes_mass_sc <- scale(dataset$log_testes_mass)
      dataset$log_mbm_testes_sc  <- scale(dataset$log_mbm_testes)
    }
    
    #Subdataset, with only senescent species
    senesc <- subset(dataset,dataset$senescence!="0") #only senescent species
    
    # Removing NAs from dataset
    {
      senesc2<- senesc[!is.na(senesc$mass),] 
      senesc2<- senesc2[!is.na(senesc2$afr),] 
      senesc2<- senesc2[!is.na(senesc2$hunting_status),] 
      senesc2<- senesc2[!is.na(senesc2$data_type),] 
      senesc2<- senesc2[!is.na(senesc2$age_range),] 
      senesc2<- senesc2[!is.na(senesc2$onset),]  #removing age classes with no onset or rate
    }
  }
}


# PGLMMs ------------------------------------------------------------------------

# DETECTION OF SENESCENCE ------------------------------------------------------

#Number of offspring as reference trait 
dataset$trait_std <- relevel(dataset$trait_std, ref = "nb_off")

#Dredge function for pglmms
dredge_pglmm(formulaRE = "senescence ~ 1 + (1 | Species__) + (1 | population)",
                  fixed = c("log_afr_sc", "hunting_status","trait_std","data_type","age_range_sc"),
                  data = dataset,
                  rank="AICc",
                  family="binomial",
                  cov_ranef=list(Species = tree_ultra),
                  estimate = T,
                  std.err=T,
                  round=2)

# Selected model ---------------------------------------------------------------
(mod1_detec <- pglmm(senescence ~ age_range  +  (1|Species__) + (1|population), 
                     data = dataset, family = "binomial", cov_ranef = list(Species = tree_ultra),
                     REML = TRUE, verbose = F, s2.init = .1))


R2(mod1_detec) #R2 value of selected model

#Part of variance explained by phylogeny
var <- ranef(mod1_detec)
(h2_phylo <- (var["1|Species__",]$Variance + var["1|Species",]$Variance) / sum(var$Variance)) #heritability of phylogeny
(h2_pop <- (var["1|population",]$Variance) / sum(var$Variance)) #heritability of population

# Subsample with the number of observations ------------------------------------

dataset_nobs <- dataset[!is.na(dataset$nb_obs),]
dataset_nobs$nb_obs <- as.numeric(dataset_nobs$nb_obs)
dataset_nobs$log_nb_obs <- log(dataset_nobs$nb_obs)
dataset_nobs$log_nb_obs_sc <- scale(dataset_nobs$log_nb_obs)

dredge_pglmm(formulaRE = "senescence ~ 1 + (1 | Species__) + (1 | population)",
                  fixed = c("log_afr_sc", "hunting_status","trait_std","data_type","age_range_sc","log_nb_obs_sc"),
                  data = dataset_nobs,
                  rank="AICc",
                  family="binomial",
                  cov_ranef=list(Species = tree_ultra),
                  estimate = T,
                  std.err=T,
                  round=2)


# Selected model 
(mod2_detec <- pglmm(senescence ~  log_nb_obs + age_range +(1|Species__) + (1|population), 
                     data = dataset_nobs, family = "binomial", cov_ranef = list(Species = tree_ultra),
                     REML = TRUE, verbose = F, s2.init = .1))

R2(mod2_detec) #0.60

#Part of variance explained by phylogeny
var <- ranef(mod2_detec)
(h2_phylo <- (var["1|Species__",]$Variance + var["1|Species",]$Variance) / sum(var$Variance)) 
(h2_pop <- (var["1|population",]$Variance) / sum(var$Variance)) 

# ONSET OF SENESCENCE ----------------------------------------------------------

# h2 constant model, to assess the impact phylogenetic relatedness has when no factor is included

(onset_cst <- pglmm(log_onset ~ 1 + (1|Species__) + (1|population), 
                    data = senesc2, family = "gaussian", cov_ranef = list(Species = tree_ultra),
                    REML = TRUE, verbose = F, s2.init = .1))

#Part of variance explained by phylogeny
var <- ranef(onset_cst)
(h2_phylo <- (var["1|Species__",]$Variance + var["1|Species",]$Variance) / sum(var$Variance)) #Species
(h2_pop <- (var["1|Species__",]$Variance) / sum(var$Variance)) #Phylogeny
(h2_pop <- (var["1|population",]$Variance) / sum(var$Variance)) #Population
(h2_resid <- (var["residual",]$Variance) / sum(var$Variance)) #Residual


# With SSD and Mating system ---------------------------------------------------

senesc2$trait_std <- relevel(senesc2$trait_std, ref = "nb_off") #Number of offspring as reference trait

dredge_pglmm(formulaRE = "log_onset ~ 1 + (1 | Species__) + (1 | population)",
             fixed = c("log_afr_sc", "log_mass_sc","data_type", "hunting_status", "trait_std","ssd_sc","social_monogamy", "log_afr_sc*hunting_status","log_afr_sc*log_mass_sc"),
             data = senesc2,
             rank="AICc",
             family="gaussian",
             cov_ranef=list(Species = tree_ultra),
             estimate = T,
             std.err=T,
             round=2)

#Selected model

(mod1_onset <- pglmm(log_onset ~ log_afr + (1|Species__) + (1|population), 
                     data = senesc2, family = "gaussian", cov_ranef = list(Species = tree_ultra),
                     REML = TRUE, verbose = F, s2.init = .1))

#Part of variance explained by phylogeny
var <- ranef(mod1_onset)
(h2_spe <- (var["1|Species__",]$Variance + var["1|Species",]$Variance) / sum(var$Variance)) #Species
(h2_phy <- (var["1|Species__",]$Variance) / sum(var$Variance)) #Phylogeny
(h2_pop <- (var["1|population",]$Variance) / sum(var$Variance)) #Population
(h2_resid <- (var["residual",]$Variance) / sum(var$Variance)) #Residual

# Relative testes mass ---------------------------------------------------------
senesc3<- subset(senesc2,senesc2$log_testes_mass_sc!="NA") 
senesc3<- subset(senesc3,senesc3$log_mbm_testes_sc!="NA") 

dredge_pglmm(formulaRE = "log_onset ~ 1 + (1 | Species__) + (1 | population)",
                  fixed = c("log_afr_sc", "data_type","hunting_status", "log_testes_mass_sc ","log_mbm_testes_sc ", "trait_std",  "log_afr_sc*hunting_status"),
                  data = senesc3,
                  rank="AICc",
                  family="gaussian",
                  cov_ranef=list(Species = tree_ultra),
                  estimate = T,
                  std.err=T,
                  round=2)

# Analyses on traits related to reproductive success only ----------------------

# Restricting dataset to reproductive traits only
data_RS <- senesc2[!senesc2$trait_std %in% c("cop", "mate"), ]
data_RS$trait_std <- droplevels(data_RS$trait_std)

# With ssd and mating system 
dredge_pglmm(formulaRE = "log_onset ~ 1 + (1 | Species__) + (1 | population)",
                  fixed = c("log_afr_sc","hunting_status","log_mass_sc","ssd","social_monogamy","log_afr_sc*log_mass_sc","log_afr_sc*hunting_status","trait_std"),
                  data = data_RS, 
                  rank="AICc",
                  family="gaussian",
                  cov_ranef=list(Species = tree_ultra),
                  estimate = T,
                  std.err=T,
                  round=2)


# With Relative Testes Mass
data_RS2<- subset(data_RS,data_RS$log_testes_mass_sc!="NA") 
data_RS2<- subset(data_RS2,data_RS2$log_mbm_testes_sc!="NA") 

dredge_pglmm(formulaRE = "log_onset ~ 1 + (1 | Species__) + (1 | population)",
                  fixed = c("log_afr_sc", "data_type","hunting_status", "log_testes_mass_sc ","log_mbm_testes_sc ", "trait_std",  "log_afr_sc*hunting_status"),
                  data = data_RS2,
                  rank="AICc",
                  family="gaussian",
                  cov_ranef=list(Species = tree_ultra),
                  estimate = T,
                  std.err=T,
                  round=2)


# RATE OF SENESCENCE -----------------------------------------------------------

# Number of offspring ----------------------------------------------------------
senesc_no <- senesc2[senesc2$trait_std %in% c("nb_off"), ]

dredge_pglmm(formulaRE = "log_abs_rate ~ 1 + (1 | Species__) + (1 | population)",
                  fixed = c("log_afr","log_mass","ssd","social_monogamy","log_afr*log_mass"),
                  data = senesc_no, 
                  rank="AICc",
                  family="gaussian",
                  cov_ranef=list(Species = tree_ultra),
                  estimate = T,
                  std.err=T,
                  round=2)


# Selected Model 
(mod2_no <- pglmm(log_abs_rate ~  log_afr*log_mass + ssd + (1|Species__) + (1|population), 
                  data = senesc_no, family = "gaussian", cov_ranef = list(Species = tree_ultra),
                  REML = TRUE, verbose = F, s2.init = .1))

R2(mod2_no) 

#Part of variance explained by phylogeny
var <- ranef(mod1_no)
(h2_phylo <- (var["1|Species__",]$Variance + var["1|Species",]$Variance) / sum(var$Variance)) 
(h2_pop <- (var["1|population",]$Variance) / sum(var$Variance)) 
(h2_resid <- (var["residual",]$Variance) / sum(var$Variance)) 

# Relative testes mass

dredge_pglmm(formulaRE = "log_abs_rate ~ 1 + (1 | Species__) + (1 | population)",
                  fixed = c("log_afr_sc","log_testes_mass_sc","log_mbm_testes_sc"),
                  data = senesc_no, 
                  rank="AICc",
                  family="gaussian",
                  cov_ranef=list(Species = tree_ultra),
                  estimate = T,
                  std.err=T,
                  round=2)

# SCALING TRAITS (on reproductive traits) --------------------------------------

#Remove traits probability(mating) and number of copulations 
senesc4 <- senesc2[senesc2$trait_std %in% c("nb_off", "patern","nb_off_breeders"), ]

#Rates of senescence centered and scaled for each trait
senesc_rate <- senesc4 %>%
  group_by(trait_std) %>%
  mutate_at(vars(log_abs_rate,log_abs_rate_sc), scale)

# h2 constant model

(rate_cst <- pglmm(log_abs_rate ~ 1 + (1|Species__) + (1|population), 
                   data = senesc_rate, family = "gaussian", cov_ranef = list(Species = tree_ultra),
                   REML = TRUE, verbose = F, s2.init = .1))

#Part of variance explained by phylogeny
var <- ranef(rate_cst)
(h2_spe <- (var["1|Species__",]$Variance + var["1|Species",]$Variance) / sum(var$Variance)) #Species
(h2_phy <- (var["1|Species__",]$Variance) / sum(var$Variance)) #Phylogeny
(h2_pop <- (var["1|population",]$Variance) / sum(var$Variance)) #Population
(h2_resid <- (var["residual",]$Variance) / sum(var$Variance)) #Residual


# With SSD and Mating system ---------------------------------------------------
senesc_rate$trait_std <- droplevels(senesc_rate$trait_std)

dredge_pglmm(formulaRE = "log_abs_rate ~ 1 + (1 | Species__) + (1 | population)",
                  fixed = c("log_afr_sc","hunting_status","log_mass_sc","ssd","social_monogamy","log_afr_sc*log_mass_sc","log_afr_sc*hunting_status","trait_std"),
                  data = senesc_rate, 
                  rank="AICc",
                  family="gaussian",
                  cov_ranef=list(Species = tree_ultra),
                  estimate = T,
                  std.err=T,
                  round=2)


# Selected model 
(mod2_rate <- pglmm(log_abs_rate ~  log_afr + ssd + (1|Species__) + (1|population), 
                    data = senesc_rate, family = "gaussian", cov_ranef = list(Species = tree_ultra),
                    REML = TRUE, verbose = F, s2.init = .1))


R2(mod2_rate) 
#Part of variance explained by phylogeny
var <- ranef(mod2_rate)
(h2_spe <- (var["1|Species__",]$Variance + var["1|Species",]$Variance) / sum(var$Variance)) #Species
(h2_phy <- (var["1|Species__",]$Variance) / sum(var$Variance)) #Phylogeny
(h2_pop <- (var["1|population",]$Variance) / sum(var$Variance)) #Population
(h2_resid <- (var["residual",]$Variance) / sum(var$Variance)) #Residual

# Relative testes mass ---------------------------------------------------------

senesc_rate2<- subset(senesc_rate,senesc_rate$log_testes_mass_sc!="NA") 
senesc_rate2<- subset(senesc_rate2,senesc_rate2$log_mbm_testes_sc!="NA") 

dredge_pglmm(formulaRE = "log_abs_rate ~ 1 + (1 | Species__) + (1 | population)",
                  fixed = c("log_afr_sc", "hunting_status", "log_afr_sc*hunting_status","log_testes_mass_sc","log_mbm_testes_sc"),
                  data = senesc_rate2,
                  rank="AICc",
                  family="gaussian",
                  cov_ranef=list(Species = tree_ultra),
                  estimate = T,
                  std.err=T,
                  round=2)

# TRAIT TIMING -----------------------------------------------------------------
# Analyses checking if adding the timing of sampling of offspring would improve our previously selected models

# Senescence detection ---------------------------------------------------------

data_timing<- dataset[!is.na(dataset$timing),] 
data_timing$timing <- droplevels(data_timing$timing)

dredge_pglmm(formulaRE = "senescence ~ 1 + (1 | Species__) + (1 | population)",
                  fixed = c("age_range","timing"),
                  data = data_timing, 
                  rank="AICc",
                  family="gaussian",
                  cov_ranef=list(Species = tree_ultra),
                  estimate = T,
                  std.err=T,
                  round=2)

(mod_detec_timing <- pglmm(senescence ~ age_range  + timing +  (1|Species__) + (1|population), 
                     data = data_timing, family = "binomial", cov_ranef = list(Species = tree_ultra),
                     REML = TRUE, verbose = F, s2.init = .1))


# ONSET OF SENESCENCE ----------------------------------------------------------

senesc_timing<- senesc2[!is.na(senesc2$timing),] 
senesc_timing$timing <- droplevels(senesc_timing$timing)
senesc_timing$timing <- relevel(senesc_timing$timing, ref = "copulatory")

dredge_pglmm(formulaRE = "log_onset ~ 1 + (1 | Species__) + (1 | population)",
                  fixed = c("log_afr","timing"),
                  data = senesc_timing, 
                  rank="AICc",
                  family="gaussian",
                  cov_ranef=list(Species = tree_ultra),
                  estimate = T,
                  std.err=T,
                  round=2)

(mod_o_timing <- pglmm(log_onset ~ log_afr +timing + (1|Species__) + (1|population), 
                       data = senesc_timing, family = "gaussian", cov_ranef = list(Species = tree_ultra),
                       REML = TRUE, verbose = F, s2.init = .1))


# RATE OF SENESCENCE -----------------------------------------------------------

rate_timing<- senesc_rate[!is.na(senesc_rate$timing),] 
rate_timing$timing <- droplevels(rate_timing$timing)

dredge_pglmm(formulaRE = "log_abs_rate ~ 1 + (1 | Species__) + (1 | population)",
                  fixed = c("log_afr","ssd","timing"),
                  data = rate_timing, 
                  rank="AICc",
                  family="gaussian",
                  cov_ranef=list(Species = tree_ultra),
                  estimate = T,
                  std.err=T,
                  round=2)

# Selected model 
(mod_rate_timing <- pglmm(log_abs_rate ~  log_afr + ssd + timing + (1|Species__) + (1|population), 
                    data = rate_timing, family = "gaussian", cov_ranef = list(Species = tree_ultra),
                    REML = TRUE, verbose = F, s2.init = .1))


# BAYESIAN ANALYSES ------------------------------------------------------------
#Bayesian analyses (using MCMCglmm) to add the uncertainty of the onset and the rate in the analyses, in the previously selected models

# Onset of senescence ----------------------------------------------------------

# Bayesian MCMCglmm

#Phylogenetic tree
tree_ultra$edge.length[tree_ultra$edge.length < 1e-8] <- 1e-8 #replacing branches at 0 by very small branches to ensure stability during the matrix inversion

#Setting priors
a <- 1000
prior <- list(R = list(V = 1, nu = 0.001),G = list(G1 = list(V = 1, nu = 0.001),
                                                   G2 = list(V = 1, nu = 0.001),G3 = list(V = 1, nu = 0.001)))


inv.phylo<-inverseA(tree_ultra,nodes="TIPS" #or ALL for a larger phylogeny more efficient 
                    ,scale=TRUE)
senesc2$Species2 <- senesc2$Species #need to create a duplicate for the phylogenetic and non-phylogenetic related species effects

# Testing MCMCglmm models without standard errors of onsets to check the values and compare with frequentist models
model_no_SE_onset <- MCMCglmm(log_onset ~ log_afr # model formula here null model
                              , random = ~Species+Species2+population # add here random effect of species and phylogeny
                              , family = "gaussian" # distribution of trait
                              , ginverse = list(Species = inv.phylo$Ainv) # include a custom matrix for argument phylo
                              , prior = prior
                              , data = senesc2
                              , nitt = 4e+05 # number of iteration after burnin
                              , burnin = 1000 # number of iteration before beginning sample
                              , thin = 50 # nb of iteration between sample
                              , pr=TRUE) #save random posterior distribution

plot(model_no_SE_onset) #check convergence
summary(model_no_SE_onset)

#Adding SE for first age models

min(senesc2$onset_sd, na.rm = TRUE) #0.642021
senesc2$onset_sd[is.na(senesc2$onset_sd)] <- 0.0641749061969476
senesc2$onset_sd <- as.numeric(senesc2$onset_sd)


#formula only works for log tranformed data (as we use log values of onset we need log values of the variance)
onset_var <- senesc2$onset_sd^2/senesc2$onset^2

model_SE_onset <- MCMCglmm(log_onset ~ log_afr# model formula here null model
                           , random = ~Species+Species2+population # add here random effect of species and phylogeny
                           , family = "gaussian" # distribution of trait
                           , mev = onset_var # error variance associated to each data point
                           , ginverse = list(Species = inv.phylo$Ainv) # include a custom matrix for argument phylo
                           , prior = prior
                           , data = senesc2
                           , nitt = 4e+05 # number of iteration after burnin
                           , burnin = 1000 # number of iteration before beginning sample
                           , thin = 50 # nb of iteration between sample
                           , pr=TRUE) #save random posterior distribution

plot(model_SE_onset)
summary(model_SE_onset)

# Rate of senescence -----------------------------------------------------------

# Number of offspring --------------------------------------------------------
# Bayesian
a <- 1000
prior <- list(R = list(V = 1, nu = 0.001),G = list(G1 = list(V = 1, nu = 0.001),
                                                   G2 = list(V = 1, nu = 0.001),G3 = list(V = 1, nu = 0.001)))
inv.phylo<-inverseA(tree_ultra,nodes="TIPS" #or ALL for a larger phylogeny more efficient 
                    ,scale=TRUE)

senesc_no$Species2 <- senesc_no$Species

# Testing MCMCglmm models without standard errors of rates to check the values and compare with frequentist models
model_NO_no_SE <- MCMCglmm(log_abs_rate ~ log_afr*log_mass + ssd# model formula here null model
                        , random = ~Species+Species2+population # add here random effect of species and phylogeny
                        , family = "gaussian" # distribution of trait
                        , ginverse = list(Species = inv.phylo$Ainv) # include a custom matrix for argument phylo
                        , prior = prior
                        , data = senesc_no
                        , nitt = 4e+05 # number of iteration after burnin
                        , burnin = 1000 # number of iteration before beginning sample
                        , thin = 50 # nb of iteration between sample
                        , pr=TRUE) #save random posterior distribution
plot(model_NO_no_SE)
summary(model_NO_no_SE)

#formula only works for log tranformed data
slope_var <- senesc_no$rate_sd^2/senesc_no$abs_rate^2
model_NO_SE <- MCMCglmm(log_abs_rate ~ log_afr*log_mass + ssd # model formula here null model
                     , random = ~Species+Species2+population # add here random effect of species and phylogeny
                     , family = "gaussian" # distribution of trait
                     , mev = slope_var # error variance associated to each data point
                     , ginverse = list(Species = inv.phylo$Ainv) # include a custom matrix for argument phylo
                     , prior = prior
                     , data = senesc_no
                     , nitt = 4e+05 # number of iteration after burnin
                     , burnin = 1000 # number of iteration before beginning sample
                     , thin = 50 # nb of iteration between sample
                     , pr=TRUE) #save random posterior distribution
plot(model_SE)
summary(model_NO_SE)


# Scaled traits ---------------------------------------------------------

senesc_rate$Species2 <- senesc_rate$Species #need to create a duplicate for the phylogenetic and non-phylogenetic related species effects

#Reload data, because bug
{
  # write.csv(senesc_rate, file = "dataset_rate_mcmcglmm.csv", row.names = FALSE) #need to extract data with scaled traits otherwise bugging in MCMCglmm
  dataset_rate <- read.table("dataset_rate_mcmcglmm.csv",sep=",",dec=".",header=TRUE) 
  #Phylogenetic tree
  {
    tree_sen2=read.tree("phylogeny_final_timetree.nwk") 
    tree_sen2=as(tree_sen2, "phylo4")
    art <- as.data.frame(unique(dataset_rate$Species))
    names(art)<-"Species"
    dim(art)
    row.names(art)=art[,1]
    a2=which(!(tree_sen2@label %in% row.names (art)))
    toto2=prune(tree_sen2,tip = a2)
    phylo_data2=phylo4d(toto2,tip.data = art)
    tree_sen2 <-as(phylo_data2,"phylo")
  }
  
  # Ultrametric tree
  is.ultrametric(tree_sen2)
  tree_ultra <- force.ultrametric(tree_sen2, method=c("nnls","extend"))
  is.ultrametric(tree_ultra)
  tree_ultra$edge.length[tree_ultra$edge.length < 1e-8] <- 1e-8 #replacing branches at 0 by very small branches to ensure stability during the matrix inversion
}


# Bayesian
a <- 1000
prior <- list(R = list(V = 1, nu = 0.001),G = list(G1 = list(V = 1, nu = 0.001),
                                                   G2 = list(V = 1, nu = 0.001),G3 = list(V = 1, nu = 0.001)))
inv.phylo<-inverseA(tree_ultra,nodes="TIPS" #or ALL for a larger phylogeny more efficient 
                    ,scale=TRUE)

# Testing MCMCglmm models without standard errors of rates to check the values and compare with frequentist models
model_rate_no_SE <- MCMCglmm(log_abs_rate ~ log_afr + ssd# model formula here null model
                        , random = ~Species+ Species2 +population # add here random effect of species and phylogeny
                        , family = "gaussian" # distribution of trait
                        , ginverse = list(Species = inv.phylo$Ainv) # include a custom matrix for argument phylo
                        , prior = prior
                        , data = dataset_rate
                        , nitt = 4e+05 # number of iteration after burnin
                        , burnin = 1000 # number of iteration before beginning sample
                        , thin = 50 # nb of iteration between sample
                        , pr=TRUE) #save random posterior distribution

plot(model_rate_no_SE)
summary(model_rate_no_SE)

#formula only works for log tranformed data
slope_var <- dataset_rate$rate_sd^2/dataset_rate$abs_rate^2

model_rate_SE <- MCMCglmm(log_abs_rate ~ log_afr+ ssd # model formula here null model
                     , random = ~Species+Species2+population # add here random effect of species and phylogeny
                     , family = "gaussian" # distribution of trait
                     , mev = slope_var # error variance associated to each data point
                     , ginverse = list(Species = inv.phylo$Ainv) # include a custom matrix for argument phylo
                     , prior = prior
                     , data = dataset_rate
                     , nitt = 4e+05 # number of iteration after burnin
                     , burnin = 1000 # number of iteration before beginning sample
                     , thin = 50 # nb of iteration between sample
                     , pr=TRUE) #save random posterior distribution
plot(model_rate_SE)
summary(model_rate_SE)