Start

rm(list=ls())
setwd("~/Documents/PhD Thesis/Heritability meta-analysis")


library(readxl)
library(tidyverse); theme_set(theme_light())
library(metafor)
library(tidytext)
library(shades)
library(ggridges)
library(ggplotify)
require(ggnewscale)
library(grid)
library(gridExtra)
library(gt)
library(cowplot)
library(RColorBrewer)
library(viridis)
source('Functions/useful_functions.R')
source('Functions/mlm.variance.distribution.R')

# Load complete data:
load("Data/heritability_estimates_processed.RData")
data <- data.full %>% 
    filter(trait != "gamete contribution") %>% 
    mutate(trait = fct_drop(trait))
dat_s <- data %>% 
    filter(!(trait %in% c("photochemistry", "immune response"))) %>% 
    mutate(trait = factor(trait))
dat_sh <- data %>% 
    filter(!(trait %in% c("photochemistry", "nutrient content", "immune response"))) %>% 
    mutate(trait = factor(trait))
dat_shg <- data %>% 
    filter(trait %in% c("symbiont community", "survival", "nutrient content", "growth", "bleaching")) %>% 
    filter(!(growth.form %in% c("columnar", "encrusting"))) %>%
    mutate(trait = factor(trait), growth.form = factor(growth.form))
dat_tm <- data %>% 
    filter(!is.na(temp.diff)) %>%
    filter(!(trait %in% c("symbiont community", "morphology", "gene expression"))) %>% 
    mutate(trait = factor(trait),
           temp.manip = temp.manip.plus.others,
           temp.s = sqrt(temp.diff))
dat_tm_nz <- dat_tm %>% filter(temp.diff >0) %>%
    mutate(trait = factor(trait), 
           stage = factor(stage), 
           h2 = factor(h2),
           growth.form = factor(growth.form))

# Load final models:
load("Models/final.model.Rdata")
load("Models/final.model.A.Rdata")
load("Models/final.model.B.Rdata")
load("Models/final.model.C.Rdata")
load("Models/final.model.Di.Rdata")
load("Models/final.model.Dii.Rdata")
load("Models/final.model.Diii.Rdata")

# Objects for plotting:
mod.coef <- update(final.model, val.log ~ trait-1, data=data.full) %>% coef() # includes gamete contribution
trait.colors <- gg_color_hue(length(levels(data.full$trait)))
trait.levels <- gsub("trait","",names(mod.coef))[order(mod.coef)]

Supplementary tables and figures

Fig. S1

funnel(final.model, yaxis="seinv")

Figure S1. Funnel plot for overall analysis

Table S1-2, Fig. S2 - Overall analysis: main effects only (excluding temperature due to incompleteness across all data)

x <- read.csv(file="Suppl Tables/TableS1.csv")
y <- read.csv(file="Suppl Tables/TableS2.csv") %>% head(5)
data.subset = "the complete dataset"
beyond.opt.mod="trait + life stage + heritability type + growth form"

n=1

make_modsel_table(x,y,data.subset, beyond.opt.mod,n=n); n=n+1

Table S1. Model selection for random and fixed effect structures using the complete dataset, ordered by lowest AICc.

		df	AICc	ΔAICc
Top random effects^*
1	Estimate ID nested in study ID	15	112.2	0.00
2	Estimate ID only	15	118.2	5.97
3	Estimate ID nested in species	15	121.6	9.36
4	Study ID only	15	190.3	78.12
5	Species only	15	311.9	199.72
6	No random effect	15	353.6	241.37
Top fixed effects^†
1	trait	8	100.8	0.00
2	trait + heritability type	9	101.1	0.37
3	trait + life stage	10	101.9	1.09
4	trait + heritability type + life stage	11	103.4	2.61
5	trait + growth form	12	109.0	8.25
^* Models fit via REML using a 'beyond optimal' fixed effects structure of: trait + life stage + heritability type + growth form. ^† Models fit via ML using the top random effects structure, seen above.

make_coef_table(final.model, y, data.subset=data.subset, n=n); n=n+1

Figure S2. Model coefficients for the top fixed effect structures of the overall analysis

Table S2. Model coefficients (estimated on the log[h2+0.2]-scale) for the top 3 models within 2 ΔAICc, using the complete dataset.^*

	trait model		trait + heritability type model		trait + life stage model
	est ± SE (95%CI)	test statistic	est ± SE (95%CI)	test statistic	est ± SE (95%CI)	test statistic
intrcpt	−1.15±0.22 (−1.58, −0.71)	t = −5.27 p = 0.00	−1.07±0.22 (−1.50, −0.63)	t = −4.90 p = 0.00	−1.15±0.26 (−1.66, −0.65)	t = −4.52 p = 0.00
photochemistry	0.36±0.18 (−0.00, 0.73)	t = 1.98 p = 0.05	0.35±0.18 (−0.02, 0.72)	t = 1.90 p = 0.06	0.36±0.18 (−0.00, 0.73)	t = 1.97 p = 0.05
growth	0.39±0.21 (−0.04, 0.81)	t = 1.81 p = 0.07	0.38±0.21 (−0.05, 0.80)	t = 1.76 p = 0.08	0.41±0.21 (−0.02, 0.84)	t = 1.91 p = 0.06
nutrientcontent	0.53±0.24 (0.06, 1.01)	t = 2.25 p = 0.03	0.52±0.24 (0.05, 0.99)	t = 2.18 p = 0.03	0.54±0.24 (0.07, 1.01)	t = 2.27 p = 0.03
bleaching	0.57±0.24 (0.08, 1.05)	t = 2.33 p = 0.02	0.60±0.24 (0.12, 1.08)	t = 2.46 p = 0.02	0.58±0.24 (0.10, 1.06)	t = 2.41 p = 0.02
symbiontcommunity	0.67±0.31 (0.06, 1.29)	t = 2.18 p = 0.03	0.78±0.31 (0.16, 1.39)	t = 2.53 p = 0.01	0.70±0.31 (0.09, 1.31)	t = 2.29 p = 0.02
morphology	0.79±0.31 (0.18, 1.40)	t = 2.57 p = 0.01	0.79±0.30 (0.19, 1.39)	t = 2.63 p = 0.01	0.84±0.30 (0.24, 1.44)	t = 2.77 p = 0.01
immuneresponse	0.94±0.25 (0.44, 1.45)	t = 3.73 p = 0.00	0.92±0.25 (0.42, 1.43)	t = 3.63 p = 0.00	0.94±0.25 (0.43, 1.44)	t = 3.70 p = 0.00
survival	1.13±0.22 (0.70, 1.57)	t = 5.16 p = 0.00	1.13±0.22 (0.69, 1.57)	t = 5.15 p = 0.00	1.15±0.22 (0.71, 1.59)	t = 5.21 p = 0.00
hnarrow	± (, )	t = p =	−0.24±0.17 (−0.58, 0.09)	t = −1.43 p = 0.16	± (, )	t = p =
juvenile	± (, )	t = p =	± (, )	t = p =	−0.28±0.23 (−0.75, 0.18)	t = −1.21 p = 0.23
adult	± (, )	t = p =	± (, )	t = p =	0.10±0.20 (−0.30, 0.50)	t = 0.50 p = 0.62
^* Models fit via REML using the optimal random effects structure selected previously.

Table S3-4, Fig. S3 - Sub-analysis A: examine trait x life stage interaction

x <- read.csv(file="Suppl Tables/TableS3.csv")
y <- read.csv(file="Suppl Tables/TableS4.csv") %>% head(5)
data.subset = "a subset to examine interactions between trait type and life stage"
beyond.opt.mod="trait x life stage + heritability type + growth form"

make_modsel_table(x,y,data.subset, beyond.opt.mod,n=n); n=n+1

Table S3. Model selection for random and fixed effect structures using a subset to examine interactions between trait type and life stage, ordered by lowest AICc.

		df	AICc	ΔAICc
Top random effects^*
1	Estimate ID only	22	127.7	0.00
2	Study ID only	22	130.9	3.19
3	Estimate ID nested in study ID	22	134.0	6.26
4	Estimate ID nested in species	22	135.6	7.85
5	No random effect	22	141.8	14.04
6	Species only	22	148.1	20.35
Top fixed effects^†
1	trait x life stage + heritability type	18	80.1	0.00
2	trait x life stage	17	82.2	2.03
3	trait x life stage + growth form	21	89.0	8.84
4	trait x life stage + heritability type + growth form	22	90.1	9.98
5	trait + life stage + growth form	12	91.9	11.70
^* Models fit via REML using a 'beyond optimal' fixed effects structure of: trait x life stage + heritability type + growth form. ^† Models fit via ML using the top random effects structure, seen above.

make_coef_table(final.model.A, y, data.subset=data.subset, n=n); n=n+1

Figure S3. Model coefficients for the top fixed effect structures of sub-analysis A

Table S4. Model coefficients (estimated on the log[h2+0.2]-scale) for the top model (no other models within 2 ΔAICc), using a subset to examine interactions between trait type and life stage.^*

	trait x life stage + heritability type model
	est ± SE (95%CI)	test statistic
intrcpt	−0.27±0.53 (−1.34, 0.79)	t = −0.51 p = 0.61
growth	0.44±0.17 (0.10, 0.78)	t = 2.60 p = 0.01
nutrientcontent	0.06±0.56 (−1.07, 1.18)	t = 0.10 p = 0.92
bleaching	−0.32±0.58 (−1.48, 0.85)	t = −0.55 p = 0.59
symbiontcommunity	0.18±0.56 (−0.94, 1.30)	t = 0.32 p = 0.75
morphology	0.17±0.58 (−0.98, 1.33)	t = 0.30 p = 0.77
survival	0.26±0.54 (−0.82, 1.33)	t = 0.48 p = 0.63
juvenile	0.16±0.17 (−0.18, 0.51)	t = 0.94 p = 0.35
adult	−0.66±0.55 (−1.77, 0.45)	t = −1.19 p = 0.24
hnarrow	−0.27±0.12 (−0.51, −0.03)	t = −2.27 p = 0.03
growthjuvenile	−1.44±0.59 (−2.62, −0.25)	t = −2.44 p = 0.02
bleachingjuvenile	−0.67±0.35 (−1.38, 0.04)	t = −1.90 p = 0.06
symbiontcommunityjuvenile	−0.51±0.69 (−1.89, 0.87)	t = −0.74 p = 0.46
morphologyjuvenile	−0.47±0.35 (−1.17, 0.22)	t = −1.36 p = 0.18
nutrientcontentadult	−0.16±0.59 (−1.35, 1.03)	t = −0.27 p = 0.79
bleachingadult	1.16±0.62 (−0.09, 2.40)	t = 1.86 p = 0.07
symbiontcommunityadult	0.68±0.61 (−0.55, 1.91)	t = 1.11 p = 0.27
morphologyadult	0.36±0.65 (−0.94, 1.67)	t = 0.56 p = 0.58
survivaladult	0.33±0.57 (−0.81, 1.47)	t = 0.57 p = 0.57
^* Models fit via REML using the optimal random effects structure selected previously.

Table S5-6, Fig. S4 - Sub-analysis B: examine trait x heritability interaction

x <- read.csv(file="Suppl Tables/TableS5.csv")
y <- read.csv(file="Suppl Tables/TableS6.csv") %>% head(5)
data.subset = "a subset to examine interactions between trait type and heritability type"
beyond.opt.mod="trait x heritability type + trait x life stage + growth form"

make_modsel_table(x,y,data.subset, beyond.opt.mod,n=n); n=n+1

Table S5. Model selection for random and fixed effect structures using a subset to examine interactions between trait type and heritability type, ordered by lowest AICc.

		df	AICc	ΔAICc
Top random effects^*
1	Study ID only	23	161.1	0.00
2	Estimate ID only	23	161.2	0.04
3	No random effect	23	168.1	6.98
4	Estimate ID nested in study ID	23	172.8	11.66
5	Estimate ID nested in species	23	174.5	13.32
6	Species only	23	178.7	17.57
Top fixed effects^†
1	trait x life stage	15	78.2	0.00
2	trait x life stage + heritability type	16	80.8	2.68
3	trait x life stage + trait x heritability type	19	84.4	6.28
4	trait + heritability type	11	87.7	9.54
5	trait x heritability type + life stage	13	89.2	11.03
^* Models fit via REML using a 'beyond optimal' fixed effects structure of: trait x heritability type + trait x life stage + growth form. ^† Models fit via ML using the top random effects structure, seen above.

make_coef_table(final.model.B, y, data.subset=data.subset, n=n); n=n+1

Figure S4. Model coefficients for the top fixed effect structures of sub-analysis B

Table S6. Model coefficients (estimated on the log[h2+0.2]-scale) for the top model (no other models within 2 ΔAICc), using a subset to examine interactions between trait type and heritability type.^*

	trait x life stage model
	est ± SE (95%CI)	test statistic
intrcpt	−0.50±0.55 (−1.60, 0.60)	t = −0.91 p = 0.37
growth	0.58±0.16 (0.26, 0.90)	t = 3.61 p = 0.00
bleaching	0.03±0.60 (−1.18, 1.24)	t = 0.04 p = 0.96
symbiontcommunity	0.13±0.60 (−1.07, 1.33)	t = 0.21 p = 0.83
morphology	0.49±0.58 (−0.68, 1.66)	t = 0.84 p = 0.41
survival	0.34±0.56 (−0.78, 1.46)	t = 0.61 p = 0.54
juvenile	0.25±0.21 (−0.18, 0.67)	t = 1.18 p = 0.24
adult	−0.57±0.58 (−1.73, 0.60)	t = −0.98 p = 0.33
growthjuvenile	−1.50±0.59 (−2.68, −0.32)	t = −2.55 p = 0.01
bleachingjuvenile	−1.07±0.30 (−1.68, −0.47)	t = −3.55 p = 0.00
symbiontcommunityjuvenile	−0.69±0.69 (−2.08, 0.71)	t = −0.99 p = 0.33
morphologyjuvenile	−0.82±0.37 (−1.56, −0.08)	t = −2.21 p = 0.03
bleachingadult	0.89±0.64 (−0.40, 2.18)	t = 1.38 p = 0.17
symbiontcommunityadult	0.60±0.64 (−0.69, 1.89)	t = 0.93 p = 0.36
morphologyadult	0.12±0.69 (−1.27, 1.50)	t = 0.17 p = 0.87
survivaladult	0.69±0.59 (−0.50, 1.88)	t = 1.17 p = 0.25
^* Models fit via REML using the optimal random effects structure selected previously.

Table S7-8, Fig. S5 - Sub-analysis C: examine trait x growth form interaction

x <- read.csv(file="Suppl Tables/TableS7.csv")
y <- read.csv(file="Suppl Tables/TableS8.csv") %>% head(5)
data.subset = "a subset to examine interactions between trait type and growth form"
beyond.opt.mod="trait x growth form + trait x life stage + heritability type"

make_modsel_table(x,y,data.subset, beyond.opt.mod,n=n); n=n+1

Table S7. Model selection for random and fixed effect structures using a subset to examine interactions between trait type and growth form, ordered by lowest AICc.

		df	AICc	ΔAICc
Top random effects^*
1	Study ID only	18	120.6	0.00
2	Estimate ID only	18	121.2	0.51
3	No random effect	18	122.5	1.87
4	Estimate ID nested in study ID	18	132.8	12.12
5	Species only	18	133.1	12.42
6	Estimate ID nested in species	18	134.2	13.58
Top fixed effects^†
1	trait x life stage	12	58.0	0.00
2	trait x life stage + heritability type	13	60.5	2.57
3	trait x life stage + growth form	14	64.9	6.95
4	trait x life stage + trait x growth form + heritability type	18	65.0	7.01
5	trait x life stage + trait x growth form	17	65.3	7.32
^* Models fit via REML using a 'beyond optimal' fixed effects structure of: trait x growth form + trait x life stage + heritability type. ^† Models fit via ML using the top random effects structure, seen above.

make_coef_table(final.model.C, y, data.subset=data.subset, n=n); n=n+1

Figure S5. Model coefficients for the top fixed effect structures of sub-analysis C

Table S8. Model coefficients (estimated on the log[h2+0.2]-scale) for the top model (no other models within 2 ΔAICc), using a subset to examine interactions between trait type and growth form.^*

	trait x life stage model
	est ± SE (95%CI)	test statistic
intrcpt	−0.62±0.19 (−1.00, −0.23)	t = −3.23 p = 0.00
nutrientcontent	0.49±0.25 (−0.01, 0.99)	t = 1.99 p = 0.05
bleaching	0.15±0.29 (−0.44, 0.73)	t = 0.50 p = 0.62
symbiontcommunity	0.25±0.31 (−0.37, 0.86)	t = 0.80 p = 0.43
survival	0.46±0.13 (0.20, 0.72)	t = 3.52 p = 0.00
juvenile	−0.56±0.27 (−1.10, −0.01)	t = −2.07 p = 0.04
adult	0.14±0.20 (−0.26, 0.54)	t = 0.70 p = 0.49
bleachingjuvenile	−0.27±0.34 (−0.96, 0.42)	t = −0.79 p = 0.43
symbiontcommunityjuvenile	0.12±0.71 (−1.32, 1.56)	t = 0.17 p = 0.87
survivaljuvenile	0.81±0.16 (0.48, 1.14)	t = 4.96 p = 0.00
nutrientcontentadult	−0.74±0.24 (−1.22, −0.25)	t = −3.08 p = 0.00
bleachingadult	0.18±0.31 (−0.45, 0.81)	t = 0.57 p = 0.57
symbiontcommunityadult	−0.11±0.33 (−0.78, 0.56)	t = −0.33 p = 0.74
^* Models fit via REML using the optimal random effects structure selected previously.

Table S9-10, Fig. S6 - Sub-analysis D-i: examine trait x temperature difference interaction

x <- read.csv(file="Suppl Tables/TableS9.csv")
y <- read.csv(file="Suppl Tables/TableS10.csv") %>% head(5)
data.subset = "a subset to examine interactions between trait type and temperature difference (covariate; square-root transformed)"
beyond.opt.mod="trait x temp diff + life stage + heritability type + growth form"

make_modsel_table(x,y,data.subset, beyond.opt.mod,n=n); n=n+1

Table S9. Model selection for random and fixed effect structures using a subset to examine interactions between trait type and temperature difference (covariate; square-root transformed), ordered by lowest AICc.

		df	AICc	ΔAICc
Top random effects^*
1	Estimate ID nested in study ID	18	106.3	0.00
2	Estimate ID only	18	109.9	3.66
3	Estimate ID nested in species	18	115.8	9.52
4	Study ID only	18	122.2	15.94
5	No random effect	18	212.4	106.16
6	Species only	18	212.4	106.20
Top fixed effects^†
1	trait + heritability type	6	69.4	0.00
2	trait	5	69.9	0.50
3	trait x temp diff	11	71.3	1.89
4	trait + heritability type + temp diff	7	71.7	2.32
5	trait + temp diff	6	72.4	3.00
^* Models fit via REML using a 'beyond optimal' fixed effects structure of: trait x temp diff + life stage + heritability type + growth form. ^† Models fit via ML using the top random effects structure, seen above.

make_coef_table(final.model.Di, y, data.subset=data.subset, n=n); n=n+1

Figure S6. Model coefficients for the top fixed effect structures of sub-analysis Dii

Table S10. Model coefficients (estimated on the log[h2+0.2]-scale) for the top 3 models within 2 ΔAICc, using a subset to examine interactions between trait type and temperature difference (covariate; square-root transformed).^*

	trait + heritability type model		trait model		trait x temp diff model
	est ± SE (95%CI)	test statistic	est ± SE (95%CI)	test statistic	est ± SE (95%CI)	test statistic
intrcpt	−0.75±0.16 (−1.07, −0.43)	t = −4.70 p = 0.00	−0.85±0.16 (−1.18, −0.53)	t = −5.25 p = 0.00	−0.83±0.22 (−1.27, −0.38)	t = −3.73 p = 0.00
growth	0.05±0.15 (−0.25, 0.35)	t = 0.33 p = 0.75	0.04±0.15 (−0.26, 0.34)	t = 0.26 p = 0.80	0.04±0.23 (−0.41, 0.50)	t = 0.20 p = 0.85
nutrientcontent	0.19±0.18 (−0.16, 0.54)	t = 1.08 p = 0.29	0.19±0.18 (−0.16, 0.54)	t = 1.07 p = 0.29	0.15±0.24 (−0.33, 0.63)	t = 0.62 p = 0.53
bleaching	0.17±0.20 (−0.24, 0.58)	t = 0.84 p = 0.41	0.12±0.20 (−0.29, 0.52)	t = 0.58 p = 0.56	0.44±0.30 (−0.16, 1.04)	t = 1.46 p = 0.15
immuneresponse	0.58±0.20 (0.19, 0.98)	t = 2.97 p = 0.00	0.59±0.20 (0.19, 0.98)	t = 2.99 p = 0.00	0.15±0.26 (−0.38, 0.68)	t = 0.55 p = 0.58
survival	0.82±0.16 (0.51, 1.13)	t = 5.25 p = 0.00	0.80±0.16 (0.48, 1.11)	t = 5.10 p = 0.00	0.86±0.22 (0.42, 1.30)	t = 3.93 p = 0.00
temps	± (, )	t = p =	± (, )	t = p =	−0.01±0.13 (−0.28, 0.26)	t = −0.06 p = 0.95
growthtemps	± (, )	t = p =	± (, )	t = p =	−0.02±0.16 (−0.33, 0.30)	t = −0.10 p = 0.92
nutrientcontenttemps	± (, )	t = p =	± (, )	t = p =	0.05±0.17 (−0.29, 0.40)	t = 0.31 p = 0.76
bleachingtemps	± (, )	t = p =	± (, )	t = p =	−0.25±0.20 (−0.65, 0.14)	t = −1.29 p = 0.20
immuneresponsetemps	± (, )	t = p =	± (, )	t = p =	0.53±0.20 (0.13, 0.94)	t = 2.64 p = 0.01
survivaltemps	± (, )	t = p =	± (, )	t = p =	−0.03±0.15 (−0.33, 0.27)	t = −0.20 p = 0.85
hnarrow	−0.39±0.23 (−0.85, 0.06)	t = −1.72 p = 0.09	± (, )	t = p =	± (, )	t = p =
^* Models fit via REML using the optimal random effects structure selected previously.

Table S11-12, Fig. S7 - Sub-analysis D-ii: examine trait x temperature difference (binary) interaction

x <- read.csv(file="Suppl Tables/TableS11.csv")
y <- read.csv(file="Suppl Tables/TableS12.csv") %>% head(5)
data.subset = "a subset to examine interactions between trait type and temperature manipulation (binary)"
beyond.opt.mod="trait x temp manip + life stage + heritability type + growth form"

make_modsel_table(x,y,data.subset, beyond.opt.mod,n=n); n=n+1

Table S11. Model selection for random and fixed effect structures using a subset to examine interactions between trait type and temperature manipulation (binary), ordered by lowest AICc.

		df	AICc	ΔAICc
Top random effects^*
1	Estimate ID nested in study ID	18	106.3	0.00
2	Estimate ID only	18	108.0	1.66
3	Estimate ID nested in species	18	113.9	7.52
4	Study ID only	18	120.0	13.67
5	Species only	18	178.1	71.74
6	No random effect	18	181.5	75.16
Top fixed effects^†
1	trait + heritability type	6	69.4	0.00
2	trait	5	69.9	0.50
3	trait + heritability type + temp manip	7	70.2	0.81
4	trait + temp manip	6	71.2	1.85
5	trait x temp manip	11	71.7	2.30
^* Models fit via REML using a 'beyond optimal' fixed effects structure of: trait x temp manip + life stage + heritability type + growth form. ^† Models fit via ML using the top random effects structure, seen above.

make_coef_table(final.model.Dii, y, data.subset=data.subset, n=n); n=n+1

Figure S7. Model coefficients for the top fixed effect structures of sub-analysis Di

Table S12. Model coefficients (estimated on the log[h2+0.2]-scale) for the top 4 models within 2 ΔAICc, using a subset to examine interactions between trait type and temperature manipulation (binary).^*

	trait + heritability type model		trait model		trait + heritability type + temp manip model		trait + temp manip model
	est ± SE (95%CI)	test statistic	est ± SE (95%CI)	test statistic	est ± SE (95%CI)	test statistic	est ± SE (95%CI)	test statistic
intrcpt	−0.75±0.16 (−1.07, −0.43)	t = −4.70 p = 0.00	−0.85±0.16 (−1.18, −0.53)	t = −5.25 p = 0.00	−0.79±0.17 (−1.14, −0.44)	t = −4.56 p = 0.00	−0.88±0.18 (−1.24, −0.52)	t = −4.89 p = 0.00
growth	0.05±0.15 (−0.25, 0.35)	t = 0.33 p = 0.75	0.04±0.15 (−0.26, 0.34)	t = 0.26 p = 0.80	0.06±0.15 (−0.24, 0.36)	t = 0.40 p = 0.69	0.04±0.15 (−0.26, 0.35)	t = 0.29 p = 0.77
nutrientcontent	0.19±0.18 (−0.16, 0.54)	t = 1.08 p = 0.29	0.19±0.18 (−0.16, 0.54)	t = 1.07 p = 0.29	0.19±0.18 (−0.17, 0.54)	t = 1.05 p = 0.30	0.19±0.18 (−0.17, 0.54)	t = 1.05 p = 0.30
bleaching	0.17±0.20 (−0.24, 0.58)	t = 0.84 p = 0.41	0.12±0.20 (−0.29, 0.52)	t = 0.58 p = 0.56	0.17±0.20 (−0.23, 0.58)	t = 0.85 p = 0.40	0.12±0.20 (−0.29, 0.52)	t = 0.59 p = 0.56
immuneresponse	0.58±0.20 (0.19, 0.98)	t = 2.97 p = 0.00	0.59±0.20 (0.19, 0.98)	t = 2.99 p = 0.00	0.58±0.20 (0.18, 0.97)	t = 2.91 p = 0.01	0.58±0.20 (0.19, 0.98)	t = 2.94 p = 0.00
survival	0.82±0.16 (0.51, 1.13)	t = 5.25 p = 0.00	0.80±0.16 (0.48, 1.11)	t = 5.10 p = 0.00	0.81±0.16 (0.50, 1.13)	t = 5.18 p = 0.00	0.79±0.16 (0.48, 1.11)	t = 5.03 p = 0.00
hnarrow	−0.39±0.23 (−0.85, 0.06)	t = −1.72 p = 0.09	± (, )	t = p =	−0.41±0.23 (−0.86, 0.05)	t = −1.79 p = 0.08	± (, )	t = p =
temps	± (, )	t = p =	± (, )	t = p =	0.03±0.05 (−0.08, 0.14)	t = 0.57 p = 0.57	0.02±0.05 (−0.09, 0.12)	t = 0.36 p = 0.72
^* Models fit via REML using the optimal random effects structure selected previously.

Table S13-14, Fig. S8 - Sub-analysis D-iii: examine trait x temperature difference interaction, excluding zero values

x <- read.csv(file="Suppl Tables/TableS13.csv")
y <- read.csv(file="Suppl Tables/TableS14.csv") %>% head(5)
data.subset = "a subset to examine interactions between trait type and temperature difference (covariate; square-root transformed, **non-zero values only**)"
beyond.opt.mod="trait x temp diff + life stage + heritability type + growth form"

make_modsel_table(x,y,data.subset, beyond.opt.mod,n=n); n=n+1

Table S13. Model selection for random and fixed effect structures using a subset to examine interactions between trait type and temperature difference (covariate; square-root transformed, non-zero values only), ordered by lowest AICc.

		df	AICc	ΔAICc
Top random effects^*
1	Estimate ID only	15	120.5	0.00
2	Study ID only	15	131.1	10.63
3	Estimate ID nested in study ID	15	135.0	14.51
4	Estimate ID nested in species	15	137.3	16.80
5	No random effect	15	146.0	25.56
6	Species only	15	150.2	29.72
Top fixed effects^†
1	trait + life stage	7	56.8	0.00
2	life stage	2	57.1	0.30
3	trait + temp diff + life stage	8	57.4	0.56
4	life stage + heritability type	3	59.4	2.53
5	life stage + temp diff	3	59.7	2.84
^* Models fit via REML using a 'beyond optimal' fixed effects structure of: trait x temp diff + life stage + heritability type + growth form. ^† Models fit via ML using the top random effects structure, seen above.

make_coef_table(final.model.Diii, y, data.subset=data.subset, n=n); n=n+1

Figure S8. Model coefficients for the top fixed effect structures of sub-analysis Diii

Table S14. Model coefficients (estimated on the log[h2+0.2]-scale) for the top 3 models within 2 ΔAICc, using a subset to examine interactions between trait type and temperature difference (covariate; square-root transformed, non-zero values only).^*

	trait + life stage model		life stage model		trait + temp diff + life stage model
	est ± SE (95%CI)	test statistic	est ± SE (95%CI)	test statistic	est ± SE (95%CI)	test statistic
intrcpt	−0.35±0.21 (−0.78, 0.09)	t = −1.62 p = 0.11	−0.26±0.15 (−0.56, 0.05)	t = −1.71 p = 0.09	0.07±0.35 (−0.64, 0.79)	t = 0.21 p = 0.83
juvenile	−0.70±0.21 (−1.12, −0.27)	t = −3.33 p = 0.00	−0.80±0.21 (−1.22, −0.38)	t = −3.82 p = 0.00	−0.77±0.21 (−1.20, −0.34)	t = −3.62 p = 0.00
adult	−0.26±0.19 (−0.64, 0.13)	t = −1.34 p = 0.19	−0.31±0.17 (−0.65, 0.04)	t = −1.80 p = 0.08	−0.35±0.20 (−0.75, 0.05)	t = −1.77 p = 0.08
growth	−0.07±0.21 (−0.50, 0.35)	t = −0.36 p = 0.72	± (, )	t = p =	−0.07±0.20 (−0.49, 0.34)	t = −0.36 p = 0.72
nutrientcontent	−0.18±0.21 (−0.61, 0.24)	t = −0.88 p = 0.38	± (, )	t = p =	−0.23±0.21 (−0.65, 0.19)	t = −1.12 p = 0.27
bleaching	−0.21±0.24 (−0.70, 0.28)	t = −0.89 p = 0.38	± (, )	t = p =	−0.26±0.24 (−0.74, 0.23)	t = −1.07 p = 0.29
immuneresponse	0.46±0.27 (−0.08, 1.01)	t = 1.72 p = 0.09	± (, )	t = p =	0.47±0.26 (−0.06, 1.01)	t = 1.80 p = 0.08
survival	0.38±0.19 (−0.01, 0.77)	t = 1.96 p = 0.06	± (, )	t = p =	0.45±0.20 (0.06, 0.85)	t = 2.31 p = 0.03
temps	± (, )	t = p =	± (, )	t = p =	−0.19±0.13 (−0.46, 0.07)	t = −1.49 p = 0.14
^* Models fit via REML using the optimal random effects structure selected previously.

Table S15

Table S15. Number of heritability estimates calculated across relatedness type (rows) and heritability model type (columns). Note that the heritabilities calculated may be either narrow-sense or broad-sense depending on the experimental design and whether or not clones were used.

require(stringr); require(knitr)
data %>%
    mutate(relatedness.info = str_remove(relatedness.info, " \\(.*\\)"),
           relatedness.info = str_replace(relatedness.info, "inferred pedigree", "pedigree"),
           model.type = str_remove(model.type, " \\(.*\\)"),
           model.type = str_replace(model.type, "animal model approximation", "animal model")) %>%
    select(relatedness.info, model.type) %>%
    table %>%
    knitr::kable()

	animal model	ANOVA	Ritland multiple regression model
clones	33	33	0
genetic markers	2	0	6
pedigree	20	0	0

Fig. S9

Plot of residuals by trait, while accounting for heritability type, life stage, and whether or not the individuals were reared in a common shared environment from larvae. Note the clear confound of all adults not being raised in controlled conditions from larvae, and almost all juveniles and larvae being raised in controlled environments.

data$resid <- resid(final.model)
pS9 <- data %>%
    ggplot(aes(trait, resid, col=reared.from.larvae.in.common.env)) +
    geom_point() +
    geom_hline(yintercept=0, linetype = 2, color="grey") +
    facet_grid(stage~h2, switch="y") +
    labs(x = "Trait type", y = "Residual value") +
    theme(#axis.text.y.left=element_blank(),
          axis.text.x=element_text(angle = 60, hjust=1),
          #axis.ticks.y.left=element_blank(),
          strip.text.y = element_text(),
          panel.spacing = unit(0.1, "lines"),
          panel.grid.major=element_blank(),
          panel.grid.minor=element_blank(),
          legend.text.align = 0.5,
          strip.background = element_rect(fill=NA, color = "grey"),
          strip.text = element_text(color="black"),
          strip.placement = "outside") +
        scale_color_discrete(name = "Reared in a common\nshared environment\nfrom larvae?")

pS9

Fig. S9. Model residual values from the overall analysis using trait type alone as a fixed effect vs. life stage (row facets), heritability type (column facets), and whether or not individuals were reared from larvae in a common shared environment (colour).

setEPS()
postscript("Figures/FigS9.eps", family="Helvetica", width=6, height=6)
pS9 # 1000 x 700
dev.off()

quartz_off_screen 
                2

Fig. S10

plot_grid(pS10)

Fig. S10. Heritability estimate by temperature model type

setEPS()
postscript("Figures/FigS10.eps", family="Helvetica", width=9, height=13)
plot_grid(pS10) # 1000 x 700
dev.off()

quartz_off_screen 
                2

Fig. S11

Similar to Fig. S9, but for the heritability model used to estimate heritability for each study.

data <- data %>%
    mutate(model.type2 = str_remove(model.type, " \\(.*\\)"),
           model.type2 = str_replace(model.type2, "animal model approximation", "animal model"))

pS11 <- data %>%
    ggplot(aes(trait, resid, col=model.type2)) +
    geom_point() +
    geom_hline(yintercept=0, linetype = 2, color="grey") +
    facet_grid(stage~h2, switch="y") +
    labs(x = "Trait type", y = "Residual value") +
    scale_color_discrete(name = "Model type",
                         labels = c("animal model", "ANOVA-method", "Ritland regression\n(all from a single study)")) +
    theme(#axis.text.y.left=element_blank(),
          axis.text.x=element_text(angle = 60, hjust=1),
          #axis.ticks.y.left=element_blank(),
          strip.text.y = element_text(),
          panel.spacing = unit(0.1, "lines"),
          panel.grid.major=element_blank(),
          panel.grid.minor=element_blank(),
          legend.text.align = 0.5,
          strip.background = element_rect(fill=NA, color = "grey"),
          strip.text = element_text(color="black"),
          strip.placement = "outside")

pS11

Fig. S11. Model residual values from the overall analysis using trait type alone as a fixed effect vs. life stage (row facets), heritability type (column facets), and the type of model used to estimate the heritability coefficient (colour).

setEPS()
postscript("Figures/FigS11.eps", family="Helvetica", width=6, height=6)
pS11 # 1000 x 700
dev.off()

quartz_off_screen 
                2

Reported heritabilities

Code to provide summary numbers for literature review results paragraph of text.

# Review results paragraph:
summarise.coverage("h2", "study", data.full)
summarise.coverage("h2", "est.id", data.full)$Y

summarise.coverage("trait", "study", data.full)
summarise.coverage("trait", "est.id", data.full)

summarise.coverage("study", "est.id", data.full)$Y

summarise.coverage("stage", "est.id", data.full)$Y
summarise.coverage("stage", "study", data.full)$Y
head(summarise.coverage("stage", "study", data.full)$X, 3)
data.full %>% filter(study == "Quigley et al. 2017") %>%
    select(study, species, trait, stage)

summarise.coverage("h2", "est.id", data.full)$Y
summarise.coverage("h2", "study", data.full)$Y
head(summarise.coverage("h2", "study", data.full)$X, 3)

summarise.coverage("growth.form", "est.id", data.full)$Y
summarise.coverage("growth.form", "study", data.full)$Y
head(summarise.coverage("growth.form", "study", data.full)$X, 3)

data.full %>% pull(study) %>% unique %>% length
data.full %>% filter(!is.na(temp.diff)) %>% pull(study) %>% unique %>% length

data.full %>% pull(est.id) %>% unique %>% length
data.full %>% filter(!is.na(temp.diff)) %>% pull(est.id) %>% unique %>% length

hist(dat_tm$temp.diff, breaks=10)

data.full %>% filter(temp.diff == 0) %>% pull(est.id) %>% unique %>% length
data.full %>% filter(temp.diff == 0) %>% pull(study) %>% unique %>% length

data.full %>% filter(temp.diff != 0) %>% pull(est.id) %>% unique %>% length
data.full %>% filter(temp.diff != 0) %>% pull(study) %>% unique %>% length

summarise.coverage("trait", "est.id", dat_tm)$Y
summarise.coverage("trait", "study", dat_tm)$Y
head(summarise.coverage("growth.form", "study", dat_tm)$X, 3)


data.full %>% group_by(species) %>% count(sort=T) # n=19 species
data.full %>% group_by(reproduction, symbiont.transmission) %>% count(sort=T)

data.full %>% 
    group_by(species, reproduction, symbiont.transmission) %>% 
    count(sort=T) %>% select(-n) %>%
    ungroup() %>%
    group_by(reproduction, symbiont.transmission) %>%
    count(sort=T)

Note: code above not run to avoid excessive output tables.

Main model results

Overall model: ~ trait type

Final model summary

# Uses nearly full data, minus one trait with only a single estimate
final.model <- rma.mv(yi=val.log ~ trait, V=sv.log, 
                      random = ~1|study/est.id, data=data, method="REML", tdist=T)
mlm.variance.distribution(final.model)

	% of total variance	I2
Level 1	8.831568	—
Level 2	34.418637	34.42
Level 3	56.749795	56.75

Model I^2 and R^2:

# Main model:

# Get estimate of I^2, residual heterogeneity:
W <- diag(1/data$sv.log) # create initial weighting matrix
X <- model.matrix(final.model)
P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W
I2 <- 100 * final.model$sigma2 / (sum(final.model$sigma2) + (final.model$k-final.model$p)/sum(diag(P)))
cat(paste0("Overall residual heterogeneity, I^2 = ",round(sum(I2),1), "%"))

Overall residual heterogeneity, I^2 = 84.2%

cat(paste0("Between-study residual heterogeneity, I^2 = ",round(I2[1],1), "%"))

Between-study residual heterogeneity, I^2 = 52.4%

cat(paste0("Within-study residual heterogeneity, I^2 = ",round(I2[2],1), "%"))

Within-study residual heterogeneity, I^2 = 31.8%

# Get estimate of R^2: 
Model1 <- rma.mv(yi=val.log ~ 1, V=sv.log, 
                      random = ~1|study/est.id, data=data, method="ML", tdist=T)
Model2 <- rma.mv(yi=val.log ~ trait, V=sv.log, 
                      random = ~1|study/est.id, data=data, method="ML", tdist=T)
(Model1$sigma2[1] - Model2$sigma2[1]) / Model1$sigma2[1] # R^2 = 0 or -0.14

[1] -0.3159911

(Model1$sigma2[2] - Model2$sigma2[2]) / Model1$sigma2[2] # R^2 = 0.56

[1] 0.5728472

(sum(Model1$sigma2) - sum(Model2$sigma2)) / sum(Model1$sigma2) # R^2 = 0.29 overall?

[1] 0.2698959

Model estimates

Model estimates ~ trait type:

overall_trait_mod_est <- update(final.model, .~.-1, ) %>% # remove intercept
    summary %>% coef %>% # extract model coefficients + CI
    rownames_to_column("trait") %>%
    mutate(trait = str_remove(trait, "trait")) %>%
    rename(est = estimate)
overall_trait_mod_est_backtrans <- overall_trait_mod_est %>%
    mutate(across(.cols = c(est, ci.lb, ci.ub), # apply an anonymous fn
                  .fns = function(x){exp(x)-0.2} )) %>% # to back-transform
    select(-se, -tval, -pval) %>%
    mutate(ci.ub = ifelse(ci.ub >1, 1, ci.ub))
write.csv(overall_trait_mod_est, file="Model Estimates/overall_trait_mod_est.csv", row.names=F)
# write.csv(overall_trait_mod_est_backtrans, file="Model Estimates/overall_trait_mod_est_backtrans.csv", row.names=F)

overall_trait_mod_est

trait	est	se	tval	pval	ci.lb	ci.ub
gene expression	-1.1454361	0.2174874	-5.2666774	0.0000010	-1.5778594	-0.7130128
photochemistry	-0.7823753	0.1524744	-5.1311923	0.0000018	-1.0855352	-0.4792154
growth	-0.7570162	0.1309643	-5.7803256	0.0000001	-1.0174083	-0.4966241
nutrient content	-0.6119458	0.1609494	-3.8021018	0.0002692	-0.9319562	-0.2919354
bleaching	-0.5764261	0.1492929	-3.8610425	0.0002195	-0.8732603	-0.2795919
symbiont community	-0.4712294	0.2267899	-2.0778239	0.0407426	-0.9221484	-0.0203104
morphology	-0.3535242	0.2334017	-1.5146599	0.1335684	-0.8175894	0.1105409
immune response	-0.2019735	0.1906822	-1.0592157	0.2925023	-0.5811007	0.1771537
survival	-0.0109606	0.1218698	-0.0899367	0.9285491	-0.2532703	0.2313492

Model estimates ~ trait type + h2, ~ trait type + life stage:

overall_trait_mod_est_addh2 <- update(final.model, .~trait + h2 -1, ) %>% # remove intercept
    summary %>% coef %>% # extract model coefficients + CI
    rownames_to_column("trait") %>%
    mutate(trait = str_remove(trait, "trait")) %>%
    rename(est = estimate)
overall_trait_mod_est_addh2_backtrans <- overall_trait_mod_est_addh2 %>%
    mutate(across(.cols = c(est, ci.lb, ci.ub), # apply an anonymous fn
                  .fns = function(x){exp(x)-0.2} )) %>% # to back-transform
    select(-se, -tval, -pval) %>%
    mutate(ci.ub = ifelse(ci.ub >1, 1, ci.ub))
write.csv(overall_trait_mod_est_addh2, file="Model Estimates/overall_trait_mod_est_addh2.csv", row.names=F)

overall_trait_mod_est_addh2

trait	est	se	tval	pval	ci.lb	ci.ub
gene expression	-1.0679538	0.2181622	-4.8952284	0.0000047	-1.5017932	-0.6341143
photochemistry	-0.7171386	0.1516333	-4.7294265	0.0000090	-1.0186781	-0.4155992
growth	-0.6915615	0.1314874	-5.2595279	0.0000011	-0.9530386	-0.4300845
nutrient content	-0.5497992	0.1632611	-3.3676060	0.0011461	-0.8744619	-0.2251365
bleaching	-0.4680659	0.1619014	-2.8910552	0.0048861	-0.7900246	-0.1461072
symbiont community	-0.2919001	0.2426445	-1.2029948	0.2323581	-0.7744254	0.1906251
morphology	-0.2747410	0.2325601	-1.1813763	0.2407872	-0.7372123	0.1877303
immune response	-0.1462565	0.1930812	-0.7574869	0.4508779	-0.5302197	0.2377067
survival	0.0611561	0.1294393	0.4724693	0.6378171	-0.1962482	0.3185604
h2narrow	-0.2413887	0.1686414	-1.4313724	0.1560337	-0.5767506	0.0939733

overall_trait_mod_est_addstage <- update(final.model, .~trait + stage -1, ) %>% # remove intercept
    summary %>% coef %>% # extract model coefficients + CI
    rownames_to_column("trait") %>%
    mutate(trait = str_remove(trait, "trait")) %>%
    rename(est = estimate)
overall_trait_mod_est_addstage_backtrans <- overall_trait_mod_est_addstage %>%
    mutate(across(.cols = c(est, ci.lb, ci.ub), # apply an anonymous fn
                  .fns = function(x){exp(x)-0.2} )) %>% # to back-transform
    select(-se, -tval, -pval) %>%
    mutate(ci.ub = ifelse(ci.ub >1, 1, ci.ub))
write.csv(overall_trait_mod_est_addstage, file="Model Estimates/overall_trait_mod_est_addstage.csv", row.names=F)

overall_trait_mod_est_addstage

trait	est	se	tval	pval	ci.lb	ci.ub
gene expression	-1.1527253	0.2551335	-4.5181250	0.0000205	-1.6601756	-0.6452749
photochemistry	-0.7910471	0.2069458	-3.8224845	0.0002545	-1.2026540	-0.3794402
growth	-0.7431715	0.1956972	-3.7975572	0.0002772	-1.1324054	-0.3539375
nutrient content	-0.6146672	0.2029390	-3.0288278	0.0032715	-1.0183047	-0.2110297
bleaching	-0.5705276	0.2016775	-2.8289105	0.0058546	-0.9716560	-0.1693992
symbiont community	-0.4511896	0.2367748	-1.9055640	0.0601691	-0.9221252	0.0197460
morphology	-0.3127279	0.2653877	-1.1783812	0.2420120	-0.8405733	0.2151176
immune response	-0.2142387	0.2328888	-0.9199184	0.3602813	-0.6774452	0.2489678
survival	-0.0049826	0.1668999	-0.0298537	0.9762554	-0.3369398	0.3269746
stagejuvenile	-0.2824011	0.2329974	-1.2120352	0.2289380	-0.7458235	0.1810214
stageadult	0.1001852	0.2008286	0.4988592	0.6191981	-0.2992547	0.4996251

These estimates are used to calculate effect sizes in separate excel spreadsheets (found at my github repo: ecolology/heritabilitymeta

Sub-analysis A: ~ trait x life stage + h2 (dat_s subset)

Final model output

# Data subset:
dat_s <- data %>% 
    filter(!(trait %in% c("photochemistry", "immune response"))) %>%
    mutate(trait = fct_drop(trait))
# table.plot("trait", "stage", Data=dat_s)

final.model.A <- rma(yi=val.log ~ trait*stage + h2, 
                     vi=sv.log, data=dat_s, method="REML", test="knha")
final.model.A


Mixed-Effects Model (k = 74; tau^2 estimator: REML)

tau^2 (estimated amount of residual heterogeneity):     0.0382 (SE = 0.0173)
tau (square root of estimated tau^2 value):             0.1953
I^2 (residual heterogeneity / unaccounted variability): 46.98%
H^2 (unaccounted variability / sampling variability):   1.89
R^2 (amount of heterogeneity accounted for):            78.27%

Test for Residual Heterogeneity:
QE(df = 55) = 96.0189, p-val = 0.0005

Test of Moderators (coefficients 2:19):
F(df1 = 18, df2 = 55) = 7.8780, p-val < .0001

Model Results:

                                       estimate      se     tval    pval 
intrcpt                                 -0.2731  0.5329  -0.5126  0.6103 
traitgrowth                              0.4395  0.1690   2.6012  0.0119 
traitnutrient content                    0.0561  0.5609   0.1001  0.9206 
traitbleaching                          -0.3175  0.5803  -0.5470  0.5866 
traitsymbiont community                  0.1795  0.5609   0.3200  0.7501 
traitmorphology                          0.1718  0.5757   0.2984  0.7665 
traitsurvival                            0.2574  0.5368   0.4794  0.6335 
stagejuvenile                            0.1620  0.1729   0.9374  0.3526 
stageadult                              -0.6607  0.5544  -1.1917  0.2385 
h2narrow                                -0.2716  0.1195  -2.2720  0.0270 
traitgrowth:stagejuvenile               -1.4359  0.5894  -2.4364  0.0181 
traitbleaching:stagejuvenile            -0.6695  0.3533  -1.8951  0.0633 
traitsymbiont community:stagejuvenile   -0.5102  0.6872  -0.7424  0.4610 
traitmorphology:stagejuvenile           -0.4727  0.3466  -1.3640  0.1781 
traitnutrient content:stageadult        -0.1588  0.5932  -0.2678  0.7899 
traitbleaching:stageadult                1.1584  0.6216   1.8636  0.0677 
traitsymbiont community:stageadult       0.6817  0.6128   1.1125  0.2708 
traitmorphology:stageadult               0.3645  0.6503   0.5605  0.5774 
traitsurvival:stageadult                 0.3266  0.5682   0.5748  0.5678 
                                         ci.lb    ci.ub 
intrcpt                                -1.3410   0.7947    
traitgrowth                             0.1009   0.7781  * 
traitnutrient content                  -1.0679   1.1802    
traitbleaching                         -1.4805   0.8456    
traitsymbiont community                -0.9446   1.3036    
traitmorphology                        -0.9819   1.3255    
traitsurvival                          -0.8184   1.3331    
stagejuvenile                          -0.1844   0.5084    
stageadult                             -1.7718   0.4504    
h2narrow                               -0.5112  -0.0320  * 
traitgrowth:stagejuvenile              -2.6170  -0.2548  * 
traitbleaching:stagejuvenile           -1.3776   0.0385  . 
traitsymbiont community:stagejuvenile  -1.8873   0.8670    
traitmorphology:stagejuvenile          -1.1672   0.2218    
traitnutrient content:stageadult       -1.3476   1.0299    
traitbleaching:stageadult              -0.0873   2.4042  . 
traitsymbiont community:stageadult     -0.5463   1.9098    
traitmorphology:stageadult             -0.9388   1.6678    
traitsurvival:stageadult               -0.8121   1.4653    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model estimates ~ trait type x life stage + heritability type:

subA_mod_est <- dat_s %>% 
    group_by(trait, h2, stage) %>%
    count() %>% ungroup %>%
    arrange(h2, stage, trait) %>%
    mutate(est = NA_real_, se = NA_real_)

for (i in 1:nrow(subA_mod_est)) {
    row.to.fill <- slice(subA_mod_est, i)
    new.data <- dat_s %>% 
        mutate(
            trait = fct_relevel(trait, as.character(row.to.fill$trait)),
            h2    = fct_relevel(h2, as.character(row.to.fill$h2)),
            stage = fct_relevel(stage, as.character(row.to.fill$stage)))
    mod <- new.data %>%
        rma(yi=val.log ~ trait * stage + h2, vi=sv.log,
               data=., method="REML", test = "knha")
    subA_mod_est[i,] <- 
        row.to.fill %>% mutate(est = coef(mod)["intrcpt"], se = ifelse(is.na(mod$se[1]), 0, mod$se[1]))
}
rm(row.to.fill, new.data, mod)

subA_mod_est <- subA_mod_est %>%
    mutate(ci.lb = est - se*qt(0.975, 100),
           ci.ub = est + se*qt(0.975, 100))
subA_mod_est_backtrans <- subA_mod_est %>%
    mutate(across(.cols = c(est, ci.lb, ci.ub), # apply an anonymous fn
                  .fns = function(x){exp(x)-0.2} )) %>% # to back-transform
    select(-se) %>%
    mutate(ci.ub = ifelse(ci.ub >1, 1, ci.ub))
write.csv(subA_mod_est, file="Model Estimates/subA_mod_est.csv", row.names=F)
# write.csv(subA_mod_est_backtrans, file="Model Estimates/subA_mod_est_backtrans.csv", row.names=F)

subA_mod_est

trait	h2	stage	n	est	se	ci.lb	ci.ub
nutrient content	broad	larvae	2	-0.2169756	0.1750709	-0.5643114	0.1303602
bleaching	broad	larvae	2	-0.5905857	0.2299316	-1.0467635	-0.1344080
morphology	broad	larvae	1	-0.1013378	0.2179227	-0.5336903	0.3310147
survival	broad	larvae	2	-0.0157691	0.1308402	-0.2753523	0.2438141
growth	broad	juvenile	2	-1.1074603	0.1355200	-1.3763281	-0.8385925
morphology	broad	juvenile	1	-0.4119960	0.2120161	-0.8326298	0.0086378
survival	broad	juvenile	1	0.1462689	0.1413821	-0.1342292	0.4267669
gene expression	broad	adult	8	-0.9338112	0.1530899	-1.2375372	-0.6300853
growth	broad	adult	15	-0.4942861	0.0715090	-0.6361579	-0.3524144
nutrient content	broad	adult	5	-1.0364919	0.1176608	-1.2699276	-0.8030562
bleaching	broad	adult	2	-0.0928555	0.1388115	-0.3682536	0.1825425
symbiont community	broad	adult	1	-0.0725707	0.1935641	-0.4565963	0.3114549
morphology	broad	adult	2	-0.3975158	0.2405139	-0.8746886	0.0796569
survival	broad	adult	6	-0.3498817	0.1061082	-0.5603973	-0.1393661
gene expression	narrow	larvae	1	-0.5447272	0.5192737	-1.5749513	0.4854970
symbiont community	narrow	larvae	2	-0.3652115	0.2120663	-0.7859450	0.0555220
survival	narrow	larvae	3	-0.2873720	0.1360182	-0.5572283	-0.0175157
growth	narrow	juvenile	6	-1.3790632	0.1333482	-1.6436222	-1.1145042
bleaching	narrow	juvenile	3	-1.3696826	0.1757371	-1.7183401	-1.0210252
symbiont community	narrow	juvenile	1	-0.7133499	0.6303829	-1.9640116	0.5373118
morphology	narrow	juvenile	1	-0.6835989	0.1969485	-1.0743391	-0.2928587
survival	narrow	juvenile	3	-0.1253340	0.1351653	-0.3934982	0.1428302
bleaching	narrow	adult	2	-0.3644584	0.1632971	-0.6884352	-0.0404816
morphology	narrow	adult	2	-0.6691187	0.2488572	-1.1628443	-0.1753931

These estimates are used to calculate effect sizes in separate excel spreadsheets (found at my github repo: ecolology/heritabilitymeta)

Post-hoc contrasts of life stage x trait type

# Create contrast matrices for testing hypotheses:
group <- paste0(dat_s$trait, "_", dat_s$stage)
group <- aggregate(model.matrix(final.model.A) ~ group, FUN=mean)
rownames(group) <- group$group
(group <- group[,-1])

	intrcpt	traitgrowth	traitnutrient content	traitbleaching	traitsymbiont community	traitmorphology	traitsurvival	stagejuvenile	stageadult	h2narrow	traitgrowth:stagejuvenile	traitbleaching:stagejuvenile	traitsymbiont community:stagejuvenile	traitmorphology:stagejuvenile	traitnutrient content:stageadult	traitbleaching:stageadult	traitsymbiont community:stageadult	traitmorphology:stageadult	traitsurvival:stageadult
bleaching_adult	1	0	0	1	0	0	0	0	1	0.50	0	0	0	0	0	1	0	0	0
bleaching_juvenile	1	0	0	1	0	0	0	1	0	1.00	0	1	0	0	0	0	0	0	0
bleaching_larvae	1	0	0	1	0	0	0	0	0	0.00	0	0	0	0	0	0	0	0	0
gene expression_adult	1	0	0	0	0	0	0	0	1	0.00	0	0	0	0	0	0	0	0	0
gene expression_larvae	1	0	0	0	0	0	0	0	0	1.00	0	0	0	0	0	0	0	0	0
growth_adult	1	1	0	0	0	0	0	0	1	0.00	0	0	0	0	0	0	0	0	0
growth_juvenile	1	1	0	0	0	0	0	1	0	0.75	1	0	0	0	0	0	0	0	0
morphology_adult	1	0	0	0	0	1	0	0	1	0.50	0	0	0	0	0	0	0	1	0
morphology_juvenile	1	0	0	0	0	1	0	1	0	0.50	0	0	0	1	0	0	0	0	0
morphology_larvae	1	0	0	0	0	1	0	0	0	0.00	0	0	0	0	0	0	0	0	0
nutrient content_adult	1	0	1	0	0	0	0	0	1	0.00	0	0	0	0	1	0	0	0	0
nutrient content_larvae	1	0	1	0	0	0	0	0	0	0.00	0	0	0	0	0	0	0	0	0
survival_adult	1	0	0	0	0	0	1	0	1	0.00	0	0	0	0	0	0	0	0	1
survival_juvenile	1	0	0	0	0	0	1	1	0	0.75	0	0	0	0	0	0	0	0	0
survival_larvae	1	0	0	0	0	0	1	0	0	0.60	0	0	0	0	0	0	0	0	0
symbiont community_adult	1	0	0	0	1	0	0	0	1	0.00	0	0	0	0	0	0	1	0	0
symbiont community_juvenile	1	0	0	0	1	0	0	1	0	1.00	0	0	1	0	0	0	0	0	0
symbiont community_larvae	1	0	0	0	1	0	0	0	0	1.00	0	0	0	0	0	0	0	0	0

# From CI plot, a number of potential differences:
#     H1: growth juvenile vs. growth adult (broad sense constant)
#     H2: bleaching juvenile vs. bleaching adult
#     H3: bleaching juvenile vs. bleaching larvae
#     H3: nutrient content larvae vs. nutrient content adult
cont.mat <- rbind(group["growth_juvenile",] - group["growth_adult",],
                  group["bleaching_juvenile",] - group["bleaching_adult",],
                  group["bleaching_juvenile",] - group["bleaching_larvae",],
                  group["nutrient content_larvae",] - group["nutrient content_adult",]) %>% as.matrix

# anova(final.model.Dii, L=cont.mat) # unadjusted
multcomp::glht(final.model.A,
                       linfct=cont.mat, 
                       test=Ftest()) %>%
    summary()


     Simultaneous Tests for General Linear Hypotheses

Fit: rma(yi = val.log ~ trait * stage + h2, vi = sv.log, data = dat_s, 
    method = "REML", test = "knha")

Linear Hypotheses:
                             Estimate Std. Error z value Pr(>|z|)    
growth_juvenile == 0          -0.8169     0.1427  -5.725  < 1e-04 ***
bleaching_juvenile == 0       -1.1410     0.2242  -5.089  < 1e-04 ***
bleaching_juvenile1 == 0      -0.7791     0.2894  -2.692 0.027577 *  
nutrient content_larvae == 0   0.8195     0.2109   3.885 0.000419 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

Effect sizes of \(h^2\) vs. \(H^2\):

beta = -coef(final.model.A)["h2narrow"]

## Equation for H2:
# H2 = e^beta * h2 + 0.2(e^beta - 1))
# h2 = 0.1, H2 = ...
exp(beta) * 0.1 + (0.2*exp(beta)-0.2)

 h2narrow 
0.1936197

# H2 = 0.2006651

# Another way to write it: H2 = m * h2 + b:
# m = 1.33555; b = 0.06711006
h2 <- seq(0.01,0.8,0.01)
H2 <- exp(beta) * h2 + (0.2*exp(beta)-0.2)
plot(H2, h2, type = 'l', xlim = c(0,1))

# other way: h2 = (H2 - 0.2(e^beta - 1))*e^(-beta)
# h2 = exp(-beta) * H2 - (0.2(e^beta - 1))*e^(-beta)
# h2 = 0.831158 * H2 -0.0337684
(0.2*exp(beta) - 0.2)*exp(-beta) # -intercept

  h2narrow 
0.04756862

# e.g., H2 = 0.2:
0.748755 * 0.2 -0.05024899

[1] 0.09950201

# h2 = 0.1

Sub-analysis Di: models of temperature difference (dat_tm subset)

Final and trait x temperature model results

Final model:

final.model.Di


Multivariate Meta-Analysis Model (k = 70; method: REML)

Variance Components:

            estim    sqrt  nlvls  fixed        factor 
sigma^2.1  0.0722  0.2688     12     no         study 
sigma^2.2  0.0616  0.2482     70     no  study/est.id 

Test for Residual Heterogeneity:
QE(df = 63) = 343.8923, p-val < .0001

Test of Moderators (coefficients 2:7):
F(df1 = 6, df2 = 63) = 8.0903, p-val < .0001

Model Results:

                       estimate      se     tval    pval    ci.lb    ci.ub 
intrcpt                 -0.7526  0.1600  -4.7042  <.0001  -1.0723  -0.4329  *** 
traitgrowth              0.0484  0.1488   0.3253  0.7460  -0.2490   0.3458      
traitnutrient content    0.1898  0.1764   1.0760  0.2860  -0.1627   0.5422      
traitbleaching           0.1699  0.2035   0.8350  0.4069  -0.2367   0.5765      
traitimmune response     0.5845  0.1969   2.9681  0.0042   0.1910   0.9780   ** 
traitsurvival            0.8209  0.1564   5.2501  <.0001   0.5084   1.1334  *** 
h2narrow                -0.3904  0.2275  -1.7159  0.0911  -0.8451   0.0643    . 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Trait x temp model:

traitxtemp.mod <- rma.mv(yi=val.log ~ trait*temp.s -1, 
                      V=sv.log, random = ~1|study/est.id,
                      data=dat_tm, method="REML", tdist=T, sigma2 = c(NA, NA))
traitxtemp.mod


Multivariate Meta-Analysis Model (k = 70; method: REML)

Variance Components:

            estim    sqrt  nlvls  fixed        factor 
sigma^2.1  0.0992  0.3149     12     no         study 
sigma^2.2  0.0405  0.2013     70     no  study/est.id 

Test for Residual Heterogeneity:
QE(df = 58) = 257.6998, p-val < .0001

Test of Moderators (coefficients 1:12):
F(df1 = 12, df2 = 58) = 8.0738, p-val < .0001

Model Results:

                              estimate      se     tval    pval    ci.lb 
traitphotochemistry            -0.8284  0.2222  -3.7288  0.0004  -1.2731 
traitgrowth                    -0.7841  0.1606  -4.8822  <.0001  -1.1056 
traitnutrient content          -0.6782  0.2031  -3.3399  0.0015  -1.0847 
traitbleaching                 -0.3909  0.2627  -1.4879  0.1422  -0.9168 
traitimmune response           -0.6819  0.2373  -2.8737  0.0057  -1.1568 
traitsurvival                   0.0307  0.1705   0.1801  0.8577  -0.3105 
temp.s                         -0.0079  0.1345  -0.0585  0.9535  -0.2771 
traitgrowth:temp.s             -0.0151  0.1586  -0.0951  0.9246  -0.3326 
traitnutrient content:temp.s    0.0526  0.1717   0.3062  0.7606  -0.2911 
traitbleaching:temp.s          -0.2549  0.1976  -1.2902  0.2021  -0.6505 
traitimmune response:temp.s     0.5337  0.2023   2.6382  0.0107   0.1288 
traitsurvival:temp.s           -0.0296  0.1519  -0.1952  0.8459  -0.3336 
                                ci.ub 
traitphotochemistry           -0.3837  *** 
traitgrowth                   -0.4626  *** 
traitnutrient content         -0.2717   ** 
traitbleaching                 0.1350      
traitimmune response          -0.2069   ** 
traitsurvival                  0.3719      
temp.s                         0.2613      
traitgrowth:temp.s             0.3024      
traitnutrient content:temp.s   0.3962      
traitbleaching:temp.s          0.1406      
traitimmune response:temp.s    0.9387    * 
traitsurvival:temp.s           0.2743      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Trait + h2 + temp model:

trait.temp.h2.mod <- rma.mv(yi=val.log ~ trait + h2 + temp.s -1, 
                      V=sv.log, random = ~1|study/est.id,
                      data=dat_tm, method="REML", tdist=T, sigma2 = c(NA, NA))
trait.temp.h2.mod


Multivariate Meta-Analysis Model (k = 70; method: REML)

Variance Components:

            estim    sqrt  nlvls  fixed        factor 
sigma^2.1  0.0706  0.2658     12     no         study 
sigma^2.2  0.0628  0.2505     70     no  study/est.id 

Test for Residual Heterogeneity:
QE(df = 62) = 337.3466, p-val < .0001

Test of Moderators (coefficients 1:8):
F(df1 = 8, df2 = 62) = 9.0449, p-val < .0001

Model Results:

                       estimate      se     tval    pval    ci.lb    ci.ub 
traitphotochemistry     -0.7899  0.1733  -4.5571  <.0001  -1.1364  -0.4434  *** 
traitgrowth             -0.7299  0.1500  -4.8666  <.0001  -1.0297  -0.4301  *** 
traitnutrient content   -0.6043  0.1877  -3.2201  0.0020  -0.9795  -0.2292   ** 
traitbleaching          -0.6157  0.2032  -3.0297  0.0036  -1.0220  -0.2095   ** 
traitimmune response    -0.2138  0.2159  -0.9901  0.3260  -0.6453   0.2178      
traitsurvival            0.0241  0.1607   0.1497  0.8815  -0.2971   0.3452      
h2narrow                -0.4083  0.2282  -1.7888  0.0785  -0.8645   0.0480    . 
temp.s                   0.0299  0.0527   0.5672  0.5726  -0.0755   0.1353      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Effect size with +1°C across low, medium, and high heritability:
x <- c(0.1, 0.5, 0.9)
# +1°C:
round(exp( log(x+0.2) +0.01*sqrt(1)) - 0.2, 3)/x

[1] 1.030000 1.014000 1.012222

round(exp( log(x+0.2) +0.03*sqrt(3)) - 0.2, 3)/x

[1] 1.160000 1.074000 1.065556

Effect sizes for temperature difference

# Effect size of heritability:
dat_tm_traitlevels <- levels(dat_tm$trait)
h2_est <- data.frame(trait=NULL, broad=NULL, narrow=NULL)
for(i in seq_along(dat_tm_traitlevels)) {
    final.model.Di %>%
    update(.~., data=mutate(dat_tm, trait = fct_relevel(trait, dat_tm_traitlevels[i]))) %>%
    summary %>% coef %>%
    slice(1) %>%
    pull(estimate) %>%
        {data.frame(trait=dat_tm_traitlevels[i], broad=round(exp(.) - 0.2,3), narrow=round(exp(.-0.38792173)-0.2,3))} %>%
        rbind(h2_est, .) -> h2_est
    }
h2_est <- h2_est %>%
    mutate(rel.diff = broad/narrow) %>% 
    arrange(-rel.diff)

write.csv(h2_est, file="Effect size calculations/subDi_trait&h2_mod_est_effsizes.csv", row.names=F)

# Effect size of temperature difference (use intercept to avoid calculating trait-specific slopes:
for(i in seq_along(levels(dat_tm$trait))) {
    new.mod <- final.model.Di %>%
        update(.~trait*temp.s, 
           data=mutate(dat_tm, 
                    trait = fct_relevel(trait, levels(dat_tm$trait)[i]))) %>%
        summary %>% coef %>%
        rownames_to_column("trait")
    temp.slope <- filter(new.mod, trait=="temp.s") %>% pull(estimate)
    temp_est_new <- new.mod %>%
        filter(trait=="intrcpt") %>%
        pull(estimate) %>% {data.frame(
            trait=dat_tm_traitlevels[i],
            h2_0C  = exp(. +temp.slope*sqrt(0)) - 0.2,
            h2_1C  = exp(. +temp.slope*sqrt(1)) - 0.2,
            h2_3C  = exp(. +temp.slope*sqrt(3)) - 0.2)}
    
    if (i == 1) {temp_est <- temp_est_new} else {temp_est <- rbind(temp_est, temp_est_new)}
        
}
temp_est

trait	h2_0C	h2_1C	h2_3C
photochemistry	0.2367350	0.2333105	0.2308207
growth	0.2565298	0.2461696	0.2387348
nutrient content	0.3075216	0.3307204	0.3483726
bleaching	0.4764549	0.3201151	0.2290853
immune response	0.3056713	0.6555619	1.0573017
survival	0.8311815	0.7932150	0.7663107

# write.csv(temp_est, file="Model Estimates/subDi_traitxtemp.s_mod_est.csv", row.names=F)

Calculate relative effect sizes by each trait, with an increase in manipulated temperature of 0C, +1C, and +3°C:

temp_est_effsizes <- temp_est %>%
    mutate(plus1_diff = ifelse(h2_1C>h2_0C, 
                               h2_1C/h2_0C, -1+h2_1C/h2_0C),
           plus3_diff = ifelse(h2_3C>h2_1C, 
                            h2_3C/h2_1C, -1+h2_3C/h2_1C))
temp_est_effsizes

trait	h2_0C	h2_1C	h2_3C	plus1_diff	plus3_diff
photochemistry	0.2367350	0.2333105	0.2308207	-0.0144653	-0.0106718
growth	0.2565298	0.2461696	0.2387348	-0.0403857	-0.0302019
nutrient content	0.3075216	0.3307204	0.3483726	1.0754378	1.0533751
bleaching	0.4764549	0.3201151	0.2290853	-0.3281313	-0.2843658
immune response	0.3056713	0.6555619	1.0573017	2.1446630	1.6128175
survival	0.8311815	0.7932150	0.7663107	-0.0456778	-0.0339180

write.csv(temp_est_effsizes, file = "Effect size calculations/subDi_traitxtemp.s_mod_est_effsizes.csv", row.names=F)

Main Figures

Fig. 1

# Fig. 1a
p1a <- data.full %>% 
    mutate(trait = fct_relevel(trait, trait.levels)) %>%
    group_by(study) %>%
    add_count(name="n_est") %>%
    group_by(study, trait, n_est) %>%
    count() %>%
    ungroup() %>%
    mutate(study = fct_reorder(study, (n_est))) %>% 
    ggplot(aes(y=study, x=n)) +
    geom_col(aes(fill=trait, group=trait)) +
    # scale_y_discrete(expand = expansion(add=c(0,1))) + 
    scale_x_continuous(expand = expansion(add=c(0,0.5))) +
    scale_fill_manual(values = rev(lightness(trait.colors, scalefac(1.1)))) +
    theme(legend.position = c(0.6, 0.35),
          legend.background = element_blank(),
          panel.grid = element_blank(),
          axis.title.y = element_blank()) +
    guides(fill = guide_legend(reverse = TRUE)) +
    labs(x = "Number of estimates", fill = "Trait type:")

# data for Fig. 1b
plot.data <-
    data.full %>% 
    mutate(trait = fct_relevel(trait, rev(trait.levels))) %>% 
    arrange(desc(trait), desc(val.log)) %>%
    mutate(est.id = factor(1:n())) %>%
    mutate(lwr = val.log - qt(0.975, 100)*se.log, upr = val.log + qt(0.975, 100)*se.log,
           lwr=ifelse(lwr<log(0+0.2),log(0+0.2),lwr), 
           upr=ifelse(upr>log(1.2),log(1.2),upr)) %>% 
    # select(study, species, est.id, val.log, trait, h2, stage, lwr, upr) %>%
    complete(trait, nesting(stage, h2, imputed.se), 
             fill = list(n = NA_real_, est = NA_real_, se = NA_real_)) %>%
    arrange(desc(trait), desc(val.log))

# Fig. 1b
p1b <- 
    plot.data %>%
    ggplot(data=., aes(y=as.numeric(est.id), x=val.log)) + 
    labs(x="Heritability estimate (± 95% CI)") + # 95% CI
    coord_cartesian(xlim =  log(c(-0.0006, 1)+0.2)) + # clip="off"
    scale_x_continuous(expand = expansion(mult=c(0.0005,0.0005)),
        breaks=log( c(0,0.10,0.25,0.5,0.75,1)+0.2),
        labels =      c(0,0.10,0.25,0.5,0.75,1)) +
    scale_y_reverse(breaks = 1:length(na.omit(plot.data$est.id)),
                    labels = na.omit(plot.data$est.id),
                    # expand = expansion(mult=c(0.0005,0.001)),
                    expand = expansion(add = c(0.01,0.01)),
                    sec.axis = 
                        sec_axis(~., name="Species", labels = na.omit(plot.data$species),
                                 breaks = 1:length(na.omit(plot.data$species)))) +
    facet_grid(trait~., scales = "free_y", space = "free_y") +
    theme(axis.text.y.left=element_blank(),
          axis.ticks.y.left=element_blank(),
          axis.text.y=element_text(size=rel(0.75)),
          axis.title.y=element_blank(),
          axis.title.x=element_text(hjust=0.5),
          strip.text.y = element_blank(),
          panel.spacing = unit(0.1, "lines"),
          panel.grid.major=element_blank(),
          panel.grid.minor=element_blank(),
          legend.text.align = 0, # expression causes centred legend labels
          axis.text.y.right = element_text(face = "italic")) + 
    
    # Add heritability legend
    geom_point(aes(x = val.log+100, fill=h2, shape=h2),
               color="white", size=4, stroke=1.5) + #shape=21,
    scale_shape_manual(name = "Heritability type:",
                       values = c(21, 21),
                       labels = c(expression("broad-sense, "*italic(H^2)),
                                  expression("narrow-sense, "*italic(h^2)))) +
    scale_fill_manual(name = "Heritability type:",
                      values = c(lightness("grey", scalefac(1.1)),
                               lightness("grey", scalefac(0.5))),
                      labels = c(expression("broad-sense, "*italic(H^2)),
                                  expression("narrow-sense, "*italic(h^2)))) +
    ggnewscale::new_scale("shape") + 
    ggnewscale::new_scale("fill") + 
    
    # Add trait legend + error bars
    geom_errorbar(
        aes(color=trait, xmin=lwr, xmax=upr, 
            linetype=imputed.se), width = 1) +
    scale_linetype_discrete(name = "Imputed 95% CI") +
    geom_point(aes(x=val.log+100, fill=trait), 
               color="white", size=3, stroke=1, shape=21) +
    scale_fill_discrete(name = "Trait type:") +
    scale_colour_discrete(name = "Trait type:") +
    ggnewscale::new_scale("fill") + 
    
    # Add actual plotted points:
    geom_point(aes(fill=interaction(trait, h2)), 
               color="white", size=2, stroke=1, shape=21) +
    scale_fill_manual(
        guide=FALSE,
        values = c(lightness(trait.colors, scalefac(1.2)),
                   lightness(trait.colors, scalefac(0.5))))

## Manually edit plot border to have trait colouration
g <- ggplot_gtable(ggplot_build(p1b))
panel_names <- which(grepl('panel-', g$layout$name))
colos <- trait.colors
k <- 1
for (i in panel_names) {
    j <- which(grepl('panel.border', g$grobs[[i]]$childrenOrder))
    g$grobs[[i]]$children[[j]]$gp$col <- colos[k]
    k <- k+1
}
# grid::grid.draw(g)
p1b <- g

# Save the plot
setEPS()
postscript("Figures/Fig1.eps", family="Helvetica", width=10, height=8)
plot_grid(p1a, p1b, ncol=2, axis="tblr", align="h", rel_widths=c(3,4)) # 1200 x 700
dev.off()

quartz_off_screen 
                2

# Plot the plot
plot_grid(p1a, p1b, ncol=2, axis="tblr", align="h", rel_widths=c(3,4)) # 1200 x 700

Fig. 1. Heritability estimates (N = 95) of various traits across 19 studies of reef-building corals. Colour indicates the specific trait type (hue) and heritability type (broad-sense H2 as lighter tint circles, narrow-sense h2 as darker shade). Left: Number of estimates reported in each study. Right: Point estimates of heritability and their associated 95% confidence/credible intervals (whiskers) on a logarithmic scale. The right-most axis indicates the species of coral. Heritability estimates closer to h2 = 1 indicate higher heritability and thus the potential for higher rates of adaptation within the population. Dashed lines represent heritability estimates where standard errors/confidence intervals were imputed.

Fig. 2

Next, we fit the model with different intercepts to accurately estimate all model coefficients, where data exists (a hacky way of avoiding using the model parameter estimates and variance-covariance matrix to estimate means + variances for all factor combinations).

missing.trait <- !(trait.levels %in% levels(data$trait))
new.levels <- trait.levels[!missing.trait]
new.colors <- trait.colors[!missing.trait]

trait.est <- final.model %>% update(.~trait-1) %>% {
    data.frame(trait = gsub("trait","", x=names(coef(.))), 
               est = coef(.), se = {.}$se, row.names=NULL)} %>%
    arrange(est) %>%
    mutate(trait = fct_reorder(trait, est))

p2 <-
    data %>% group_by(trait) %>%
    count() %>% ungroup %>%
    full_join(trait.est, by = c("trait")) %>%
    mutate(h2 = rep_len(c("broad", "narrow"),length(new.levels))) %>%
    mutate(trait = fct_relevel(trait, new.levels)) %>%
    mutate(lwr = est - se, upr = est + se,
           lwr=ifelse(lwr<log(0+0.2),log(0+0.2),lwr), 
           upr=ifelse(upr>log(1.2),log(1.2),upr)) %>%
    ggplot(aes(y=est, x=trait)) + 
    labs(y="Heritability estimate (±SE)", x=NULL) + # expression("(±"~sigma*")")
    coord_cartesian(ylim =  log(c(-0.0006, 1)+0.2)) +
    scale_y_continuous(
        expand = expansion(mult = c(0,0)),
        breaks=log( c(0,0.10,0.25,0.5,0.75,1)+0.2),
        labels =      c(0,0.10,0.25,0.5,0.75,1)) +
    theme(#axis.title.x=element_blank(),
        axis.text.x=element_text(color="white"),
        axis.ticks.x=element_blank(),
        panel.grid = element_blank()) +
    geom_errorbar(position=position_dodge(width=0.3), width=0.3,
                  aes(ymin=lwr, ymax=upr, col=trait)) +

    geom_point(aes(y=est, fill=trait), position=position_dodge(width=0.3), 
               color = "white", shape=21, size=3, stroke=1.5) +
    scale_fill_manual(name = "Trait type:", 
                        values = rev(new.colors)) +
    scale_color_manual(name = "Trait type:", 
                      values = rev(new.colors)) +
    geom_text(aes(y=lwr, x=trait, label = paste0(n)), color="grey", nudge_y = -0.05) +
    guides(fill = guide_legend(order=1, reverse=T), color = "none")

# Save it:
setEPS()
postscript("Figures/Fig2.eps", family="Helvetica", width=6, height=5)
p2
dev.off()

quartz_off_screen 
                2

# Plot it:
p2

Fig. 2. Heritability estimates ± S.E. for the overall dataset model. Traits are sorted along the spectrum according to their overall relative heritability, with values closer to 1 indicating more heritable traits. The number of estimates included in the meta-analysis for each trait type are indicated below each error bar in grey. The gamete compatibility trait type is excluded due to its reliance on only a single study/estimate.

Fig. 3

# # Load previously-saved model estimates:
# subA_mod_est <- read.csv("Model Estimates/subA_mod_est.csv")

missing.trait <- !(trait.levels %in% levels(dat_s$trait))
new.levels <- trait.levels[!missing.trait]
new.colors <- trait.colors[!missing.trait]

subA_mod_est2 <- subA_mod_est %>% # fill NA values to dodge groups appropriately
    complete(trait, nesting(stage, h2), fill = list(n=NA_real_, est=NA_real_, se=NA_real_))

p3 <-
    subA_mod_est2 %>%
    mutate(trait = fct_relevel(trait, new.levels),
           trait = fct_rev(trait),
           stage = fct_relevel(stage, "larvae", "juvenile")) %>%
    mutate(lwr = est - se, upr = est + se,
           # lwr=ifelse(lwr<log(0+0.2),log(0+0.2),lwr), 
           upr=ifelse(upr>log(1.2),log(1.2),upr)) %>%
    ggplot(aes(x=stage, y=est)) + 
    facet_wrap(trait ~ ., ncol = 2) +
    labs(y="Heritability estimate (±SE)", x=NULL) +
    coord_cartesian(ylim =  log(c(-0.0006, 1)+0.2)) +
    scale_y_continuous(
        expand = expansion(add = c(0,0.01)),
        breaks=log( c(0,0.10,0.25,0.5,0.75,1)+0.2),
        labels =      c(0,0.10,0.25,0.5,0.75,1)) +
    theme(panel.grid = element_blank(),
          strip.background = element_rect(fill="white", colour = "white"),
          strip.text = element_text(color="black"),
          legend.position = c(0.76, 0.1)
          ) +
    # geom_hline(data=trait.est,  # line for trait.est mean
    #          aes(color=trait, group=trait, yintercept=est),
    #          linetype="dashed", size=rel(0.25)) +
    
    
    # Add heritability legend:x=stage, y=est
    geom_point(aes(y=est, x=stage, fill=h2, shape=h2),
               position=position_dodge(width=0.5),
               color="white", size=3, stroke=1) + #shape=21,
    scale_shape_manual(name = "Heritability type:",
                       values = c(21, 21),
                       labels = c(expression("broad-sense, "*italic(H^2)),
                               expression("narrow-sense, "*italic(h^2)))) +
    scale_fill_manual(name = "Heritability type:",
                      values = c(lightness("grey", scalefac(1.1)),
                               lightness("grey", scalefac(0.5))),
                      labels = c(expression("broad-sense, "*italic(H^2)),
                               expression("narrow-sense, "*italic(h^2)))) +
    ggnewscale::new_scale("shape") +
    ggnewscale::new_scale("fill") +
    
    
    geom_errorbar(aes(x=stage, ymin=lwr, ymax=upr, col=trait, group=h2),
                  position=position_dodge(width=0.5), 
                  width=0.48) +
    geom_point(aes(x=stage, y=est, fill=interaction(trait, h2), shape=h2), 
               position=position_dodge(width=0.5), 
               color="white", size=3, stroke=1) +
    scale_colour_manual(name = "Trait type:", values = new.colors) +
    scale_fill_manual(name = "Trait type:",
                      values = c(lightness(new.colors,scalefac(1.2)), 
                        lightness(new.colors,scalefac(0.5)))) +
    # nutrient content is only H^2, so skip
    # Set up grey numbers in aesthetically-pleasing locations (using nudges):
    geom_text(aes(y=est, x=stage, label = paste0(n), group=h2),
              nudge_x = ifelse(subA_mod_est2$h2 == "narrow", 1, -1) *
                ifelse(subA_mod_est2$n>9,1.3,1) * 0.33,
              nudge_y = ifelse(subA_mod_est2$h2 == "broad" &
                                subA_mod_est2$trait == "survival" &
                                subA_mod_est2$stage == "juvenile", -0.1,0),
              color="grey", hjust=0.5) +

    scale_shape_manual(name = "Heritability type:", values = c(21, 21),
                       labels = c(expression("broad-sense, "*italic(H^2)),
                               expression("narrow-sense, "*italic(h^2)))) +
    guides(color = "none", fill = "none", shape="legend")
# p3



## Manually edit plot border to have trait colouration

# Generate the ggplot2 plot grob
g <- grid.force(ggplotGrob(p3))

# Get the names of grobs and their gPaths into a data.frame structure
grobs_df <- do.call(cbind.data.frame, grid.ls(g, print = FALSE))
# Build optimal gPaths that will be later used to identify grobs and edit them
grobs_df$gPath_full <- paste(grobs_df$gPath, grobs_df$name, sep = "::")
grobs_df$gPath_full <- gsub(pattern = "layout::", 
                            replacement = "", 
                            x = grobs_df$gPath_full, 
                            fixed = TRUE)
# Get the gPaths of the strip background grobs
panel_bd_gpath <- grobs_df$gPath_full[grepl(pattern = ".*panel\\.border*", 
                                            x = grobs_df$gPath_full)]
# Get the gPaths of the strip titles
strip_txt_gpath <- grobs_df$gPath_full[grepl(pattern = "strip.*titleGrob.*text.*", 
                                             x = grobs_df$gPath_full)]
# Create a vector of the colours:
nc=2 # number of facet columns (because panel_names are by row, need by col!)
nr <- ceiling(length(new.colors) / nc)
color_mat <- matrix(`length<-`(new.colors, nr * nc), nr, byrow = TRUE)
color_border <- color_mat %>% as.vector() %>% na.omit() %>% c()
color_text <- color_mat[,nc:1] %>% t() %>% as.vector %>% na.omit %>% rev

# Edit the grobs:
for (i in 1:length(panel_bd_gpath)){
  g <- editGrob(grob = g, gPath = panel_bd_gpath[i], gp = gpar(col = color_border[i]))
  g <- editGrob(grob = g, gPath = strip_txt_gpath[i], gp = gpar(col = color_text[i]))
}
# grid.newpage(); grid.draw(g)

p3 <- g

# Save it:
setEPS()
postscript("Figures/Fig3.eps", family="Helvetica", width=5, height=7)
grid.newpage(); grid.draw(p3)
dev.off()

quartz_off_screen 
                2

# Plot it:
grid.newpage(); grid.draw(p3)

Fig. 3. Heritability estimates ± S.E. across trait types with multiple life stages (x-axis) and different heritability types (lighter points = broad-sense heritability; darker points = narrow-sense heritability). Associated sample sizes (# of original estimates) are adjacent to each point in grey.

Fig. 4

It is difficult to plot the relationship of temperature across the complete dataset and across trait type (and heritability type), since temperature difference was only important when there was an interaction with trait type (when heritability is included in the model, this interaction is not supported). However, the interaction with trait type is only robust in the case of immune response (for which all estimates originate from a single study).

Thus, we report on the additive effect of trait type and for the subset (including trait type and heritability type). In all cases, the effect of trait type is negligible, as was previously seen.

# mod.temp.s <- 
#   rma.mv(yi=val.log ~ trait + sqrt(temp.diff),
#          V=sv.log, random = ~1|study/est.id, method="REML", tdist=T,
#          data=filter(data, !is.na(temp.diff)))

plot.temp.mod <- final.model.Di %>% update(. ~ trait + temp.s, data=dat_tm)

# Get values for plot:
temp.s.beta <- round(coef(plot.temp.mod)["temp.s"], 4)
temp.s.beta.se <- round(plot.temp.mod$se[names(coef(plot.temp.mod))=="temp.s"], 4)

The effect of sqrt(+1°C) on temperature difference on \(log(h^2 + 0.2)\) (±SE) across the entire dataset (while accounting for trait type, not heritability type) is: 0.0189 ± 0.053.

plot.temp.mod <- final.model.Di %>% update(. ~ trait + h2 + temp.s, data=dat_tm)

# Get values for plot:
temp.s.beta <- round(coef(plot.temp.mod)["temp.s"], 4)
temp.s.beta.se <- round(plot.temp.mod$se[names(coef(plot.temp.mod))=="temp.s"], 4)

Using the data subset, the effect of sqrt(+1°C) on temperature difference on \(log(h^2 + 0.2)\) (±SE) across the entire dataset (while accounting for heritability type and trait type) is: 0.0299 ± 0.0527.

Next, we manually predict the lines of best fit for a plot of temperature vs. hertability. We must do so manually since the final model of best fit is of class rma.mv, of which predict does not currently work with in the metafor package.

missing.trait <- !(trait.levels %in% levels(dat_tm$trait))
new.levels <- trait.levels[!missing.trait]
new.colors <- rev(trait.colors)[!missing.trait]

# create data.frame for equation of lines to plot:
abline_equation <- dat_tm %>% filter(!is.na(temp.s)) %>%
    mutate(h2_2 = factor(h2, labels = c("Broad-sense", "Narrow-sense"))) %>%
    group_by(h2_2) %>% # group by if h2 is included in the model!
    summarise(mean_h2 = mean(val.log), # quick hacky-way to get mean values (not strictly correct)
              mean_temp = (mean(temp.s))) %>%
    mutate(temp.s.slope = temp.s.beta, # slope from h2~trait+h2+temp.diff model
           temp.s.intercept = mean_h2 - mean_temp*temp.s.slope) # intercept = y_mean - slope * x_mean

p4 <-
    dat_tm %>%
    mutate(h2_2 = factor(h2, labels = c("Broad-sense", "Narrow-sense"))) %>%
    ggplot(aes(y=val.log, x=temp.s, col=trait)) + 
    facet_wrap(h2_2 ~ .) +
    guides(shape = "legend") +
    labs(y="Heritability estimate", 
         x="Temperature difference (+°C)") +
    coord_cartesian(ylim =  log(c(-0.0006, 1)+0.2)) +
    scale_y_continuous(
        expand = c(0.01, 0.01),
        breaks=log( c(0,0.10,0.25,0.5,0.75,1)+0.2),
        labels =    c(0,0.10,0.25,0.5,0.75,1)) +
    scale_x_continuous(
        limits =  sqrt(c(0,10)),
        expand = c(0.01,0.01),
        breaks=sqrt(c(0:10)),
        labels=0:10
    ) +
    theme(panel.grid = element_blank(),
          strip.background = element_rect(fill="white"),
          strip.text = element_text(colour = "black", size=rel(1.25))) +
    
    # Add fill legend:
    # geom_point(aes(y=val.log+5, fill=trait),
    #          color="white", shape=21, size=4, stroke=0) +
    scale_shape_manual(name = "Heritability:",
                       values = c(21, 21),
                       labels = c(expression("broad-sense, "*italic(H^2)),
                                  expression("narrow-sense, "*italic(h^2)))) +
    geom_point(aes(y=val.log, fill=trait, size=sqrt(1/sv.log)),
               color="white", stroke=0.1, alpha=1, shape=21) +
    scale_fill_manual(name = "Trait type:", values = new.colors, 
                      guide = guide_legend(order = 1, reverse = TRUE, override.aes = list(size = 3))) +
        scale_size_continuous(
            name = "Study precision\n(inverse of \nstandard error):",
            breaks=c(5,10,15,20),
            guide = guide_legend(order = 2, override.aes = list(fill = "grey")))
p4 <- p4 + geom_abline(data = abline_equation, 
                       color="black", linetype = "dashed",
                       aes(intercept = temp.s.intercept, 
                        slope = temp.s.slope)) # 750 x 350

# Save it:
setEPS()
postscript("Figures/Fig4.eps", family="Helvetica", width=8, height=4)
grid.newpage(); grid.draw(p4)
dev.off()

quartz_off_screen 
                2

p4

Fig. 4. Heritability vs. study temperature difference (treatment temperature relative to ambient/control temperature) for each trait type and heritability type, with the size of each point represents its relative precision (i.e., 1/sampling variance) and thus weighting in the analysis. Black lines indicate the estimated marginal mean effect of (square-root) temperature difference (sqrt[temperature difference] slope ± S.E. = 0.0295 ± 0.0528 log-(h2+0.2) heritability units), after accounting for trait type and heritability type effects.

Metafor summary statistics

In addition to \(\tau^2\), anova.rma() outputs the pseudo-R² statistic (Raudenbush, 2009), the residual heterogeneity or I² (Higgins and Thompson, 2002), and the unaccounted variability or H² (see Table below).

Statistic	Definition	Formula
\(\tau^2\)	amount of residual heterogeneity (not accounted for by model moderators)
\(R^2\)	percent of residual heterogeneity not accounted for in the reduced model that is in the full model with moderators	\[R^{2} = \frac{\hat{\tau}^{2}_{reduced}-\hat{\tau}^{2}_{full}}{\hat{\tau}^{2}_{reduced}}\]
\(I^2\)	percent of variance in the observed effects or outcomes that is residual heterogeneity (i.e., that is \(\tau^2\))	\[I^{2} = 100\% \times \frac{\hat{\tau}^{2}}{\hat{\tau}^{2}+\tilde{v}} \]
\(H^2\)	ratio of the total amount of variability in the effect size estimates to the amount of sampling variability	\[ H² = \frac{\hat{\tau}^{2} + \tilde{v}}{ \tilde{v}} \]

Both H² and I² are calculated using the within-study sampling variances, \(\tilde{v}\), which in turn is based on the inverse variance scaling of \(w_{i} = 1/v_{i}\) for the \(i^{th}\) study: \[\tilde{v}=\frac{(k-1)\sum w_{i}}{\left(\sum w_{i} \right)^{2}-\sum w^{2}_{i}} \]

Note that these values may not be accurate unless k is large (Lopez-Lopez et al., 2014). Also, note that they will depend on the estimator (ML or REML) being used and R² is only computed for mixed-effects model that include an intercept.

Explanatory factors

Heritability type
broad-sense \[H^{2}\]	the proportion of phenotypic variation \(\sigma^{2}_{P}\) explained by genetic effects (\(\sigma^{2}_{G}\)) including variation associated with additive (\(\sigma^{2}_{A}\)), dominance (\(\sigma^{2}_{D}\)), and epistatic effects (\(\sigma^{2}_{I}\)) \[H^{2} = \frac{\sigma^{2}_{G}}{\sigma^{2}_{P}}=\frac{\sigma^{2}_{A}+\sigma^{2}_{D}+\sigma^{2}_{I}}{\sigma^{2}_{P}}\]
narrow-sense \[h^{2}\]	the proportion of phenotypic variation \(\sigma^{2}_{P}\) explained by additive genetic effects (\(\sigma^{2}_{A}\)) \[h^{2} = \frac{\sigma^{2}_{A}}{\sigma^{2}_{P}}\]

Trait type
gene expression	up- or down-regulation of various genes involved in intracellular stress pathways
photochemistry	measures of symbiont photochemistry, chromoprotein content
growth	coral or corallite growth measures including calcification rates, buoyant weight change, larval areal expansion, linear extension, and new growth branches
nutrient content	total protein or carbohydrate content present in hosts or whole holobiont tissues
bleaching	symbiont cell densities or change in cell densities (during a treatment), bleaching index scores (proxy for symbiont cell density), and Chlorophyll A content (correlated to symbiont cell density)
morphology	static intraspecific corallite measurements and larval volumes upon birth
symbiont community	symbiont community indices (e.g., Leinster and Cobbald’s D), and proportion of Durusdinium symbionts
immune response	catalase and phenoloxidase activity within holobiont tissues
survival	measures of survival/mortality/settlement success including counts of settlement success or survival, percent survival/mortality at the end of a fixed period, larval survival through high temperatures, or differences in survival between control and temperature treatments
gamete compatibility	π-value, the percent larval contribution of various sires to various dams (excluded from the analysis due to only one estimate)

Coral life stage
larvae	estimates for free-swimming gamete or planula larvae stages (often asymbiotic) up to successful settlement
juvenile	estimates from post-settlement to sexually-mature adult
adult	estimates after sexual maturity or using adult coral nubbins

Temperature
temperature manipulation	binary factor of whether the study intentionally manipulated the temperatures higher than ambient conditions
temperature magnitude	difference in degrees Celsius between temperature treatment and ambient (control) conditions

Growth form
corymbose	finger-like corals
branching	arborescent; tree-like corals
massive	ball- or boulder-shaped corals
columnar	upward-growing cylindrical corals
encrusting	low, spreading corals often on hard, rocky substrates

R version 4.0.2 (2020-06-22)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.4

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] knitr_1.29         viridis_0.5.1      viridisLite_0.3.0  RColorBrewer_1.1-2
 [5] cowplot_1.1.0      gt_0.2.2           gridExtra_2.3      ggnewscale_0.4.3  
 [9] ggplotify_0.0.5    ggridges_0.5.2     shades_1.4.0       tidytext_0.2.5    
[13] metafor_2.4-0      Matrix_1.2-18      forcats_0.5.0      stringr_1.4.0     
[17] dplyr_1.0.2        purrr_0.3.4        readr_1.3.1        tidyr_1.1.2       
[21] tibble_3.0.3       ggplot2_3.3.3      tidyverse_1.3.0    readxl_1.3.1      

loaded via a namespace (and not attached):
 [1] nlme_3.1-149        fs_1.5.0            lubridate_1.7.9    
 [4] httr_1.4.2          SnowballC_0.7.0     tools_4.0.2        
 [7] backports_1.1.10    R6_2.4.1            DBI_1.1.0          
[10] colorspace_1.4-1    withr_2.2.0         tidyselect_1.1.0   
[13] compiler_4.0.2      cli_2.0.2           rvest_0.3.6        
[16] xml2_1.3.2          sandwich_2.5-1      labeling_0.3       
[19] sass_0.2.0          scales_1.1.1        checkmate_2.0.0    
[22] mvtnorm_1.1-1       commonmark_1.7      digest_0.6.25      
[25] rmarkdown_2.5       pkgconfig_2.0.3     htmltools_0.5.0    
[28] dbplyr_1.4.4        highr_0.8           rlang_0.4.7        
[31] rstudioapi_0.11     gridGraphics_0.5-0  generics_0.0.2     
[34] farver_2.0.3        zoo_1.8-8           jsonlite_1.7.1     
[37] tokenizers_0.2.1    magrittr_1.5        Rcpp_1.0.5         
[40] munsell_0.5.0       fansi_0.4.1         lifecycle_0.2.0    
[43] stringi_1.5.3       multcomp_1.4-13     yaml_2.2.1         
[46] MASS_7.3-53         plyr_1.8.6          blob_1.2.1         
[49] crayon_1.3.4        lattice_0.20-41     splines_4.0.2      
[52] haven_2.3.1         hms_0.5.3           pillar_1.4.6       
[55] codetools_0.2-16    reprex_0.3.0        glue_1.4.2         
[58] evaluate_0.14       BiocManager_1.30.10 modelr_0.1.8       
[61] vctrs_0.3.4         cellranger_1.1.0    gtable_0.3.0       
[64] assertthat_0.2.1    xfun_0.17           broom_0.7.0        
[67] janeaustenr_0.1.5   survival_3.2-3      rvcheck_0.1.8      
[70] TH.data_1.0-10      ellipsis_0.3.1

Coral Heritability Meta-analysis

Supplementary Code Documentation C: Supplementary Tables and Figures

Kevin R Bairos-Novak

Document last run on 2021-08-06