hypothesis test - Part 4
Anastasia Shavrova
11/11/2021
Sample Sizes
Data reloading
full_es <- read.csv("full_model_es2.csv")
str(full_es)## 'data.frame': 940 obs. of 29 variables:
## $ parent : int 75 75 75 75 75 75 75 75 67 67 ...
## $ node : int 25 25 25 25 25 25 25 25 18 18 ...
## $ branch.length : num 0.0204 0.0204 0.0204 0.0204 0.0204 ...
## $ spcs : chr "Agelenopsis pennsylvanica" "Agelenopsis pennsylvanica" "Agelenopsis pennsylvanica" "Agelenopsis pennsylvanica" ...
## $ genus : chr "Agelenopsis" "Agelenopsis" "Agelenopsis" "Agelenopsis" ...
## $ ID : chr "ft016" "ft016" "ft016" "ft016" ...
## $ year : int 2012 2012 2012 2012 2012 2012 2012 2012 2019 2019 ...
## $ es_ID : chr "es_102" "es_103" "es_104" "es_100" ...
## $ female_ID : chr "f_20" "f_20" "f_20" "f_20" ...
## $ male_ID : chr "m_23" "m_24" "m_25" "m_22" ...
## $ treatment_ID : chr "trt_92" "trt_93" "trt_94" "trt_90" ...
## $ cannibal : chr "female" "female" "female" "female" ...
## $ mate_status : chr "virgin" "virgin" "virgin" "virgin" ...
## $ behaviour : chr "cannibal" "cannibal" "cannibal" "cannibal" ...
## $ categories : chr NA NA NA NA ...
## $ general_response2: chr "cannibal diet" "cannibal diet" "cannibal diet" "cannibal personality" ...
## $ exprt : chr "lab" "lab" "lab" "lab" ...
## $ habitat : chr "field" "field" "field" "field" ...
## $ data : chr "graph" "graph" "graph" "graph" ...
## $ occur : chr "pre" "pre" "pre" "pre" ...
## $ yi : num -0.6848 -0.341 0.0254 -1.0066 -0.0307 ...
## $ vi : num 0.1728 0.1602 0.1551 0.037 0.0332 ...
## $ adaptive : chr "Y" "Y" "Y" "N" ...
## $ aggressive : chr "Y" "Y" "Y" "Y" ...
## $ parental : chr "N" "N" "N" "N" ...
## $ choice : chr "N" "N" "N" "N" ...
## $ es_type : chr "OR_hedge_g_est_n01" "OR_hedge_g_est_n02" "OR_hedge_g_est_n03" "SMD_estimated_sample_size01" ...
## $ precision : num 2.41 2.5 2.54 5.2 5.49 ...
## $ world : chr "new" "new" "new" "new" ...Categorical moderators: * behaviour * experiment environment * habitat origin * which sex was the cannibal
summ.mod.all <- full_es %>%
summarise( # Calculate the number of effect sizes, studies and species for the main categorical variables
`Studies` = n_distinct(ID),
`Genus` = n_distinct(genus),
`Species` = n_distinct(spcs),
`New World Species` = n_distinct(spcs[world=="new"]),
`Old World Species` = n_distinct(spcs[world=="old"]),
`Effect sizes` = n_distinct(es_ID),
`Shared treatment measurements` = n_distinct(treatment_ID),
`Shared female measurements` = n_distinct(female_ID),
`Study (cannibal is virgin)` = n_distinct(ID[mate_status=="virgin"]),
`Study (cannibal is mated)` = n_distinct(ID[mate_status=="mated"]),
`Study (unknown mating status of cannibal)` = n_distinct(ID[mate_status=="unknown"])
)
kable(summ.mod.all) %>% kable_styling("striped", position="left")| Studies | Genus | Species | New World Species | Old World Species | Effect sizes | Shared treatment measurements | Shared female measurements | Study (cannibal is virgin) | Study (cannibal is mated) | Study (unknown mating status of cannibal) |
|---|---|---|---|---|---|---|---|---|---|---|
| 111 | 32 | 50 | 20 | 30 | 911 | 829 | 348 | 88 | 30 | 23 |
# per study
summ.mod.studies <- full_es %>%
summarise(
`Studies (behaviour attack)` = n_distinct(ID[behaviour=="attack"]),
`Studies (behaviour cannibalism)` = n_distinct(ID[behaviour=="cannibal"]),
`Studies (experiment in lab)` = n_distinct(ID[exprt=="lab"]),
`Studies (experiment in field)` = n_distinct(ID[exprt=="field"]),
`Studies (reared in lab)` = n_distinct(ID[exprt=="lab"]),
`Studies (reared in the field)` = n_distinct(ID[exprt=="field"]),
`Studies (female cannibals)` = n_distinct(ID[cannibal=="female"]),
`Studies (male cannibals)` = n_distinct(ID[cannibal=="male"]),
)
kable(summ.mod.studies) %>% kable_styling("striped", position="left") | Studies (behaviour attack) | Studies (behaviour cannibalism) | Studies (experiment in lab) | Studies (experiment in field) | Studies (reared in lab) | Studies (reared in the field) | Studies (female cannibals) | Studies (male cannibals) |
|---|---|---|---|---|---|---|---|
| 9 | 104 | 97 | 12 | 97 | 12 | 107 | 6 |
#per effect sizes
summ.mod.es <- full_es %>%
summarise(
`Effect Sizes (behaviour attack)` = n_distinct(es_ID[behaviour=="attack"]),
`Effect Sizes (behaviour cannibalism)` = n_distinct(es_ID[behaviour=="cannibal"]),
`Effect Sizes (experiment in lab)` = n_distinct(es_ID[exprt=="lab"]),
`Effect Sizes (experiment in field)` = n_distinct(es_ID[exprt=="field"]),
`Effect Sizes (experiment place unknown)` = n_distinct(es_ID[exprt=="NA"]),
`Effect Sizes (reared in lab)` = n_distinct(es_ID[exprt=="lab"]),
`Effect Sizes (reared in the field)` = n_distinct(es_ID[exprt=="field"]),
`Effect Sizes (unknown rearing of spider)` = n_distinct(es_ID[exprt=="NA"]),
`Effect Sizes (female cannibals)` = n_distinct(es_ID[cannibal=="female"]),
`Effect Sizes (male cannibals)` = n_distinct(es_ID[cannibal=="male"]),
)
kable(summ.mod.es) %>% kable_styling("striped", position="left") | Effect Sizes (behaviour attack) | Effect Sizes (behaviour cannibalism) | Effect Sizes (experiment in lab) | Effect Sizes (experiment in field) | Effect Sizes (experiment place unknown) | Effect Sizes (reared in lab) | Effect Sizes (reared in the field) | Effect Sizes (unknown rearing of spider) | Effect Sizes (female cannibals) | Effect Sizes (male cannibals) |
|---|---|---|---|---|---|---|---|---|---|
| 50 | 861 | 840 | 62 | 0 | 840 | 62 | 0 | 877 | 34 |
Data Exploration
Species <- full_es %>%
count(spcs)
genus <- full_es %>%
count(genus)
world <- full_es %>%
group_by(spcs) %>%
count(world)Data
full_es <- read.csv("full_model_es2.csv")
str(full_es)## 'data.frame': 940 obs. of 29 variables:
## $ parent : int 75 75 75 75 75 75 75 75 67 67 ...
## $ node : int 25 25 25 25 25 25 25 25 18 18 ...
## $ branch.length : num 0.0204 0.0204 0.0204 0.0204 0.0204 ...
## $ spcs : chr "Agelenopsis pennsylvanica" "Agelenopsis pennsylvanica" "Agelenopsis pennsylvanica" "Agelenopsis pennsylvanica" ...
## $ genus : chr "Agelenopsis" "Agelenopsis" "Agelenopsis" "Agelenopsis" ...
## $ ID : chr "ft016" "ft016" "ft016" "ft016" ...
## $ year : int 2012 2012 2012 2012 2012 2012 2012 2012 2019 2019 ...
## $ es_ID : chr "es_102" "es_103" "es_104" "es_100" ...
## $ female_ID : chr "f_20" "f_20" "f_20" "f_20" ...
## $ male_ID : chr "m_23" "m_24" "m_25" "m_22" ...
## $ treatment_ID : chr "trt_92" "trt_93" "trt_94" "trt_90" ...
## $ cannibal : chr "female" "female" "female" "female" ...
## $ mate_status : chr "virgin" "virgin" "virgin" "virgin" ...
## $ behaviour : chr "cannibal" "cannibal" "cannibal" "cannibal" ...
## $ categories : chr NA NA NA NA ...
## $ general_response2: chr "cannibal diet" "cannibal diet" "cannibal diet" "cannibal personality" ...
## $ exprt : chr "lab" "lab" "lab" "lab" ...
## $ habitat : chr "field" "field" "field" "field" ...
## $ data : chr "graph" "graph" "graph" "graph" ...
## $ occur : chr "pre" "pre" "pre" "pre" ...
## $ yi : num -0.6848 -0.341 0.0254 -1.0066 -0.0307 ...
## $ vi : num 0.1728 0.1602 0.1551 0.037 0.0332 ...
## $ adaptive : chr "Y" "Y" "Y" "N" ...
## $ aggressive : chr "Y" "Y" "Y" "Y" ...
## $ parental : chr "N" "N" "N" "N" ...
## $ choice : chr "N" "N" "N" "N" ...
## $ es_type : chr "OR_hedge_g_est_n01" "OR_hedge_g_est_n02" "OR_hedge_g_est_n03" "SMD_estimated_sample_size01" ...
## $ precision : num 2.41 2.5 2.54 5.2 5.49 ...
## $ world : chr "new" "new" "new" "new" ...full_es$spcs <- as.factor(full_es$spcs)
resolved_name <- tnrs_match_names(names=levels(full_es$spcs), context_name="Animals")
# Tree couldn't find: Trichonephila edulis, Trichonephila fenestrata, Trichonephila senegalensis because in tree of life still under "nephila"
#over 2 matches in rotl
matches.2 <- resolved_name %>%
filter(number_matches != 1)
#Hogna Helluo is ok
tree_es <- tol_induced_subtree(ott_ids = resolved_name$ott_id, label_format="name")## Warning in collapse_single_cpp(ances = tree$edge[, 1], desc = tree$edge[, :
## printing of extremely long output is truncated## Warning in collapse_single_cpp(ances = tree$edge[, 1], desc = tree$edge[, :
## printing of extremely long output is truncated## Warning in collapse_singles(tr, show_progress): Dropping singleton nodes
## with labels: mrcaott948ott572426, mrcaott948ott362702, mrcaott948ott3559852,
## mrcaott948ott693356, mrcaott948ott3559212, mrcaott948ott3560052,
## mrcaott948ott7473, mrcaott948ott177047, mrcaott948ott162088,
## mrcaott948ott180587, mrcaott948ott522289, mrcaott948ott219965,
## mrcaott948ott45861, Salticidae, mrcaott948ott29589, mrcaott948ott68490,
## mrcaott948ott254622, mrcaott948ott29585, mrcaott948ott124151,
## mrcaott948ott223069, mrcaott948ott272539, mrcaott948ott207505,
## mrcaott948ott73705, mrcaott948ott207502, mrcaott948ott203431, Pelleninae,
## mrcaott948ott105528, Habronattus, mrcaott8673ott471324, Evarcha,
## mrcaott4522ott93785, mrcaott4522ott3544199, mrcaott4522ott3544233,
## mrcaott4522ott6171006, mrcaott4522ott7686162, mrcaott4522ott3548025,
## Dendryphantinae, Phidippus, Thomisidae, Misumena, mrcaott138984ott3554226,
## mrcaott138984ott301864, mrcaott301864ott4694689, mrcaott301864ott3553908,
## mrcaott301864ott4696272, mrcaott301864ott3553911, mrcaott301864ott389252,
## mrcaott301864ott432787, mrcaott301864ott4696542, mrcaott301864ott4694281,
## mrcaott301864ott360021, Micaria, mrcaott7880ott608495, Allocosa,
## Aglaoctenus, mrcaott62344ott63321, Pisaura, Pisaurina, mrcaott33580ott33585,
## mrcaott33585ott504334, mrcaott33585ott220920, Hololena, Agelenopsis,
## mrcaott7883ott453835, mrcaott7883ott29040, mrcaott7883ott164006, Linyphiidae,
## Erigoninae, Diplocephalus, Theridiidae, Metellina, Araneus, Cyrtophora (genus
## in Opisthokonta), Larinioides, Phonognatha, Alpaida, Larinia, Austracantha,
## Nephilengys, Trichonephilatree_es$tip_label <- strip_ott_ids(tree_es$tip.label, remove_underscores=TRUE)
is.binary(tree_es)## [1] FALSE#resolving polytomies
tree_es_2 <- multi2di(tree_es)
###Make tree binary
is.binary(tree_es_2)## [1] TRUEtree_es_2$node.label <- NULL#Computing branch lengths & Final Tree
#computing branch lengths
branches2 <- compute.brlen(tree_es_2, method="Grafen", power=1)
is.ultrametric(branches2)## [1] TRUE#saving phylogenetic matrix
phylo_vcv2 <- vcv.phylo(branches2, corr=T)
# Generate variance-covariance matrices
#Calculating variance-covariance matrices for treatment_ID (assuming 0.5 corr)
make_VCV_matrix <- function(data, V, cluster, obs, type = c("vcv", "cor"), rho = 0.5) {
type <- match.arg(type)
if (missing(data)) {
stop("Must specify dataframe via 'data' argument.")
}
if (missing(V)) {
stop("Must specify name of the variance variable via 'V' argument.")
}
if (missing(cluster)) {
stop("Must specify name of the clustering variable via 'cluster' argument.")
}
if (missing(obs)) {
obs <- 1:length(V)
}
if (missing(type)) {
type <- "vcv"
}
new_matrix <- matrix(0, nrow = dim(data)[1], ncol = dim(data)[1]) #make empty matrix of the same size as data length
rownames(new_matrix) <- data[, obs]
colnames(new_matrix) <- data[, obs]
# find start and end coordinates for the subsets
shared_coord <- which(data[, cluster] %in% data[duplicated(data[, cluster]),
cluster] == TRUE)
# matrix of combinations of coordinates for each experiment with shared control
combinations <- do.call("rbind", tapply(shared_coord, data[shared_coord, cluster],
function(x) t(utils::combn(x, 2))))
if (type == "vcv") {
# calculate covariance values between values at the positions in shared_list and
# place them on the matrix
for (i in 1:dim(combinations)[1]) {
p1 <- combinations[i, 1]
p2 <- combinations[i, 2]
p1_p2_cov <- rho * sqrt(data[p1, V]) * sqrt(data[p2, V])
new_matrix[p1, p2] <- p1_p2_cov
new_matrix[p2, p1] <- p1_p2_cov
}
diag(new_matrix) <- data[, V] #add the diagonal
}
if (type == "cor") {
# calculate covariance values between values at the positions in shared_list and
# place them on the matrix
for (i in 1:dim(combinations)[1]) {
p1 <- combinations[i, 1]
p2 <- combinations[i, 2]
p1_p2_cov <- rho
new_matrix[p1, p2] <- p1_p2_cov
new_matrix[p2, p1] <- p1_p2_cov
}
diag(new_matrix) <- 1 #add the diagonal of 1
}
return(new_matrix)
}Vmat_all <- make_VCV_matrix(data = full_es, V = "vi", cluster = "treatment_ID", obs = "es_ID")Effect size is Hedges g
Full Intercept model with all effect sizes and random effects
full_mod_int <- rma.mv(yi=yi~1, #intercept model
V=Vmat_all, #covariance matrix with treatment_ID
method="REML", #restricted maximum likelihood
random=list(~1|female_ID, #when females were reused
~1|male_ID, #when males were reused
~1|es_type, #type of effect size converted to hedges g
~1|es_ID), # all possible random effects
R = list(spcs = phylo_vcv2),
data=full_es,
sparse=TRUE)## Warning: Argument 'R' specified, but list name(s) not in 'random'.summary(full_mod_int)##
## Multivariate Meta-Analysis Model (k = 940; method: REML)
##
## logLik Deviance AIC BIC AICc ​
## -1869.3593 3738.7186 3748.7186 3772.9427 3748.7830
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 1.0290 1.0144 348 no female_ID
## sigma^2.2 0.0000 0.0003 381 no male_ID
## sigma^2.3 1.4316 1.1965 940 no es_type
## sigma^2.4 0.2128 0.4613 911 no es_ID
##
## Test for Heterogeneity:
## Q(df = 939) = 8346.6313, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub ​
## 0.0591 0.0792 0.7458 0.4558 -0.0962 0.2144
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1results_int <- mod_results(full_mod_int, group="es_ID", data=full_es)Total Heterogeneity and heterogeneity by random effects
i2_ml(full_mod_int)## I2_Total I2_female_ID I2_male_ID I2_es_type I2_es_ID
## 9.485555e+01 3.651102e+01 3.794058e-06 5.079533e+01 7.549208e+00Plot Intercept model
orchard_plot(full_mod_int, xlab="yi", group="es_ID", data=full_es) +
annotate(geom = "text",x = 1, y = 30, label = paste0("italic(I)^{2} == ", round(i2_ml(full_mod_int)[1] * 100, 2)), color = "black", parse = TRUE, size = 5) + geom_errorbarh(aes(xmin = results_int$lowerPR, xmax = results_int$upperPR), height = 0, show.legend = FALSE, size = 0.8, alpha = 0.6) +
geom_errorbarh(aes(xmin = results_int$lowerCL, xmax = results_int$upperCL), height = 0.04, show.legend = FALSE, size = 1) +
scale_colour_manual(values = "darkred")+ # change colour of data points
scale_fill_manual(values="darkorange1") + # change colour of the mean estimate +
theme(text = element_text(size = 20, colour = "black", hjust = 0.5),
axis.text.y = element_text(size = 15,colour = "black", hjust = 0.5),
legend.title = element_text(size=16)) +
scale_y_discrete(labels = "Overall mean")+ # rename y label
scale_size_continuous(range = c(1, 9))+ # change scaling of data points
ylab("Effect Size")## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
Plotting full model with all moderators
Not sure if this is the model to use, or to use the separate models done below
full_mod <- rma.mv(yi=yi, #intercept model
V=Vmat_all, #covariance matrix with treatment_ID
method="REML", #restricted maximum likelihood
random=list(~1|female_ID, #when females were reused
~1|male_ID, #when males were reused
~1|es_type, #type of effect size converted to hedges g
~1|es_ID), # all possible random effects
R = list(spcs = phylo_vcv2),
mod=~cannibal + mate_status + behaviour + general_response2 + exprt + habitat + occur + world,
data=full_es,
sparse=TRUE)## Warning: Argument 'R' specified, but list name(s) not in 'random'.## Warning: Rows with NAs omitted from model fitting.summary(full_mod) # too big a model? can't plot... ##
## Multivariate Meta-Analysis Model (k = 912; method: REML)
##
## logLik Deviance AIC BIC AICc ​
## -1402.2967 2804.5934 2872.5934 3035.1880 2875.4034
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.5403 0.7351 336 no female_ID
## sigma^2.2 0.0000 0.0002 366 no male_ID
## sigma^2.3 0.3956 0.6289 912 no es_type
## sigma^2.4 0.1306 0.3614 883 no es_ID
##
## Test for Residual Heterogeneity:
## QE(df = 882) = 5248.9235, p-val < .0001
##
## Test of Moderators (coefficients 2:30):
## QM(df = 29) = 219.5230, p-val < .0001
##
## Model Results:
##
## estimate se zval​
## intrcpt -1.7060 0.4717 -3.6167
## cannibalmale -0.9219 0.2777 -3.3198
## mate_statusunknown 0.4627 0.2093 2.2111
## mate_statusvirgin 0.3313 0.1507 2.1993
## behaviourcannibal -0.5073 0.2268 -2.2372
## general_response2cannibal age 1.9460 0.3757 5.1796
## general_response2cannibal condition, non-diet 1.7569 0.3528 4.9800
## general_response2cannibal diet 1.3799 0.3527 3.9129
## general_response2cannibal personality 1.8872 0.4249 4.4416
## general_response2cannibal size 2.0135 0.3472 5.7999
## general_response2copulation duration 1.9224 0.3498 5.4954
## general_response2courtship behaviour 2.2743 0.4461 5.0984
## general_response2courtship duration 1.7878 0.4090 4.3716
## general_response2fecundity 2.0472 0.3462 5.9134
## general_response2mate availability 1.3167 0.3676 3.5816
## general_response2offspring fitness 1.9977 0.3482 5.7365
## general_response2paternity 1.6568 0.3999 4.1427
## general_response2relatedness 1.7844 0.5661 3.1522
## general_response2SSD 1.9504 0.3818 5.1089
## general_response2victim age 1.6122 0.3687 4.3725
## general_response2victim condition 1.5105 0.3616 4.1771
## general_response2victim mating status 1.1320 0.3769 3.0030
## general_response2victim size 1.5387 0.3445 4.4666
## exprtlab 0.5587 0.1881 2.9704
## exprtunknown -0.0091 0.5109 -0.0177
## habitatlab -0.1003 0.1234 -0.8127
## habitatunknown 10.6877 1.2222 8.7445
## occurpre -0.2845 0.1109 -2.5661
## occurunknown 1.8363 1.3261 1.3847
## worldold -0.2253 0.1114 -2.0230
## pval ci.lb ci.ub
## intrcpt 0.0003 -2.6305 -0.7815 ***
## cannibalmale 0.0009 -1.4662 -0.3776 ***
## mate_statusunknown 0.0270 0.0526 0.8729 *
## mate_statusvirgin 0.0279 0.0361 0.6266 *
## behaviourcannibal 0.0253 -0.9518 -0.0629 *
## general_response2cannibal age <.0001 1.2096 2.6824 ***
## general_response2cannibal condition, non-diet <.0001 1.0654 2.4484 ***
## general_response2cannibal diet <.0001 0.6887 2.0712 ***
## general_response2cannibal personality <.0001 1.0544 2.7200 ***
## general_response2cannibal size <.0001 1.3331 2.6940 ***
## general_response2copulation duration <.0001 1.2367 2.6080 ***
## general_response2courtship behaviour <.0001 1.4000 3.1485 ***
## general_response2courtship duration <.0001 0.9863 2.5894 ***
## general_response2fecundity <.0001 1.3687 2.7258 ***
## general_response2mate availability 0.0003 0.5962 2.0372 ***
## general_response2offspring fitness <.0001 1.3151 2.6802 ***
## general_response2paternity <.0001 0.8729 2.4406 ***
## general_response2relatedness 0.0016 0.6749 2.8939 **
## general_response2SSD <.0001 1.2022 2.6987 ***
## general_response2victim age <.0001 0.8896 2.3349 ***
## general_response2victim condition <.0001 0.8018 2.2193 ***
## general_response2victim mating status 0.0027 0.3932 1.8708 **
## general_response2victim size <.0001 0.8635 2.2138 ***
## exprtlab 0.0030 0.1901 0.9274 **
## exprtunknown 0.9859 -1.0104 0.9923
## habitatlab 0.4164 -0.3421 0.1416
## habitatunknown <.0001 8.2922 13.0832 ***
## occurpre 0.0103 -0.5018 -0.0672 *
## occurunknown 0.1661 -0.7628 4.4354
## worldold 0.0431 -0.4436 -0.0070 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Total Heterogeneity and heterogeneity by random effects
i2_ml(full_mod)## I2_Total I2_female_ID I2_male_ID I2_es_type I2_es_ID
## 8.811908e+01 4.464567e+01 3.985771e-06 3.268416e+01 1.078924e+01Plot model - can’t model too big?
_____________________________By each moderator____________________
Moderator - Habitat
mod_habitat <- rma.mv(yi=yi, #intercept model
V=Vmat_all, #covariance matrix with treatment_ID
method="REML", #restricted maximum likelihood
random=list(~1|female_ID, #when females were reused
~1|male_ID, #when males were reused
~1|es_type, #type of effect size converted to hedges g
~1|es_ID), # all possible random effects
R = list(spcs = phylo_vcv2),
mod= ~habitat,
data=full_es,
sparse=TRUE)## Warning: Argument 'R' specified, but list name(s) not in 'random'.summary(mod_habitat)##
## Multivariate Meta-Analysis Model (k = 940; method: REML)
##
## logLik Deviance AIC BIC AICc ​
## -1845.6857 3691.3714 3705.3714 3739.2702 3705.4919
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.8606 0.9277 348 no female_ID
## sigma^2.2 0.0000 0.0004 381 no male_ID
## sigma^2.3 1.3754 1.1728 940 no es_type
## sigma^2.4 0.2139 0.4625 911 no es_ID
##
## Test for Residual Heterogeneity:
## QE(df = 937) = 7747.8697, p-val < .0001
##
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 42.7504, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub ​
## intrcpt 0.0650 0.0880 0.7381 0.4604 -0.1075 0.2374
## habitatlab -0.0855 0.1519 -0.5628 0.5736 -0.3833 0.2123
## habitatunknown 10.5600 1.6257 6.4957 <.0001 7.3737 13.7464 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1results_habitat <- mod_results(mod_habitat, mod="habitat", group="es_ID", data=full_es)
orchard_plot(results_habitat, mod="habitat", xlab="yi", group="es_ID", data=full_es)
Moderator - Experimental Habitat
mod_habitat_exp <- rma.mv(yi=yi, #intercept model
V=Vmat_all, #covariance matrix with treatment_ID
method="REML", #restricted maximum likelihood
random=list(~1|female_ID, #when females were reused
~1|male_ID, #when males were reused
~1|es_type, #type of effect size converted to hedges g
~1|es_ID), # all possible random effects
R = list(spcs = phylo_vcv2),
mod= ~exprt,
data=full_es,
sparse=TRUE)## Warning: Argument 'R' specified, but list name(s) not in 'random'.summary(mod_habitat_exp)##
## Multivariate Meta-Analysis Model (k = 940; method: REML)
##
## logLik Deviance AIC BIC AICc ​
## -1862.7564 3725.5129 3739.5129 3773.4116 3739.6334
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 1.0392 1.0194 348 no female_ID
## sigma^2.2 0.0000 0.0003 381 no male_ID
## sigma^2.3 1.4126 1.1885 940 no es_type
## sigma^2.4 0.2105 0.4588 911 no es_ID
##
## Test for Residual Heterogeneity:
## QE(df = 937) = 8277.1268, p-val < .0001
##
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 6.7203, p-val = 0.0347
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub ​
## intrcpt -0.1760 0.2473 -0.7117 0.4767 -0.6607 0.3087
## exprtlab 0.2233 0.2488 0.8977 0.3694 -0.2643 0.7109
## exprtunknown 1.7514 0.6756 2.5924 0.0095 0.4272 3.0756 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1results_habitat_exprt <- mod_results(mod_habitat_exp, mod="exprt", group="es_ID", data=full_es)
orchard_plot(results_habitat_exprt, mod="exprt", xlab="yi", group="es_ID", data=full_es)
Moderator - cannibal
mod_cannibal <- rma.mv(yi=yi, #intercept model
V=Vmat_all, #covariance matrix with treatment_ID
method="REML", #restricted maximum likelihood
random=list(~1|female_ID, #when females were reused
~1|male_ID, #when males were reused
~1|es_type, #type of effect size converted to hedges g
~1|es_ID), # all possible random effects
R = list(spcs = phylo_vcv2),
mod= ~cannibal,
data=full_es,
sparse=TRUE)## Warning: Argument 'R' specified, but list name(s) not in 'random'.summary(mod_cannibal)##
## Multivariate Meta-Analysis Model (k = 940; method: REML)
##
## logLik Deviance AIC BIC AICc ​
## -1864.9373 3729.8746 3741.8746 3770.9371 3741.9648
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.9900 0.9950 348 no female_ID
## sigma^2.2 0.0000 0.0003 381 no male_ID
## sigma^2.3 1.4470 1.2029 940 no es_type
## sigma^2.4 0.1975 0.4444 911 no es_ID
##
## Test for Residual Heterogeneity:
## QE(df = 938) = 8252.2082, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 5.7105, p-val = 0.0169
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub ​
## intrcpt 0.1049 0.0806 1.3014 0.1931 -0.0531 0.2629
## cannibalmale -0.8187 0.3426 -2.3897 0.0169 -1.4903 -0.1472 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1results_cannibal <- mod_results(mod_cannibal, mod="cannibal", group="es_ID", data=full_es)
orchard_plot(results_cannibal, mod="cannibal", xlab="yi", group="es_ID", data=full_es)
Moderator - mate_status
mod_status <- rma.mv(yi=yi,
V=Vmat_all, #covariance matrix with treatment_ID
method="REML", #restricted maximum likelihood
random=list(~1|female_ID, #when females were reused
~1|male_ID, #when males were reused
~1|es_type, #type of effect size converted to hedges g
~1|es_ID), # all possible random effects
R = list(spcs = phylo_vcv2),
mod= ~mate_status,
data=full_es,
sparse=TRUE)## Warning: Argument 'R' specified, but list name(s) not in 'random'.summary(mod_status)##
## Multivariate Meta-Analysis Model (k = 940; method: REML)
##
## logLik Deviance AIC BIC AICc ​
## -1865.7884 3731.5768 3745.5768 3779.4756 3745.6974
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 1.0369 1.0183 348 no female_ID
## sigma^2.2 0.0000 0.0004 381 no male_ID
## sigma^2.3 1.4277 1.1949 940 no es_type
## sigma^2.4 0.2191 0.4681 911 no es_ID
##
## Test for Residual Heterogeneity:
## QE(df = 937) = 8318.2542, p-val < .0001
##
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 0.6005, p-val = 0.7406
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub ​
## intrcpt 0.1217 0.1839 0.6619 0.5081 -0.2387 0.4820
## mate_statusunknown 0.0386 0.2451 0.1576 0.8748 -0.4418 0.5190
## mate_statusvirgin -0.1008 0.1983 -0.5082 0.6113 -0.4895 0.2879
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1results_status <- mod_results(mod_status, mod="mate_status", group="es_ID", data=full_es)
orchard_plot(results_status, mod="mate_status", xlab="yi", group="es_ID", data=full_es)
Moderator - occur
mod_occur <- rma.mv(yi=yi,
V=Vmat_all, #covariance matrix with treatment_ID
method="REML", #restricted maximum likelihood
random=list(~1|female_ID, #when females were reused
~1|male_ID, #when males were reused
~1|es_type, #type of effect size converted to hedges g
~1|es_ID), # all possible random effects
R = list(spcs = phylo_vcv2),
mod= ~occur,
data=full_es,
sparse=TRUE)## Warning: Argument 'R' specified, but list name(s) not in 'random'.summary(mod_occur)##
## Multivariate Meta-Analysis Model (k = 940; method: REML)
##
## logLik Deviance AIC BIC AICc ​
## -1864.9041 3729.8083 3743.8083 3777.7071 3743.9288
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 1.0520 1.0257 348 no female_ID
## sigma^2.2 0.0000 0.0003 381 no male_ID
## sigma^2.3 1.4331 1.1971 940 no es_type
## sigma^2.4 0.1988 0.4458 911 no es_ID
##
## Test for Residual Heterogeneity:
## QE(df = 937) = 8343.8524, p-val < .0001
##
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 2.6279, p-val = 0.2688
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub ​
## intrcpt 0.1330 0.0957 1.3899 0.1646 -0.0546 0.3206
## occurpre -0.2109 0.1458 -1.4464 0.1481 -0.4966 0.0749
## occurunknown 1.2610 1.8253 0.6908 0.4897 -2.3166 4.8385
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1results_occur <- mod_results(mod_occur, mod="occur", group="es_ID", data=full_es)
orchard_plot(results_occur, mod="occur", xlab="yi", group="es_ID", data=full_es)
Moderator - world
mod_world <- rma.mv(yi=yi,
V=Vmat_all, #covariance matrix with treatment_ID
method="REML", #restricted maximum likelihood
random=list(~1|female_ID, #when females were reused
~1|male_ID, #when males were reused
~1|es_type, #type of effect size converted to hedges g
~1|es_ID), # all possible random effects
R = list(spcs = phylo_vcv2),
mod= ~world,
data=full_es,
sparse=TRUE)## Warning: Argument 'R' specified, but list name(s) not in 'random'.summary(mod_world)##
## Multivariate Meta-Analysis Model (k = 940; method: REML)
##
## logLik Deviance AIC BIC AICc ​
## -1867.6785 3735.3569 3747.3569 3776.4194 3747.4472
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 1.0358 1.0178 348 no female_ID
## sigma^2.2 0.0000 0.0003 381 no male_ID
## sigma^2.3 1.4363 1.1985 940 no es_type
## sigma^2.4 0.2080 0.4561 911 no es_ID
##
## Test for Residual Heterogeneity:
## QE(df = 938) = 8339.2814, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 0.0311, p-val = 0.8601
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub ​
## intrcpt 0.0421 0.1250 0.3364 0.7366 -0.2030 0.2871
## worldold 0.0259 0.1471 0.1763 0.8601 -0.2623 0.3142
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1results_world <- mod_results(mod_world, mod="world", group="es_ID", data=full_es)
orchard_plot(results_world, mod="world", xlab="yi", group="es_ID", data=full_es)
Moderator - behaviour
mod_behaviour <- rma.mv(yi=yi,
V=Vmat_all, #covariance matrix with treatment_ID
method="REML", #restricted maximum likelihood
random=list(~1|female_ID, #when females were reused
~1|male_ID, #when males were reused
~1|es_type, #type of effect size converted to hedges g
~1|es_ID), # all possible random effects
R = list(spcs = phylo_vcv2),
mod= ~behaviour,
data=full_es,
sparse=TRUE)## Warning: Argument 'R' specified, but list name(s) not in 'random'.summary(mod_behaviour)##
## Multivariate Meta-Analysis Model (k = 940; method: REML)
##
## logLik Deviance AIC BIC AICc ​
## -1866.6872 3733.3745 3745.3745 3774.4370 3745.4647
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 1.0386 1.0191 348 no female_ID
## sigma^2.2 0.0000 0.0003 381 no male_ID
## sigma^2.3 1.4164 1.1901 940 no es_type
## sigma^2.4 0.2213 0.4705 911 no es_ID
##
## Test for Residual Heterogeneity:
## QE(df = 938) = 8314.2772, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 2.1296, p-val = 0.1445
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub ​
## intrcpt 0.4487 0.2789 1.6089 0.1076 -0.0979 0.9952
## behaviourcannibal -0.4231 0.2899 -1.4593 0.1445 -0.9914 0.1452
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1results_behaviour <- mod_results(mod_behaviour, mod="behaviour", group="es_ID", data=full_es)
orchard_plot(results_behaviour, mod="behaviour", xlab="yi", group="es_ID", data=full_es)
Moderator - general_response2
- check which are spread over two areas
full_es_na_gen <- full_es %>%
drop_na(general_response2)
Vmat_na_gen <- make_VCV_matrix(data = full_es_na_gen, V = "vi", cluster = "treatment_ID", obs = "es_ID")
mod_response <- rma.mv(yi=yi,
V=Vmat_na_gen, #covariance matrix with treatment_ID and no NA's
method="REML", #restricted maximum likelihood
random=list(~1|female_ID, #when females were reused
~1|male_ID, #when males were reused
~1|es_type, #type of effect size converted to hedges g
~1|es_ID), # all possible random effects
R = list(spcs = phylo_vcv2),
mod= ~general_response2,
data=full_es_na_gen,
sparse=TRUE)## Warning: Argument 'R' specified, but list name(s) not in 'random'.summary(mod_response)##
## Multivariate Meta-Analysis Model (k = 912; method: REML)
##
## logLik Deviance AIC BIC AICc ​
## -1472.2841 2944.5683 2990.5683 3100.8438 2991.8387
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.7637 0.8739 336 no female_ID
## sigma^2.2 0.0066 0.0813 366 no male_ID
## sigma^2.3 0.4479 0.6693 912 no es_type
## sigma^2.4 0.1330 0.3646 883 no es_ID
##
## Test for Residual Heterogeneity:
## QE(df = 893) = 5955.4592, p-val < .0001
##
## Test of Moderators (coefficients 2:19):
## QM(df = 18) = 86.4875, p-val < .0001
##
## Model Results:
##
## estimate se zval​
## intrcpt -1.6796 0.3450 -4.8686
## general_response2cannibal age 1.9998 0.3998 5.0026
## general_response2cannibal condition, non-diet 1.7262 0.3727 4.6320
## general_response2cannibal diet 1.3998 0.3752 3.7305
## general_response2cannibal personality 1.9315 0.4490 4.3013
## general_response2cannibal size 2.0706 0.3677 5.6316
## general_response2copulation duration 2.0172 0.3715 5.4303
## general_response2courtship behaviour 2.3645 0.4784 4.9423
## general_response2courtship duration 1.8679 0.4297 4.3469
## general_response2fecundity 2.1046 0.3684 5.7127
## general_response2mate availability 1.2766 0.3932 3.2464
## general_response2offspring fitness 2.0112 0.3698 5.4380
## general_response2paternity 1.6412 0.4243 3.8682
## general_response2relatedness 1.7643 0.6075 2.9044
## general_response2SSD 1.9352 0.4017 4.8180
## general_response2victim age 1.6092 0.3941 4.0834
## general_response2victim condition 1.6055 0.3829 4.1932
## general_response2victim mating status 1.0714 0.3962 2.7040
## general_response2victim size 1.5643 0.3647 4.2888
## pval ci.lb ci.ub
## intrcpt <.0001 -2.3558 -1.0034 ***
## general_response2cannibal age <.0001 1.2163 2.7833 ***
## general_response2cannibal condition, non-diet <.0001 0.9958 2.4567 ***
## general_response2cannibal diet 0.0002 0.6643 2.1352 ***
## general_response2cannibal personality <.0001 1.0514 2.8116 ***
## general_response2cannibal size <.0001 1.3500 2.7913 ***
## general_response2copulation duration <.0001 1.2892 2.7453 ***
## general_response2courtship behaviour <.0001 1.4268 3.3022 ***
## general_response2courtship duration <.0001 1.0257 2.7101 ***
## general_response2fecundity <.0001 1.3825 2.8266 ***
## general_response2mate availability 0.0012 0.5059 2.0474 **
## general_response2offspring fitness <.0001 1.2863 2.7361 ***
## general_response2paternity 0.0001 0.8096 2.4728 ***
## general_response2relatedness 0.0037 0.5737 2.9549 **
## general_response2SSD <.0001 1.1480 2.7224 ***
## general_response2victim age <.0001 0.8368 2.3815 ***
## general_response2victim condition <.0001 0.8551 2.3559 ***
## general_response2victim mating status 0.0069 0.2948 1.8479 **
## general_response2victim size <.0001 0.8494 2.2791 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1results_response <- mod_results(mod_response, mod="general_response2", group="es_ID", data=full_es_na_gen)
orchard_plot(results_response, mod="general_response2", xlab="yi", group="es_ID", data=full_es)