#R script for running all models presented in the paper pacman::p_load(lme4, MCMCglmm, dplyr, DescTools, metafor, matrixcalc) #Import data dat <- read.csv("DataFile.csv") #convert file to csv before running this code #preparing the data---- # response dat$zr_f <- dat$zr_fem*dat$good_fem dat$zr_m <- dat$zr_male*dat$good_male # rename - we cannot have "trait" in the data set dat <- rename(dat, c("trait2" = "trait")) # missingness is not allowed in error so we assume males and females will have similar errors when this info is missing for one of the sexes dat$v_fem <- ifelse(is.na(dat$zr_se_fem), dat$zr_se_male^2, dat$zr_se_fem^2) dat$v_male <- ifelse(is.na(dat$zr_se_male), dat$zr_se_fem^2, dat$zr_se_male^2) #phylogenetic data tree <- read.nexus(file = "output.nex") phylo <- tolower(dat$scientific_names) phylo <- StrCap(phylo, method=c("first")) phylo <- gsub(" ", "_", phylo) dat[, "phylo"] <- phylo dat$phylo[dat$phylo=="Cyanistes_caeruleus"] <-"Parus_caeruleus" dat$phylo[dat$phylo=="Spinus_tristis"] <-"Carduelis_tristis" dat$phylo[dat$phylo=="Aquila_pennata"] <-"Hieraaetus_pennatus" dat$phylo[dat$phylo=="Setophaga_petechia"] <-"Dendroica_petechia" dat$animal <- dat$phylo #analysing sexual dimorphism using Cohen's d effect sizes data.dimorphism <- subset(dat, dimorphism_Cohend!="NA") # Specify data set (no missing values) median(data.dimorphism$dimorphism_Cohend) range(data.dimorphism$dimorphism_Cohend) vc_matrix <- matrix(0,nrow = dim(data.dimorphism)[1],ncol = dim(data.dimorphism)[1]) diag(vc_matrix) <- (data.dimorphism$dimorphism_CohendSE)^2 check1 <- is.positive.definite(vc_matrix) # needs to be positive invVC <- solve(vc_matrix) invVC <- as(invVC,"dgCMatrix") rownames(invVC)=(1:length(data.dimorphism$dimorphism_CohendSE)) colnames(invVC)=(1:length(data.dimorphism$dimorphism_CohendSE)) data.dimorphism$obs<-(1:length(data.dimorphism$dimorphism_CohendSE)) prior <- list(R = list (V = 1,n = 0.002), G = list (G1= list(V=1, nu=1, alpha.mu = 0, alpha.V = 1000), G2= list(V=1, nu=1, alpha.mu = 0, alpha.V = 1000), G3= list(V=1,fix=1)) ) thin <- 100 burn.in <- 1000 n.itt <- thin*1000+burn.in modelos <- NULL # setting up an empty list to store models class(modelos) <- 'list' #loop to control for phylogenetic effects for (i in 1:50){ AinvTA.2<-inverseA(tree[[i]])$Ainv mc <- MCMCglmm(fixed=dimorphism_Cohend~1, random=~scientific_names+animal+obs, ginverse=list(animal = AinvTA.2, obs = invVC), data=data.dimorphism, nitt=n.itt, thin=thin, burnin=burn.in, verbose=T, prior=prior, pr=T) modelos[[i]] <- mc print(i) } summary(modelos[[1]]) #output of the first model (example) #processing multimodel results---- #combining posterior distributions across models all.m.vcv <- modelos[[1]]$VCV all.m.sol <- modelos[[1]]$Sol for (i in 2:50){ all.m.vcv <- rbind(all.m.vcv, modelos[[i]]$VCV) all.m.sol <- rbind(all.m.sol, modelos[[i]]$Sol) } #converting into mcmc objects all.m.vcv <- as.mcmc(all.m.vcv) all.m.sol <- as.mcmc(all.m.sol) #combine output into a list---- #this includes: #(a) a list with all 50 models #(b) the combined posterior distribution of all sol in (a) #(c) the combined posterior distribution of all vcv in (a) #*************************************************************************************** output.model.dimorphism <- (list(models=modelos, sol=all.m.sol, vcv=all.m.vcv)) # SPECIFY OUTPUT NAME #*************************************************************************************** #remove temporary files rm(modelos, all.m.sol, all.m.vcv) #summarising results---- #*************************************************************************************** mod <- output.model.dimorphism # SPECIFY NAME OF THE MODEL THAT WE WANT TO SUMMARISE #*************************************************************************************** #fixed effects #posterior means for the intercept (i.e. meta-analytical mean for the degree of sexual dimorphism) mean(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) #95%CIs HPDinterval(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) #random effects #posterior means colMeans(mod$vcv) #95%CIs HPDinterval(mod$vcv) #computing effective sample sizes effectiveSize(mod$vcv) effectiveSize(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) #Bivariate meta-analysis---- #MODEL 0 - MODEL INDLUDING SEX AND STUDY TYPE (correlational/experimental)---- # prior prior <- list(R = list(V = diag(2), nu = 2), G = list(G1 = list(V = diag(2) , nu = 2, alpha.mu=c(0,0), alpha.V=diag(2)*25^2), G2 = list(V = diag(2) , nu = 2, alpha.mu=c(0,0), alpha.V=diag(2)*25^2))) modelos <- NULL # setting up an empty list to store models class(modelos) <- 'list' #******************************** datos <- dat # SPECIFY DATA SET #******************************** #run loop that iterates model across 50 phylogenies to incorporate phylogenetic uncertainty. for (i in 1:50){ AinvTA.2<-inverseA(tree[[i]])$Ainv m1 <- MCMCglmm(cbind(zr_f, zr_m) ~ trait + study_type + scale(year, scale= FALSE) - 1, random = ~ us(trait):animal + us(trait):scientific_names, rcov = ~ us(trait):units, mev = c(datos$v_fem, datos$v_male), ginverse=list(animal = AinvTA.2), data = datos, family = rep('gaussian', 2), nitt=11000, thin=100, burnin=1000, pl = TRUE, pr = TRUE, prior = prior, verbose=FALSE) modelos[[i]] <- m1 print(i) } summary(modelos[[1]]) #eg output of the first model #processing multimodel results---- #combining posterior distributions across models all.m.vcv <- modelos[[1]]$VCV all.m.sol <- modelos[[1]]$Sol for (i in 2:50){ all.m.vcv <- rbind(all.m.vcv, modelos[[i]]$VCV) all.m.sol <- rbind(all.m.sol, modelos[[i]]$Sol) } #converting into mcmc objects all.m.vcv <- as.mcmc(all.m.vcv) all.m.sol <- as.mcmc(all.m.sol) #combine output into a list---- #this includes: #(a) a list with all 50 models #(b) the combined posterior distribution of all sol in (a) #(c) the combined posterior distribution of all vcv in (a) #*************************************************************************************** output.model0 <- (list(models=modelos, sol=all.m.sol, vcv=all.m.vcv)) # SPECIFY OUTPUT NAME #*************************************************************************************** #remove temporary files rm(modelos, all.m.sol, all.m.vcv) #summarising results---- #*************************************************************************************** mod <- output.model0 # SPECIFY NAME OF MODEL WE WANT TO SUMMARISE #*************************************************************************************** #fixed effects #posterior means colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) #95%CIs HPDinterval(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) #random effects #posterior means colMeans(mod$vcv) #95%CIs HPDinterval(mod$vcv) #computing effective sample sizes effectiveSize(mod$vcv) effectiveSize(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) results.table <- rbind( #sol results data.frame(parameter=names(colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl])), post.mean=colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl]), HPDinterval(mod$sol[,1:mod$models[[1]]$Fixed$nfl]), eff.sample.size=effectiveSize(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) ), #vcv results data.frame(parameter=names(colMeans(mod$vcv)), post.mean=colMeans(mod$vcv), HPDinterval(mod$vcv), eff.sample.size=effectiveSize(mod$vcv)) ) results.table <- subset(results.table, parameter!="sqrt(mev):sqrt(mev).meta") # eliminates row for measurement error rownames(results.table) <- NULL #********************************************************************** table.0 <- results.table # give results.table a name #********************************************************************** #MODEL 1 - GENERAL MODEL WITH SEX AS COVARIATE---- # prior prior <- list(R = list(V = diag(2), nu = 2), G = list(G1 = list(V = diag(2) , nu = 2, alpha.mu=c(0,0), alpha.V=diag(2)*25^2), G2 = list(V = diag(2) , nu = 2, alpha.mu=c(0,0), alpha.V=diag(2)*25^2))) modelos <- NULL # setting up an empty list to store models class(modelos) <- 'list' #******************************** datos <- dat # SPECIFY DATA SET #******************************** #run loop. Model with only sex as covariate for (i in 1:50){ AinvTA.2<-inverseA(tree[[i]])$Ainv m1 <- MCMCglmm(cbind(zr_f, zr_m) ~ trait + scale(year, scale= FALSE) - 1, random = ~ us(trait):animal + us(trait):scientific_names, rcov = ~ us(trait):units, mev = c(datos$v_fem, datos$v_male), ginverse=list(animal = AinvTA.2), data = datos, family = rep('gaussian', 2), nitt=11000, thin=100, burnin=1000, pl = TRUE, pr = TRUE, prior = prior, verbose=FALSE) modelos[[i]] <- m1 print(i) } summary(modelos[[1]]) #processing multimodel results---- #combining posterior distributions across models all.m.vcv <- modelos[[1]]$VCV all.m.sol <- modelos[[1]]$Sol for (i in 2:50){ all.m.vcv <- rbind(all.m.vcv, modelos[[i]]$VCV) all.m.sol <- rbind(all.m.sol, modelos[[i]]$Sol) } #converting into mcmc objects all.m.vcv <- as.mcmc(all.m.vcv) all.m.sol <- as.mcmc(all.m.sol) #combine output into a list---- #this includes: #(a) a list with all 50 models #(b) the combined posterior distribution of all sol in (a) #(c) the combined posterior distribution of all vcv in (a) #*************************************************************************************** output.model1 <- (list(models=modelos, sol=all.m.sol, vcv=all.m.vcv)) # SPECIFY OUTPUT NAME #*************************************************************************************** #remove temporary files rm(modelos, all.m.sol, all.m.vcv) #summarising results---- #*************************************************************************************** mod <- output.model1 # SPECIFY NAME OF MODEL THAT WE WANT TO SUMMARISE #*************************************************************************************** #Calculating Heterogeneity I2 # first we get typical sampling variance colnames(mod$vcv) #For females sampling_variance_fem <- sum(1 / dat$v_fem) * (length(dat$v_fem) - 1) / (sum(1 / dat$v_fem)^2 - sum((1 / dat$v_fem)^2)) #Total heterogeneity for females: total_fem <- (mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10])/ (mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10] + sampling_variance_fem) mean(total_fem) HPDinterval(total_fem) #Variance for phylogeny given that VCV[,1] is phylogeny phylo_fem <- mod$vcv[,1]/ (mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10] + sampling_variance_fem) mean(phylo_fem) HPDinterval(phylo_fem) #For species ID spp_fem <- mod$vcv[,5]/ (mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10] + sampling_variance_fem) mean(spp_fem) HPDinterval(spp_fem) #Residual res_fem<-mod$vcv[,10]/ (mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10] + sampling_variance_fem) mean(res_fem) HPDinterval(res_fem) #For males sampling_variance_male <- sum(1 / dat$v_male) * (length(dat$v_male) - 1) / (sum(1 / dat$v_male)^2 - sum((1 / dat$v_male)^2)) #Total heterogeneity for males: total_male<-(mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13])/ (mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13] + sampling_variance_male) mean(total_male) HPDinterval(total_male) #I2 for phylogeny phylo_male <-mod$vcv[,4]/ (mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13] + sampling_variance_male) mean(phylo_male) HPDinterval(phylo_male) #For species ID spp_male<-mod$vcv[,8]/ (mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13] + sampling_variance_male) mean(spp_male) HPDinterval(spp_male) #Residual res_male<-mod$vcv[,13]/ (mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13] + sampling_variance_male) mean(res_male) HPDinterval(res_male) #fixed effects #posterior means colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) #95%CIs HPDinterval(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) #random effects #posterior means colMeans(mod$vcv) #95%CIs HPDinterval(mod$vcv) #computing effective sample sizes effectiveSize(mod$vcv) effectiveSize(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) results.table <- rbind( #sol results data.frame(parameter=names(colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl])), post.mean=colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl]), HPDinterval(mod$sol[,1:mod$models[[1]]$Fixed$nfl]), eff.sample.size=effectiveSize(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) ), #vcv results data.frame(parameter=names(colMeans(mod$vcv)), post.mean=colMeans(mod$vcv), HPDinterval(mod$vcv), eff.sample.size=effectiveSize(mod$vcv)) ) results.table <- subset(results.table, parameter!="sqrt(mev):sqrt(mev).meta") # eliminates row for measurement error rownames(results.table) <- NULL #********************************************************************** table.1 <- results.table # give results.table a name #********************************************************************** #Differences between sexes colMeans(mod$sol[,1:2]) sex.diff <- mod$sol[,1]-mod$sol[,2] mean(sex.diff) HPDinterval(sex.diff) #Proportion of the variance explained by phylogeny and species id # phylogeny prop.phyl.var.fem <- mod$vcv[,1]/ (mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10]) #females mean(prop.phyl.var.fem) HPDinterval(prop.phyl.var.fem) prop.phyl.var.male <- mod$vcv[,4]/ (mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13]) #males mean(prop.phyl.var.male) HPDinterval(prop.phyl.var.male) # species id prop.spp.var.fem <- mod$vcv[,5]/ (mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10]) #females mean(prop.spp.var.fem) HPDinterval(prop.spp.var.fem) prop.spp.var.male <- mod$vcv[,8]/ (mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13]) #males mean(prop.spp.var.male) HPDinterval(prop.spp.var.male) #residual prop.residual.fem <- mod$vcv[,10]/(mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10]) mean(prop.residual.fem) HPDinterval(prop.residual.fem) prop.residual.male <- mod$vcv[,13]/(mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13]) mean(prop.residual.male) HPDinterval(prop.residual.male) #correlations for phylogeny (transformed from covariances and sex specific variances) #getting correlations - (correlation = cov/(SD*SD)) spp_mat <- mod$vcv[,c(1, 3, 4)] # covariance components attr(spp_mat, "dimnames")[[2]]<- c("f_V", "fm_cov", "m_V") # female variance, covariance, male variance spp_cor <- spp_mat[, "fm_cov"]/(sqrt(spp_mat[, "f_V"]*spp_mat[, "m_V"])) phylo.cor <- spp_cor mean(phylo.cor) HPDinterval(phylo.cor) summary(phylo.cor) #slope for phylogeny #computed as covariance divided by male variance slope_phylo <- spp_mat[,2]/spp_mat[,3] summary(slope_phylo) #correlations for species id #getting correlations - spp (correlation = cov/(SD*SD)) spp_mat <- mod$vcv[,c(5, 7, 8)] attr(spp_mat, "dimnames")[[2]]<- c("f_V", "fm_cov", "m_V") # female, male variance and cov spp_cor <- spp_mat[, "fm_cov"]/(sqrt(spp_mat[, "f_V"]*spp_mat[, "m_V"])) sp.id.cor <- spp_cor mean(sp.id.cor) HPDinterval(sp.id.cor) summary(sp.id.cor) #slope for species #computed as covariance divided by male variance slope_species <- spp_mat[,2]/spp_mat[,3] summary(slope_species) #correlations for residual #getting correlations - spp (correlation = cov/(SD*SD)) spp_mat <- mod$vcv[,c(10, 12, 13)] attr(spp_mat, "dimnames")[[2]]<- c("f_V", "fm_cov", "m_V") # female, male variance and cov spp_cor <- spp_mat[, "fm_cov"]/(sqrt(spp_mat[, "f_V"]*spp_mat[, "m_V"])) res.cor <- spp_cor mean(res.cor) HPDinterval(res.cor) summary(res.cor) #slope for units #computed as covariance divided by male variance slope_units <- spp_mat[,2]/spp_mat[,3] summary(slope_units) #Publication bias---- # predictions for fixed effects pred_fix <- as.vector(apply(mod$sol,2,mean)[1:mod$models[[1]]$Fixed$nfl] %*% t(mod$models[[1]]$X)) # predictions for random effects pred_ran <- as.vector(apply(mod$sol,2,mean)[-(1:mod$models[[1]]$Fixed$nfl)] %*% t(mod$models[[1]]$Z)) # predictions for sampling variance raw_se <- c(datos$zr_se_fem, datos$zr_se_male) #pred_mev <- apply(mod$sol,2,mean)[(ncol(mod$sol)-nrow(datos)*2 + 1):ncol(mod$sol)]*raw_se # observed data raw_data <- c(datos$zr_f, datos$zr_m) # residuals pred_res <- raw_data - (pred_fix + pred_ran) # precision precision <- c(1/datos$zr_se_fem, 1/datos$zr_se_male) #funnel plots - the first 737 for females and the rest for males #females plot(pred_res[1:737], precision[1:737]) #males plot(pred_res[738:1474], precision[738:1474]) #female residuals residuales <- pred_res[1:737] precis <- precision[1:737] #Trim and fill procedures metaR <- rma(yi= residuales, sei = 1/precis) TF <-trimfill(metaR,estimator="R0") TF funnel(TF, 'seinv') #Egger's test regtest(metaR, model="lm", predictor="sei") #male residuals residuales <- pred_res[738:1474] precis <- precision[738:1474] #Trim and fill procedures metaR <- rma(yi= residuales, sei = 1/precis) TF <-trimfill(metaR,estimator="R0") TF funnel(TF, 'seinv') #Egger's test regtest(metaR, model="lm", predictor="sei") #MODEL 2 - MODEL WITH INTERACTION BETWEEN SEX AND MAJOR PARAMETER (CONDITION/FITNESS) ---- # prior prior <- list(R = list(V = diag(2), nu = 2), G = list(G1 = list(V = diag(2) , nu = 2, alpha.mu=c(0,0), alpha.V=diag(2)*25^2), G2 = list(V = diag(2) , nu = 2, alpha.mu=c(0,0), alpha.V=diag(2)*25^2))) modelos <- NULL class(modelos) <- 'list' #******************************** datos <- dat # SPECIFY DATA SET #******************************** #run loop for (i in 1:50){ AinvTA.2<-inverseA(tree[[i]])$Ainv m1 <- MCMCglmm(cbind(zr_f, zr_m) ~ trait + overall_parameter + trait:overall_parameter + scale(year, scale= FALSE) -1, random = ~ us(trait):animal + us(trait):scientific_names, rcov = ~ us(trait):units, mev = c(datos$v_fem, datos$v_male), ginverse=list(animal = AinvTA.2), data = datos, family = rep('gaussian', 2), nitt=11000, thin=100, burnin=1000, pl = TRUE, pr = TRUE, prior = prior, verbose=FALSE) modelos[[i]] <- m1 print(i) } summary(modelos[[1]]) #eg output of the first model #processing multimodel results---- #combining posterior distributions across models all.m.vcv <- modelos[[1]]$VCV all.m.sol <- modelos[[1]]$Sol for (i in 2:50){ all.m.vcv <- rbind(all.m.vcv, modelos[[i]]$VCV) all.m.sol <- rbind(all.m.sol, modelos[[i]]$Sol) } #converting into mcmc objects all.m.vcv <- as.mcmc(all.m.vcv) all.m.sol <- as.mcmc(all.m.sol) #combine output into a list---- #this includes: #(a) a list with all 50 models #(b) the combined posterior distribution of all sol in (a) #(c) the combined posterior distribution of all vcv in (a) #*************************************************************************************** output.model2 <- (list(models=modelos, sol=all.m.sol, vcv=all.m.vcv)) # SPECIFY OUTPUT NAME #*************************************************************************************** #remove temporary files rm(modelos, all.m.sol, all.m.vcv) #summarising results---- #*************************************************************************************** mod <- output.model2 # SPECIFY NAME OF MODEL THAT WE WANT TO SUMMARISE #*************************************************************************************** #fixed effects #posterior means colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) #95%CIs HPDinterval(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) #random effects #posterior means colMeans(mod$vcv) #95%CIs HPDinterval(mod$vcv) #computing effective sample sizes effectiveSize(mod$vcv) effectiveSize(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) results.table <- rbind( #sol results data.frame(parameter=names(colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl])), post.mean=colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl]), HPDinterval(mod$sol[,1:mod$models[[1]]$Fixed$nfl]), eff.sample.size=effectiveSize(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) ), #vcv results data.frame(parameter=names(colMeans(mod$vcv)), post.mean=colMeans(mod$vcv), HPDinterval(mod$vcv), eff.sample.size=effectiveSize(mod$vcv)) ) results.table <- subset(results.table, parameter!="sqrt(mev):sqrt(mev).meta") # eliminates row for measurement error rownames(results.table) <- NULL #********************************************************************** table.2 <- results.table #********************************************************************** # 95% CIs parameter calculations per sex and testing for differences between sexes #Differences between sexes #Condition colMeans(mod$sol[,1:2]) HPDinterval(mod$sol[,1:2]) #95% CIs per sex sex.diff <- mod$sol[,1]-mod$sol[,2] #fem - male mean(sex.diff) HPDinterval(sex.diff) #Fitness fem.fit <- mod$sol[,1]+mod$sol[,3] #posterior distribution for female fitness male.fit <- mod$sol[,2]+mod$sol[,3]+mod$sol[,5] #posterior distribution for male fitness colMeans(cbind(fem.fit, male.fit)) HPDinterval(fem.fit) #95% CIs female HPDinterval(male.fit) #95% CIs male sex.diff <- fem.fit - male.fit mean(sex.diff) HPDinterval(sex.diff) #Differences between condition and fitness associations per sex fem.diff <- mod$sol[,1] - fem.fit #for females condition - fitness mean(fem.diff) HPDinterval(fem.diff) male.diff <- mod$sol[,2] - male.fit #for males condition - fitness mean(male.diff) HPDinterval(male.diff) #Proportion of the variance explained by phylogeny and species id # phylogeny prop.phyl.var.fem <- mod$vcv[,1]/ (mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10]) #females mean(prop.phyl.var.fem) HPDinterval(prop.phyl.var.fem) prop.phyl.var.male <- mod$vcv[,4]/ (mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13]) #males mean(prop.phyl.var.male) HPDinterval(prop.phyl.var.male) # species id prop.spp.var.fem <- mod$vcv[,5]/ (mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10]) #females mean(prop.spp.var.fem) HPDinterval(prop.spp.var.fem) prop.spp.var.male <- mod$vcv[,8]/ (mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13]) #males mean(prop.spp.var.male) HPDinterval(prop.spp.var.male) #residual prop.residual.fem <- mod$vcv[,10]/(mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10]) #females mean(prop.residual.fem) HPDinterval(prop.residual.fem) prop.residual.male <- mod$vcv[,13]/(mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13]) #males mean(prop.residual.male) HPDinterval(prop.residual.male) #correlations for phylogeny # getting correlations - (correlation = cov/(SD*SD)) spp_mat <- mod$vcv[,c(1, 3, 4)] attr(spp_mat, "dimnames")[[2]]<- c("f_V", "fm_cov", "m_V") # female variance, covariance, male variance spp_cor <- spp_mat[, "fm_cov"]/(sqrt(spp_mat[, "f_V"]*spp_mat[, "m_V"])) phylo.cor <- spp_cor mean(phylo.cor) HPDinterval(phylo.cor) summary(phylo.cor) #slope for phylogeny #computed as covariance divided by male variance slope_phylo <- spp_mat[,2]/spp_mat[,3] summary(slope_phylo) #correlations for species id # getting correlations - spp (correlation = cov/(SD*SD)) spp_mat <- mod$vcv[,c(5, 7, 8)] attr(spp_mat, "dimnames")[[2]]<- c("f_V", "fm_cov", "m_V") # female, male variance and cov spp_cor <- spp_mat[, "fm_cov"]/(sqrt(spp_mat[, "f_V"]*spp_mat[, "m_V"])) sp.id.cor <- spp_cor mean(sp.id.cor) HPDinterval(sp.id.cor) summary(sp.id.cor) #slope for species #computed as covariance divided by male variance slope_species <- spp_mat[,2]/spp_mat[,3] summary(slope_species) #correlations for residual # getting correlations - spp (correlation = cov/(SD*SD)) spp_mat <- mod$vcv[,c(10, 12, 13)] attr(spp_mat, "dimnames")[[2]]<- c("f_V", "fm_cov", "m_V") # female, male variance and cov spp_cor <- spp_mat[, "fm_cov"]/(sqrt(spp_mat[, "f_V"]*spp_mat[, "m_V"])) res.cor <- spp_cor mean(res.cor) HPDinterval(res.cor) summary(res.cor) #slope for units #computed as covariance divided by male variance slope_units <- spp_mat[,2]/spp_mat[,3] summary(slope_units) #Publication bias---- # predictions for fixed effects pred_fix <- as.vector(apply(mod$sol,2,mean)[1:mod$models[[1]]$Fixed$nfl] %*% t(mod$models[[1]]$X)) # predictions for random effects (including mev) pred_ran <- as.vector(apply(mod$sol,2,mean)[-(1:mod$models[[1]]$Fixed$nfl)] %*% t(mod$models[[1]]$Z)) # predictions for sampling variance raw_se <- c(datos$zr_se_fem, datos$zr_se_male) # observed data raw_data <- c(datos$zr_f, datos$zr_m) # residuals - the first 737 for females and the next 737 for males pred_res <- raw_data - (pred_fix + pred_ran) # precision precision <- c(1/datos$zr_se_fem, 1/datos$zr_se_male) #funnel plot #females length(datos$overall_parameter) #getting the number of observations for the subsetted data plot(pred_res[1:737], precision[1:737]) #males plot(pred_res[738:1476], precision[738:1476]) #female residuals residuales <- pred_res[1:737] precis <- precision[1:737] #Trim and fill procedures metaR <- rma(yi= residuales, sei = 1/precis) TF <-trimfill(metaR,estimator="R0") TF funnel(TF, 'seinv') #Egger's test regtest(metaR, model="lm", predictor="sei") #male residuals residuales <- pred_res[738:1476] precis <- precision[738:1476] #Trim and fill procedures metaR <- rma(yi= residuales, sei = 1/precis) TF <-trimfill(metaR,estimator="R0") TF funnel(TF, 'seinv') regtest(metaR, model="lm", predictor="sei") #MODEL 3 - MODEL FOR ONLY CONDITION SPECIFIC PARAMETERS (INTERACTIONS WITH SEX) ---- # prior prior <- list(R = list(V = diag(2), nu = 2), G = list(G1 = list(V = diag(2) , nu = 2, alpha.mu=c(0,0), alpha.V=diag(2)*25^2), G2 = list(V = diag(2) , nu = 2, alpha.mu=c(0,0), alpha.V=diag(2)*25^2))) modelos <- NULL # setting up an empty list to store models class(modelos) <- 'list' #******************************** datos <- subset(dat, overall_parameter=="Body Condition") # SPECIFY DATA SET #******************************** #run loop for (i in 1:50){ AinvTA.2<-inverseA(tree[[i]])$Ainv m1 <- MCMCglmm(cbind(zr_f, zr_m) ~ trait + specific_parameter + trait:specific_parameter + scale(year, scale= FALSE) -1, random = ~ us(trait):animal + us(trait):scientific_names, rcov = ~ us(trait):units, mev = c(datos$v_fem, datos$v_male), ginverse=list(animal = AinvTA.2), data = datos, family = rep('gaussian', 2), nitt=11000, thin=100, burnin=1000, pl = TRUE, pr = TRUE, prior = prior, verbose=FALSE) modelos[[i]] <- m1 print(i) } summary(modelos[[1]]) #eg output of the first model #processing multimodel results---- #combining posterior distributions across models all.m.vcv <- modelos[[1]]$VCV all.m.sol <- modelos[[1]]$Sol for (i in 2:50){ all.m.vcv <- rbind(all.m.vcv, modelos[[i]]$VCV) all.m.sol <- rbind(all.m.sol, modelos[[i]]$Sol) } #converting into mcmc objects all.m.vcv <- as.mcmc(all.m.vcv) all.m.sol <- as.mcmc(all.m.sol) #combine output into a list---- #this includes: #(a) a list with all 50 models #(b) the combined posterior distribution of all sol in (a) #(c) the combined posterior distribution of all vcv in (a) #*************************************************************************************** output.model3 <- (list(models=modelos, sol=all.m.sol, vcv=all.m.vcv)) # SPECIFY OUTPUT NAME #*************************************************************************************** #remove temporary files rm(modelos, all.m.sol, all.m.vcv) #summarising results---- #*************************************************************************************** mod <- output.model3 # SPECIFY NAME OF MODEL THAT WE WANT TO SUMMARISE #*************************************************************************************** #fixed effects #posterior means colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) #95%CIs HPDinterval(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) #random effects #posterior means colMeans(mod$vcv) #95%CIs HPDinterval(mod$vcv) #computing effective sample sizes effectiveSize(mod$vcv) effectiveSize(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) results.table <- rbind( #sol results data.frame(parameter=names(colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl])), post.mean=colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl]), HPDinterval(mod$sol[,1:mod$models[[1]]$Fixed$nfl]), eff.sample.size=effectiveSize(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) ), #vcv results data.frame(parameter=names(colMeans(mod$vcv)), post.mean=colMeans(mod$vcv), HPDinterval(mod$vcv), eff.sample.size=effectiveSize(mod$vcv)) ) results.table <- subset(results.table, parameter!="sqrt(mev):sqrt(mev).meta") # eliminates row for measurement error rownames(results.table) <- NULL #********************************************************************** table.3 <- results.table # give results.table a name #********************************************************************** # Differences between sexes (95% CIs parameter calculations per sex) #body condition fem.bc <- mod$sol[,1] #posterior distribution for female body condition male.bc <- mod$sol[,2] #posterior distribution for male body condition colMeans(cbind(fem.bc, male.bc)) HPDinterval(fem.bc) #95% CIs female HPDinterval(male.bc) #95% CIs male sex.diff <- fem.bc - male.bc HPDinterval(sex.diff) #body size fem.size <- mod$sol[,1]+mod$sol[,3] #posterior distribution for female body size male.size <- mod$sol[,2]+mod$sol[,3]+mod$sol[,9] #posterior distribution for male body size colMeans(cbind(fem.size, male.size)) HPDinterval(fem.size) #95% CIs female HPDinterval(male.size) #95% CIs male sex.diff <- fem.size - male.size HPDinterval(sex.diff) #immunity fem.imm <- mod$sol[,1]+mod$sol[,5] #posterior distribution for female immunity male.imm <- mod$sol[,2]+mod$sol[,5]+mod$sol[,11] #posterior distribution for male immunity colMeans(cbind(fem.imm, male.imm)) HPDinterval(fem.imm) #95% CIs female HPDinterval(male.imm) #95% CIs male sex.diff <- fem.imm - male.imm HPDinterval(sex.diff) #stress fem.str <- mod$sol[,1]+mod$sol[,7] #posterior distribution for female immunity male.str <- mod$sol[,2]+mod$sol[,7]+mod$sol[,13] #posterior distribution for male immunity colMeans(cbind(fem.str, male.str)) HPDinterval(fem.str) #95% CIs female HPDinterval(male.str) #95% CIs male sex.diff <- fem.str - male.str HPDinterval(sex.diff) #environment fem.env <- mod$sol[,1]+mod$sol[,4] #posterior distribution for female environment male.env <- mod$sol[,2]+mod$sol[,4]+mod$sol[,10] #posterior distribution for male environment colMeans(cbind(fem.env, male.env)) HPDinterval(fem.env) #95% CIs female HPDinterval(male.env) #95% CIs male sex.diff <- fem.env - male.env HPDinterval(sex.diff) #parasites fem.par <- mod$sol[,1]+mod$sol[,6] #posterior distribution for female immunity male.par <- mod$sol[,2]+mod$sol[,6]+mod$sol[,12] #posterior distribution for male immunity colMeans(cbind(fem.par, male.par)) HPDinterval(fem.par) #95% CIs female HPDinterval(male.par) #95% CIs male sex.diff <- fem.par - male.par HPDinterval(sex.diff) #Proportion of the variance explained by phylogeny and species id # phylogeny colnames(mod$vcv) prop.phyl.var.fem <- mod$vcv[,1]/(mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10]) #females mean(prop.phyl.var.fem) HPDinterval(prop.phyl.var.fem) prop.phyl.var.male <- mod$vcv[,4]/ (mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13]) #males mean(prop.phyl.var.male) HPDinterval(prop.phyl.var.male) # species id prop.spp.var.fem <- mod$vcv[,5]/ (mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10]) #females mean(prop.spp.var.fem) HPDinterval(prop.spp.var.fem) prop.spp.var.male <- mod$vcv[,8]/ (mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13]) #males mean(prop.spp.var.male) HPDinterval(prop.spp.var.male) #residual prop.residual.fem <- mod$vcv[,10]/(mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10]) #females mean(prop.residual.fem) HPDinterval(prop.residual.fem) prop.residual.male <- mod$vcv[,13]/(mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13]) #males mean(prop.residual.male) HPDinterval(prop.residual.male) #correlations for phylogeny # getting correlations - (correlation = cov/(SD*SD)) spp_mat <- mod$vcv[,c(1, 3, 4)] attr(spp_mat, "dimnames")[[2]]<- c("f_V", "fm_cov", "m_V") # female variance, covariance, male variance spp_cor <- spp_mat[, "fm_cov"]/(sqrt(spp_mat[, "f_V"]*spp_mat[, "m_V"])) phylo.cor <- spp_cor mean(phylo.cor) HPDinterval(phylo.cor) summary(phylo.cor) #slope for phylogeny #computed as covariance divided by male variance slope_phylo <- spp_mat[,2]/spp_mat[,3] summary(slope_phylo) #correlations for species id # getting correlations - spp (correlation = cov/(SD*SD)) spp_mat <- mod$vcv[,c(5, 7, 8)] attr(spp_mat, "dimnames")[[2]]<- c("f_V", "fm_cov", "m_V") # female, male variance and cov spp_cor <- spp_mat[, "fm_cov"]/(sqrt(spp_mat[, "f_V"]*spp_mat[, "m_V"])) sp.id.cor <- spp_cor mean(sp.id.cor) HPDinterval(sp.id.cor) summary(sp.id.cor) #slope for species #computed as covariance divided by male variance slope_species <- spp_mat[,2]/spp_mat[,3] summary(slope_species) #correlations for residual # getting correlations - spp (correlation = cov/(SD*SD)) spp_mat <- mod$vcv[,c(10, 12, 13)] attr(spp_mat, "dimnames")[[2]]<- c("f_V", "fm_cov", "m_V") # female, male variance and cov spp_cor <- spp_mat[, "fm_cov"]/(sqrt(spp_mat[, "f_V"]*spp_mat[, "m_V"])) res.cor <- spp_cor mean(res.cor) HPDinterval(res.cor) summary(res.cor) #slope for units #computed as covariance divided by male variance slope_units <- spp_mat[,2]/spp_mat[,3] summary(slope_units) #Publication bias---- # predictions for fixed effects pred_fix <- as.vector(apply(mod$sol,2,mean)[1:mod$models[[1]]$Fixed$nfl] %*% t(mod$models[[1]]$X)) # predictions for random effects (including mev) pred_ran <- as.vector(apply(mod$sol,2,mean)[-(1:mod$models[[1]]$Fixed$nfl)] %*% t(mod$models[[1]]$Z)) # predictions for sampling variance raw_se <- c(datos$zr_se_fem, datos$zr_se_male) # observed data raw_data <- c(datos$zr_f, datos$zr_m) # residuals pred_res <- raw_data - (pred_fix + pred_ran) # precision precision <- c(1/datos$zr_se_fem, 1/datos$zr_se_male) #funnel plots- the first 304 for females and the rest for males length(datos$overall_parameter) #females plot(pred_res[1:304], precision[1:304]) #males plot(pred_res[305:608], precision[305:608]) #female residuals residuales <- pred_res[1:304] precis <- precision[1:304] #Trim and fill procedures metaR <- rma(yi= residuales, sei = 1/precis) TF <-trimfill(metaR,estimator="R0") TF funnel(TF, 'seinv') #Egger's test regtest(metaR, model="lm", predictor="sei") #male residuals residuales <- pred_res[305:608] precis <- precision[305:608] #Trim and fill procedures metaR <- rma(yi= residuales, sei = 1/precis) TF <-trimfill(metaR,estimator="R0") TF funnel(TF, 'seinv') regtest(metaR, model="lm", predictor="sei") #MODEL 4 - MODEL FOR ONLY FITNESS SPECIFIC PARAMETERS (INTERACTIONS WITH SEX) ---- # prior prior <- list(R = list(V = diag(2), nu = 2), G = list(G1 = list(V = diag(2) , nu = 2, alpha.mu=c(0,0), alpha.V=diag(2)*25^2), G2 = list(V = diag(2) , nu = 2, alpha.mu=c(0,0), alpha.V=diag(2)*25^2))) modelos <- NULL # setting up an empty list to store models class(modelos) <- 'list' #******************************** datos <- subset(dat, overall_parameter=="Fitness") # SPECIFY DATA SET #******************************** #run loop for (i in 1:50){ AinvTA.2<-inverseA(tree[[i]])$Ainv m1 <- MCMCglmm(cbind(zr_f, zr_m) ~ trait + specific_parameter + trait:specific_parameter + scale(year, scale= FALSE) -1, random = ~ us(trait):animal + us(trait):scientific_names, rcov = ~ us(trait):units, mev = c(datos$v_fem, datos$v_male), ginverse=list(animal = AinvTA.2), data = datos, family = rep('gaussian', 2), nitt=11000, thin=100, burnin=1000, pl = TRUE, pr = TRUE, prior = prior, verbose=FALSE) modelos[[i]] <- m1 print(i) } summary(modelos[[1]]) #eg output of the first model #processing multimodel results---- #combining posterior distributions across models all.m.vcv <- modelos[[1]]$VCV all.m.sol <- modelos[[1]]$Sol for (i in 2:50){ all.m.vcv <- rbind(all.m.vcv, modelos[[i]]$VCV) all.m.sol <- rbind(all.m.sol, modelos[[i]]$Sol) } #converting into mcmc objects all.m.vcv <- as.mcmc(all.m.vcv) all.m.sol <- as.mcmc(all.m.sol) #combine output into a list---- #this includes: #(a) a list with all 50 models #(b) the combined posterior distribution of all sol in (a) #(c) the combined posterior distribution of all vcv in (a) #*************************************************************************************** output.model4 <- (list(models=modelos, sol=all.m.sol, vcv=all.m.vcv)) # SPECIFY OUTPUT NAME #*************************************************************************************** #remove temporary files rm(modelos, all.m.sol, all.m.vcv) #summarising results---- #*************************************************************************************** mod <- output.model4 # SPECIFY NAME OF MODEL THAT WE WANT TO SUMMARISE #*************************************************************************************** #fixed effects #posterior means colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) #95%CIs HPDinterval(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) #random effects #posterior means colMeans(mod$vcv) #95%CIs HPDinterval(mod$vcv) #computing effective sample sizes effectiveSize(mod$vcv) effectiveSize(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) results.table <- rbind( #sol results data.frame(parameter=names(colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl])), post.mean=colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl]), HPDinterval(mod$sol[,1:mod$models[[1]]$Fixed$nfl]), eff.sample.size=effectiveSize(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) ), #vcv results data.frame(parameter=names(colMeans(mod$vcv)), post.mean=colMeans(mod$vcv), HPDinterval(mod$vcv), eff.sample.size=effectiveSize(mod$vcv)) ) results.table <- subset(results.table, parameter!="sqrt(mev):sqrt(mev).meta") # eliminates row for measurement error rownames(results.table) <- NULL #********************************************************************** table.4 <- results.table # give results.table a name #********************************************************************** # Differences between sexes (95% CIs parameter calculations per sex) #reproductive success fem.rep <- mod$sol[,1]+mod$sol[,5] #posterior distribution for female reproductive success male.rep <- mod$sol[,2]+mod$sol[,5]+mod$sol[,10]#posterior distribution for male reproductive success colMeans(cbind(fem.rep, male.rep)) HPDinterval(fem.rep) #95% CIs female HPDinterval(male.rep) #95% CIs male sex.diff <- fem.rep - male.rep mean(sex.diff) HPDinterval(sex.diff) #offspring quality or condition fem.off <- mod$sol[,1]+mod$sol[,3]#posterior distribution for female offspring quality male.off <- mod$sol[,2]+mod$sol[,3]+mod$sol[,8] #posterior distribution for male offspring quality colMeans(cbind(fem.off, male.off)) HPDinterval(fem.off) #95% CIs female HPDinterval(male.off) #95% CIs male sex.diff <- fem.off - male.off HPDinterval(sex.diff) #parental quality fem.par <- mod$sol[,1]+mod$sol[,4] #posterior distribution for female parental quality male.par <- mod$sol[,2]+mod$sol[,4]+mod$sol[,9] #posterior distribution for male parental quality colMeans(cbind(fem.par, male.par)) HPDinterval(fem.par) #95% CIs female HPDinterval(male.par) #95% CIs male sex.diff <- fem.par - male.par HPDinterval(sex.diff) #timing of breeding fem.tim <- mod$sol[,1] #posterior distribution for female timing of breeding male.tim <- mod$sol[,2] #posterior distribution for male timing of breeding colMeans(cbind(fem.tim, male.tim)) HPDinterval(fem.tim) #95% CIs female HPDinterval(male.tim) #95% CIs male sex.diff <- fem.tim - male.tim HPDinterval(sex.diff) #survival fem.sur <- mod$sol[,1]+mod$sol[,6] #posterior distribution for female survival male.sur <- mod$sol[,2]+mod$sol[,6]+mod$sol[,11] #posterior distribution for male survival colMeans(cbind(fem.sur, male.sur)) HPDinterval(fem.sur) #95% CIs female HPDinterval(male.sur) #95% CIs male sex.diff <- fem.sur - male.sur HPDinterval(sex.diff) #Proportion of the variance explained by phylogeny and species id # phylogeny prop.phyl.var.fem <- mod$vcv[,1]/(mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10]) #females mean(prop.phyl.var.fem) HPDinterval(prop.phyl.var.fem) prop.phyl.var.male <- mod$vcv[,4]/ (mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13]) #males mean(prop.phyl.var.male) HPDinterval(prop.phyl.var.male) # species id prop.spp.var.fem <- mod$vcv[,5]/ (mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10]) #females mean(prop.spp.var.fem) HPDinterval(prop.spp.var.fem) prop.spp.var.male <- mod$vcv[,8]/ (mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13]) #males mean(prop.spp.var.male) HPDinterval(prop.spp.var.male) #residual prop.residual.fem <- mod$vcv[,10]/(mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10]) #females mean(prop.residual.fem) HPDinterval(prop.residual.fem) prop.residual.male <- mod$vcv[,13]/(mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13]) #males mean(prop.residual.male) HPDinterval(prop.residual.male) #correlations for phylogeny # getting correlations - (correlation = cov/(SD*SD)) spp_mat <- mod$vcv[,c(1, 3, 4)] attr(spp_mat, "dimnames")[[2]]<- c("f_V", "fm_cov", "m_V") # female variance, covariance, male variance spp_cor <- spp_mat[, "fm_cov"]/(sqrt(spp_mat[, "f_V"]*spp_mat[, "m_V"])) phylo.cor <- spp_cor mean(phylo.cor) HPDinterval(phylo.cor) summary(phylo.cor) #slope for phylogeny #computed as covariance divided by male variance slope_phylo <- spp_mat[,2]/spp_mat[,3] summary(slope_phylo) #correlations for species id # getting correlations - spp (correlation = cov/(SD*SD)) spp_mat <- mod$vcv[,c(5, 7, 8)] attr(spp_mat, "dimnames")[[2]]<- c("f_V", "fm_cov", "m_V") # female, male variance and cov spp_cor <- spp_mat[, "fm_cov"]/(sqrt(spp_mat[, "f_V"]*spp_mat[, "m_V"])) sp.id.cor <- spp_cor mean(sp.id.cor) HPDinterval(sp.id.cor) summary(sp.id.cor) #slope for species #computed as covariance divided by male variance slope_species <- spp_mat[,2]/spp_mat[,3] summary(slope_species) #correlations for residual # getting correlations - spp (correlation = cov/(SD*SD)) spp_mat <- mod$vcv[,c(10, 12, 13)] attr(spp_mat, "dimnames")[[2]]<- c("f_V", "fm_cov", "m_V") # female, male variance and cov spp_cor <- spp_mat[, "fm_cov"]/(sqrt(spp_mat[, "f_V"]*spp_mat[, "m_V"])) res.cor <- spp_cor mean(res.cor) HPDinterval(res.cor) summary(res.cor) #slope for units #computed as covariance divided by male variance slope_units <- spp_mat[,2]/spp_mat[,3] summary(slope_units) #Publication bias---- # predictions for fixed effects pred_fix <- as.vector(apply(mod$sol,2,mean)[1:mod$models[[1]]$Fixed$nfl] %*% t(mod$models[[1]]$X)) # predictions for random effects (including mev) pred_ran <- as.vector(apply(mod$sol,2,mean)[-(1:mod$models[[1]]$Fixed$nfl)] %*% t(mod$models[[1]]$Z)) # predictions for sampling variance raw_se <- c(datos$zr_se_fem, datos$zr_se_male) #pred_mev <- apply(mod$sol,2,mean)[(ncol(mod$sol)-nrow(datos)*2 + 1):ncol(mod$sol)]*raw_se # observed data raw_data <- c(datos$zr_f, datos$zr_m) # residuals - the first 737 for females and the next 737 for males pred_res <- raw_data - (pred_fix + pred_ran) # precision precision <- c(1/datos$zr_se_fem, 1/datos$zr_se_male) #funnel plots #females length(datos$overall_parameter) #getting the number of observations for the subsetted data plot(pred_res[1:433], precision[1:433]) #males plot(pred_res[434:866], precision[434:866]) #female residuals residuales <- pred_res[1:433] precis <- precision[1:433] #Trim and fill procedures metaR <- rma(yi= residuales, sei = 1/precis) TF <-trimfill(metaR,estimator="R0") TF funnel(TF, 'seinv') #Egger's test regtest(metaR, model="lm", predictor="sei") #male residuals residuales <- pred_res[434:866] precis <- precision[434:866] #Trim and fill procedures metaR <- rma(yi= residuales, sei = 1/precis) TF <-trimfill(metaR,estimator="R0") TF funnel(TF, 'seinv') regtest(metaR, model="lm", predictor="sei") #MODEL 5 - MODEL INCORPORATING SEXUAL DIMORPHISM FOR CONDITION DATA ONLY (Excluding right outliers)---- # prior prior <- list(R = list(V = diag(2), nu = 2), G = list(G1 = list(V = diag(2) , nu = 2, alpha.mu=c(0,0), alpha.V=diag(2)*25^2), G2 = list(V = diag(2) , nu = 2, alpha.mu=c(0,0), alpha.V=diag(2)*25^2))) modelos <- NULL # setting up an empty list to store models class(modelos) <- 'list' #******************************** datos <- subset(dat, dimorphism_Cohend!="NA" & dimorphism_Cohend < 11 & overall_parameter=="Body Condition") # SPECIFY DATA SET (no missing values) #******************************** #run loop. for (i in 1:50){ AinvTA.2<-inverseA(tree[[i]])$Ainv m1 <- MCMCglmm(cbind(zr_f, zr_m) ~ trait + dimorphism_Cohend + trait:dimorphism_Cohend + scale(year, scale= FALSE) -1, random = ~ us(trait):animal + us(trait):scientific_names, rcov = ~ us(trait):units, mev = c(datos$v_fem, datos$v_male), ginverse=list(animal = AinvTA.2), data = datos, family = rep('gaussian', 2), nitt=11000, thin=100, burnin=1000, pl = TRUE, pr = TRUE, prior = prior, verbose=FALSE) modelos[[i]] <- m1 print(i) } summary(modelos[[1]]) #eg output of the first model #processing multimodel results---- #combining posterior distributions across models all.m.vcv <- modelos[[1]]$VCV all.m.sol <- modelos[[1]]$Sol for (i in 2:50){ all.m.vcv <- rbind(all.m.vcv, modelos[[i]]$VCV) all.m.sol <- rbind(all.m.sol, modelos[[i]]$Sol) } #converting into mcmc objects all.m.vcv <- as.mcmc(all.m.vcv) all.m.sol <- as.mcmc(all.m.sol) #combine output into a list---- #this includes: #(a) a list with all 50 models #(b) the combined posterior distribution of all sol in (a) #(c) the combined posterior distribution of all vcv in (a) #*************************************************************************************** output.model5 <- (list(models=modelos, sol=all.m.sol, vcv=all.m.vcv)) # SPECIFY OUTPUT NAME #*************************************************************************************** #remove temporary files rm(modelos, all.m.sol, all.m.vcv) #summarising results---- #*************************************************************************************** mod <- output.model5 # SPECIFY NAME OF MODEL THAT WE WANT TO SUMMARISE #*************************************************************************************** #fixed effects #posterior means colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) #95%CIs HPDinterval(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) #random effects #posterior means colMeans(mod$vcv) #95%CIs HPDinterval(mod$vcv) #computing effective sample sizes effectiveSize(mod$vcv) effectiveSize(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) results.table <- rbind( #sol results data.frame(parameter=names(colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl])), post.mean=colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl]), HPDinterval(mod$sol[,1:mod$models[[1]]$Fixed$nfl]), eff.sample.size=effectiveSize(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) ), #vcv results data.frame(parameter=names(colMeans(mod$vcv)), post.mean=colMeans(mod$vcv), HPDinterval(mod$vcv), eff.sample.size=effectiveSize(mod$vcv)) ) results.table <- subset(results.table, parameter!="sqrt(mev):sqrt(mev).meta") # eliminates row for measurement error rownames(results.table) <- NULL #********************************************************************** table.5 <- results.table # give results.table a name #********************************************************************** #Proportion of the variance explained by phylogeny and species id # phylogeny prop.phyl.var.fem <- mod$vcv[,1]/(mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10]) #females mean(prop.phyl.var.fem) HPDinterval(prop.phyl.var.fem) prop.phyl.var.male <- mod$vcv[,4]/ (mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13]) #males mean(prop.phyl.var.male) HPDinterval(prop.phyl.var.male) # species id prop.spp.var.fem <- mod$vcv[,5]/ (mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10]) #females mean(prop.spp.var.fem) HPDinterval(prop.spp.var.fem) prop.spp.var.male <- mod$vcv[,8]/ (mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13]) #males mean(prop.spp.var.male) HPDinterval(prop.spp.var.male) #residual prop.residual.fem <- mod$vcv[,10]/(mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10]) #females mean(prop.residual.fem) HPDinterval(prop.residual.fem) prop.residual.male <- mod$vcv[,13]/(mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13]) #males mean(prop.residual.male) HPDinterval(prop.residual.male) #correlations for phylogeny # getting correlations - (correlation = cov/(SD*SD)) spp_mat <- mod$vcv[,c(1, 3, 4)] attr(spp_mat, "dimnames")[[2]]<- c("f_V", "fm_cov", "m_V") # female variance, covariance, male variance spp_cor <- spp_mat[, "fm_cov"]/(sqrt(spp_mat[, "f_V"]*spp_mat[, "m_V"])) phylo.cor <- spp_cor mean(phylo.cor) HPDinterval(phylo.cor) summary(phylo.cor) #slope for phylogeny #computed as covariance divided by male variance slope_phylo <- spp_mat[,2]/spp_mat[,3] summary(slope_phylo) #correlations for species id # getting correlations - spp (correlation = cov/(SD*SD)) spp_mat <- mod$vcv[,c(5, 7, 8)] attr(spp_mat, "dimnames")[[2]]<- c("f_V", "fm_cov", "m_V") # female, male variance and cov spp_cor <- spp_mat[, "fm_cov"]/(sqrt(spp_mat[, "f_V"]*spp_mat[, "m_V"])) sp.id.cor <- spp_cor mean(sp.id.cor) HPDinterval(sp.id.cor) summary(sp.id.cor) #slope for species #computed as covariance divided by male variance slope_species <- spp_mat[,2]/spp_mat[,3] summary(slope_species) #correlations for residual # getting correlations - spp (correlation = cov/(SD*SD)) spp_mat <- mod$vcv[,c(10, 12, 13)] attr(spp_mat, "dimnames")[[2]]<- c("f_V", "fm_cov", "m_V") # female, male variance and cov spp_cor <- spp_mat[, "fm_cov"]/(sqrt(spp_mat[, "f_V"]*spp_mat[, "m_V"])) res.cor <- spp_cor mean(res.cor) HPDinterval(res.cor) summary(res.cor) #slope for units #computed as covariance divided by male variance slope_units <- spp_mat[,2]/spp_mat[,3] summary(slope_units) #Publication bias---- # predictions for fixed effects pred_fix <- as.vector(apply(mod$sol,2,mean)[1:mod$models[[1]]$Fixed$nfl] %*% t(mod$models[[1]]$X)) # predictions for random effects (including mev) pred_ran <- as.vector(apply(mod$sol,2,mean)[-(1:mod$models[[1]]$Fixed$nfl)] %*% t(mod$models[[1]]$Z)) # predictions for sampling variance raw_se <- c(datos$zr_se_fem, datos$zr_se_male) #pred_mev <- apply(mod$sol,2,mean)[(ncol(mod$sol)-nrow(datos)*2 + 1):ncol(mod$sol)]*raw_se # observed data raw_data <- c(datos$zr_f, datos$zr_m) # residuals - the first 737 for females and the next 737 for males pred_res <- raw_data - (pred_fix + pred_ran) # precision precision <- c(1/datos$zr_se_fem, 1/datos$zr_se_male) #funnel plots #females length(datos$overall_parameter) #getting the numbers of observations for the subsetted data plot(pred_res[1:105], precision[1:105]) #males plot(pred_res[106:210], precision[106:210]) #female residuals residuales <- pred_res[1:105] precis <- precision[1:105] #Trim and fill procedures metaR <- rma(yi= residuales, sei = 1/precis) TF <-trimfill(metaR,estimator="R0") TF funnel(TF, 'seinv') #Egger's test regtest(metaR, model="lm", predictor="sei") #male residuals residuales <- pred_res[106:210] precis <- precision[106:210] #Trim and fill procedures metaR <- rma(yi= residuales, sei = 1/precis) TF <-trimfill(metaR,estimator="R0") TF funnel(TF, 'seinv') regtest(metaR, model="lm", predictor="sei") #MODEL 6 - MODEL INCORPORATING SEXUAL DIMORPHISM FOR FITNESS DATA ONLY---- # prior prior <- list(R = list(V = diag(2), nu = 2), G = list(G1 = list(V = diag(2) , nu = 2, alpha.mu=c(0,0), alpha.V=diag(2)*25^2), G2 = list(V = diag(2) , nu = 2, alpha.mu=c(0,0), alpha.V=diag(2)*25^2))) modelos <- NULL # setting up an empty list to store models class(modelos) <- 'list' #******************************** datos <- subset(dat, dimorphism_Cohend!="NA" & overall_parameter=="Fitness") # SPECIFY DATA SET (no missing values) #******************************** #run loop. sex*dimorphism model for (i in 1:50){ AinvTA.2<-inverseA(tree[[i]])$Ainv m1 <- MCMCglmm(cbind(zr_f, zr_m) ~ trait + dimorphism_Cohend + trait:dimorphism_Cohend + scale(year, scale= FALSE) -1, random = ~ us(trait):animal + us(trait):scientific_names, rcov = ~ us(trait):units, mev = c(datos$v_fem, datos$v_male), ginverse=list(animal = AinvTA.2), data = datos, family = rep('gaussian', 2), nitt=11000, thin=100, burnin=1000, pl = TRUE, pr = TRUE, prior = prior, verbose=FALSE) modelos[[i]] <- m1 print(i) } summary(modelos[[1]]) #eg output of the first model #processing multimodel results---- #combining posterior distributions across models all.m.vcv <- modelos[[1]]$VCV all.m.sol <- modelos[[1]]$Sol for (i in 2:50){ all.m.vcv <- rbind(all.m.vcv, modelos[[i]]$VCV) all.m.sol <- rbind(all.m.sol, modelos[[i]]$Sol) } #converting into mcmc objects all.m.vcv <- as.mcmc(all.m.vcv) all.m.sol <- as.mcmc(all.m.sol) #combine output into a list---- #this includes: #(a) a list with all 50 models #(b) the combined posterior distribution of all sol in (a) #(c) the combined posterior distribution of all vcv in (a) #*************************************************************************************** output.model6 <- (list(models=modelos, sol=all.m.sol, vcv=all.m.vcv)) # SPECIFY OUTPUT NAME #*************************************************************************************** #remove temporary files rm(modelos, all.m.sol, all.m.vcv) #summarising results---- #*************************************************************************************** mod <- output.model6 # SPECIFY NAME OF MODEL THAT WE WANT TO SUMMARISE #*************************************************************************************** #fixed effects #posterior means colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) #95%CIs HPDinterval(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) #random effects #posterior means colMeans(mod$vcv) #95%CIs HPDinterval(mod$vcv) #computing effective sample sizes effectiveSize(mod$vcv) effectiveSize(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) results.table <- rbind( #sol results data.frame(parameter=names(colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl])), post.mean=colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl]), HPDinterval(mod$sol[,1:mod$models[[1]]$Fixed$nfl]), eff.sample.size=effectiveSize(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) ), #vcv results data.frame(parameter=names(colMeans(mod$vcv)), post.mean=colMeans(mod$vcv), HPDinterval(mod$vcv), eff.sample.size=effectiveSize(mod$vcv)) ) results.table <- subset(results.table, parameter!="sqrt(mev):sqrt(mev).meta") # eliminates row for measurement error rownames(results.table) <- NULL #********************************************************************** table.6 <- results.table # give results.table a name #********************************************************************** #Proportion of the variance explained by phylogeny and species id # phylogeny prop.phyl.var.fem <- mod$vcv[,1]/(mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10]) #females mean(prop.phyl.var.fem) HPDinterval(prop.phyl.var.fem) prop.phyl.var.male <- mod$vcv[,4]/ (mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13]) #males mean(prop.phyl.var.male) HPDinterval(prop.phyl.var.male) # species id prop.spp.var.fem <- mod$vcv[,5]/ (mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10]) #females mean(prop.spp.var.fem) HPDinterval(prop.spp.var.fem) prop.spp.var.male <- mod$vcv[,8]/ (mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13]) #males mean(prop.spp.var.male) HPDinterval(prop.spp.var.male) #residual prop.residual.fem <- mod$vcv[,10]/(mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10]) #females mean(prop.residual.fem) HPDinterval(prop.residual.fem) prop.residual.male <- mod$vcv[,13]/(mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13]) #males mean(prop.residual.male) HPDinterval(prop.residual.male) #correlations for phylogeny # getting correlations - (correlation = cov/(SD*SD)) spp_mat <- mod$vcv[,c(1, 3, 4)] attr(spp_mat, "dimnames")[[2]]<- c("f_V", "fm_cov", "m_V") # female variance, covariance, male variance spp_cor <- spp_mat[, "fm_cov"]/(sqrt(spp_mat[, "f_V"]*spp_mat[, "m_V"])) phylo.cor <- spp_cor mean(phylo.cor) HPDinterval(phylo.cor) summary(phylo.cor) #slope for phylogeny #computed as covariance divided by male variance slope_phylo <- spp_mat[,2]/spp_mat[,3] summary(slope_phylo) #correlations for species id # getting correlations - spp (correlation = cov/(SD*SD)) spp_mat <- mod$vcv[,c(5, 7, 8)] attr(spp_mat, "dimnames")[[2]]<- c("f_V", "fm_cov", "m_V") # female, male variance and cov spp_cor <- spp_mat[, "fm_cov"]/(sqrt(spp_mat[, "f_V"]*spp_mat[, "m_V"])) sp.id.cor <- spp_cor mean(sp.id.cor) HPDinterval(sp.id.cor) summary(sp.id.cor) #slope for species #computed as covariance divided by male variance slope_species <- spp_mat[,2]/spp_mat[,3] summary(slope_species) #correlations for residual # getting correlations - spp (correlation = cov/(SD*SD)) spp_mat <- mod$vcv[,c(10, 12, 13)] attr(spp_mat, "dimnames")[[2]]<- c("f_V", "fm_cov", "m_V") # female, male variance and cov spp_cor <- spp_mat[, "fm_cov"]/(sqrt(spp_mat[, "f_V"]*spp_mat[, "m_V"])) res.cor <- spp_cor mean(res.cor) HPDinterval(res.cor) summary(res.cor) #slope for units #computed as covariance divided by male variance slope_units <- spp_mat[,2]/spp_mat[,3] summary(slope_units) #Publication bias---- # predictions for fixed effects pred_fix <- as.vector(apply(mod$sol,2,mean)[1:mod$models[[1]]$Fixed$nfl] %*% t(mod$models[[1]]$X)) # predictions for random effects (including mev) pred_ran <- as.vector(apply(mod$sol,2,mean)[-(1:mod$models[[1]]$Fixed$nfl)] %*% t(mod$models[[1]]$Z)) # predictions for sampling variance raw_se <- c(datos$zr_se_fem, datos$zr_se_male) #pred_mev <- apply(mod$sol,2,mean)[(ncol(mod$sol)-nrow(datos)*2 + 1):ncol(mod$sol)]*raw_se # observed data raw_data <- c(datos$zr_f, datos$zr_m) # residuals pred_res <- raw_data - (pred_fix + pred_ran) # precision precision <- c(1/datos$zr_se_fem, 1/datos$zr_se_male) #funnel plots #females length(datos$overall_parameter) #getting the numbers of observations for the subsetted data plot(pred_res[1:109], precision[1:109]) #males plot(pred_res[110:218], precision[110:218]) #female residuals residuales <- pred_res[1:109] precis <- precision[1:109] #Trim and fill procedures metaR <- rma(yi= residuales, sei = 1/precis) TF <-trimfill(metaR,estimator="R0") TF funnel(TF, 'seinv') #Egger's test regtest(metaR, model="lm", predictor="sei") #male residuals residuales <- pred_res[110:218] precis <- precision[110:218] #Trim and fill procedures metaR <- rma(yi= residuales, sei = 1/precis) TF <-trimfill(metaR,estimator="R0") TF funnel(TF, 'seinv') regtest(metaR, model="lm", predictor="sei") #MODEL 7 - MODEL INCORPORATING SEXUAL DIMORPHISM FOR CONDITION DATA ONLY (Without removing outliers)---- # prior prior <- list(R = list(V = diag(2), nu = 2), G = list(G1 = list(V = diag(2) , nu = 2, alpha.mu=c(0,0), alpha.V=diag(2)*25^2), G2 = list(V = diag(2) , nu = 2, alpha.mu=c(0,0), alpha.V=diag(2)*25^2))) modelos <- NULL # setting up an empty list to store models class(modelos) <- 'list' #******************************** datos <- subset(dat, dimorphism_Cohend!="NA" & overall_parameter=="Body Condition") # SPECIFY DATA SET (no missing values) #******************************** #run loop. for (i in 1:50){ AinvTA.2<-inverseA(tree[[i]])$Ainv m1 <- MCMCglmm(cbind(zr_f, zr_m) ~ trait + dimorphism_Cohend + trait:dimorphism_Cohend + scale(year, scale= FALSE) -1, random = ~ us(trait):animal + us(trait):scientific_names, rcov = ~ us(trait):units, mev = c(datos$v_fem, datos$v_male), ginverse=list(animal = AinvTA.2), data = datos, family = rep('gaussian', 2), nitt=11000, thin=100, burnin=1000, pl = TRUE, pr = TRUE, prior = prior, verbose=FALSE) modelos[[i]] <- m1 print(i) } summary(modelos[[1]]) #eg output of the first model #processing multimodel results---- #combining posterior distributions across models all.m.vcv <- modelos[[1]]$VCV all.m.sol <- modelos[[1]]$Sol for (i in 2:50){ all.m.vcv <- rbind(all.m.vcv, modelos[[i]]$VCV) all.m.sol <- rbind(all.m.sol, modelos[[i]]$Sol) } #converting into mcmc objects all.m.vcv <- as.mcmc(all.m.vcv) all.m.sol <- as.mcmc(all.m.sol) #combine output into a list---- #this includes: #(a) a list with all 50 models #(b) the combined posterior distribution of all sol in (a) #(c) the combined posterior distribution of all vcv in (a) #*************************************************************************************** output.model7 <- (list(models=modelos, sol=all.m.sol, vcv=all.m.vcv)) # SPECIFY OUTPUT NAME #*************************************************************************************** #remove temporary files rm(modelos, all.m.sol, all.m.vcv) #summarising results---- #*************************************************************************************** mod <- output.model7 # SPECIFY NAME OF MODEL THAT WE WANT TO SUMMARISE #*************************************************************************************** #fixed effects #posterior means colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) #95%CIs HPDinterval(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) #random effects #posterior means colMeans(mod$vcv) #95%CIs HPDinterval(mod$vcv) #computing effective sample sizes effectiveSize(mod$vcv) effectiveSize(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) results.table <- rbind( #sol results data.frame(parameter=names(colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl])), post.mean=colMeans(mod$sol[,1:mod$models[[1]]$Fixed$nfl]), HPDinterval(mod$sol[,1:mod$models[[1]]$Fixed$nfl]), eff.sample.size=effectiveSize(mod$sol[,1:mod$models[[1]]$Fixed$nfl]) ), #vcv results data.frame(parameter=names(colMeans(mod$vcv)), post.mean=colMeans(mod$vcv), HPDinterval(mod$vcv), eff.sample.size=effectiveSize(mod$vcv)) ) results.table <- subset(results.table, parameter!="sqrt(mev):sqrt(mev).meta") # eliminates row for measurement error rownames(results.table) <- NULL #********************************************************************** table.7 <- results.table # give results.table a name #********************************************************************** #Proportion of the variance explained by phylogeny and species id # phylogeny prop.phyl.var.fem <- mod$vcv[,1]/(mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10]) #females mean(prop.phyl.var.fem) HPDinterval(prop.phyl.var.fem) prop.phyl.var.male <- mod$vcv[,4]/ (mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13]) #males mean(prop.phyl.var.male) HPDinterval(prop.phyl.var.male) # species id prop.spp.var.fem <- mod$vcv[,5]/ (mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10]) #females mean(prop.spp.var.fem) HPDinterval(prop.spp.var.fem) prop.spp.var.male <- mod$vcv[,8]/ (mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13]) #males mean(prop.spp.var.male) HPDinterval(prop.spp.var.male) #residual prop.residual.fem <- mod$vcv[,10]/(mod$vcv[,1] + mod$vcv[,5] + mod$vcv[,10]) #females mean(prop.residual.fem) HPDinterval(prop.residual.fem) prop.residual.male <- mod$vcv[,13]/(mod$vcv[,4] + mod$vcv[,8] + mod$vcv[,13]) #males mean(prop.residual.male) HPDinterval(prop.residual.male) #correlations for phylogeny # getting correlations - (correlation = cov/(SD*SD)) spp_mat <- mod$vcv[,c(1, 3, 4)] attr(spp_mat, "dimnames")[[2]]<- c("f_V", "fm_cov", "m_V") # female variance, covariance, male variance spp_cor <- spp_mat[, "fm_cov"]/(sqrt(spp_mat[, "f_V"]*spp_mat[, "m_V"])) phylo.cor <- spp_cor mean(phylo.cor) HPDinterval(phylo.cor) summary(phylo.cor) #slope for phylogeny #computed as covariance divided by male variance slope_phylo <- spp_mat[,2]/spp_mat[,3] summary(slope_phylo) #correlations for species id # getting correlations - spp (correlation = cov/(SD*SD)) spp_mat <- mod$vcv[,c(5, 7, 8)] attr(spp_mat, "dimnames")[[2]]<- c("f_V", "fm_cov", "m_V") # female, male variance and cov spp_cor <- spp_mat[, "fm_cov"]/(sqrt(spp_mat[, "f_V"]*spp_mat[, "m_V"])) sp.id.cor <- spp_cor mean(sp.id.cor) HPDinterval(sp.id.cor) summary(sp.id.cor) #slope for species #computed as covariance divided by male variance slope_species <- spp_mat[,2]/spp_mat[,3] summary(slope_species) #correlations for residual # getting correlations - spp (correlation = cov/(SD*SD)) spp_mat <- mod$vcv[,c(10, 12, 13)] attr(spp_mat, "dimnames")[[2]]<- c("f_V", "fm_cov", "m_V") # female, male variance and cov spp_cor <- spp_mat[, "fm_cov"]/(sqrt(spp_mat[, "f_V"]*spp_mat[, "m_V"])) res.cor <- spp_cor mean(res.cor) HPDinterval(res.cor) summary(res.cor) #slope for units #computed as covariance divided by male variance slope_units <- spp_mat[,2]/spp_mat[,3] summary(slope_units) #Publication bias---- # predictions for fixed effects pred_fix <- as.vector(apply(mod$sol,2,mean)[1:mod$models[[1]]$Fixed$nfl] %*% t(mod$models[[1]]$X)) # predictions for random effects (including mev) pred_ran <- as.vector(apply(mod$sol,2,mean)[-(1:mod$models[[1]]$Fixed$nfl)] %*% t(mod$models[[1]]$Z)) # predictions for sampling variance raw_se <- c(datos$zr_se_fem, datos$zr_se_male) #pred_mev <- apply(mod$sol,2,mean)[(ncol(mod$sol)-nrow(datos)*2 + 1):ncol(mod$sol)]*raw_se # observed data raw_data <- c(datos$zr_f, datos$zr_m) # residuals pred_res <- raw_data - (pred_fix + pred_ran) # precision precision <- c(1/datos$zr_se_fem, 1/datos$zr_se_male) #funnel plots #females length(datos$overall_parameter) #getting the numbers of observations for the subsetted data plot(pred_res[1:107], precision[1:107]) #males plot(pred_res[108:214], precision[108:214]) #female residuals residuales <- pred_res[1:107] precis <- precision[1:107] #Trim and fill procedures metaR <- rma(yi= residuales, sei = 1/precis) TF <-trimfill(metaR,estimator="R0") TF funnel(TF, 'seinv') #Egger's test regtest(metaR, model="lm", predictor="sei") #male residuals residuales <- pred_res[108:214] precis <- precision[108:214] #Trim and fill procedures metaR <- rma(yi= residuales, sei = 1/precis) TF <-trimfill(metaR,estimator="R0") TF funnel(TF, 'seinv') regtest(metaR, model="lm", predictor="sei")