Code
library(ggplot2)
library(here)
library(lme4)
library(lmerTest)
library(glmmTMB)
library(car)
library(tidyverse)
library(emmeans)
library(asreml)
library(metafor)
library(gridExtra)C.elegans - re-analysis
1 Data wrangling
Restructuring the data to calculate fitness difference (effect size).
First load the data:
Code
alldata <- read.table(here("Data", "all.csv"), sep = ";", head = T, stringsAsFactors = T)
ancdata <- read.table(here("Data", "anc.csv"), sep = ";", head = T, stringsAsFactors = T)
gens <- read.table(here("Data", "generations.csv"), sep = ";", head = T, stringsAsFactors = T)Summarise the non-ancestral data to have means in each of 4 repeats:
Code
alldata_2 <- alldata %>%
filter(anc_pop != "X") %>%
group_by(block_unique) %>%
summarise(mean_Nstart = mean(propN_start),
mean_N1 = mean(propN_F1),
SD_Nstart = sd(propN_start),
SD_N1 = sd(propN_F1),
population = unique(population),
isoline = unique(isoline),
temperature = unique(temperature),
repr.type = unique(repr.type),
anc_pop = unique(anc_pop),
.groups = "drop")Summarise the ancestral data to have means in each of 4 repeats:
Code
ancdata_2 <- ancdata %>%
group_by(block_unique) %>%
summarise(mean_Nstart_anc = mean(propN_start),
mean_N1_anc = mean(propN_F1),
SD_Nstart_anc = sd(propN_start),
SD_N1_anc = sd(propN_F1))Merge actual data and ancestral data:
Code
alldata_2 <- alldata_2 %>%
left_join(ancdata_2, by = c("anc_pop" = "block_unique")) %>%
select(-anc_pop)
glimpse(alldata_2)Rows: 80
Columns: 13
$ block_unique <fct> I.I_E02_Iz6, I.I_E05_Iz6, I.I_E06_Iz6, I.I_E08_Iz6, I.…
$ mean_Nstart <dbl> 0.5818966, 0.5750916, 0.7437037, 0.6898734, 0.6875000,…
$ mean_N1 <dbl> 0.3776727, 0.4294794, 0.5099959, 0.3424359, 0.4647014,…
$ SD_Nstart <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ SD_N1 <dbl> 0.05545831, 0.05415434, 0.02281169, 0.02699229, 0.0166…
$ population <fct> E02, E05, E06, E08, E09, E73, E74, E75, K11, K32, K60,…
$ isoline <fct> Iz6, Iz6, Iz6, Iz6, Iz6, Iz6, Iz6, Iz6, Iz6, Iz6, Iz6,…
$ temperature <fct> 24C, 24C, 24C, 24C, 24C, 24C, 24C, 24C, 20C, 20C, 20C,…
$ repr.type <fct> FOG, FOG, FOG, FOG, FOG, WT, WT, WT, WT, WT, FOG, FOG,…
$ mean_Nstart_anc <dbl> 0.5939742, 0.5939742, 0.5939742, 0.5939742, 0.5939742,…
$ mean_N1_anc <dbl> 0.3680913, 0.3680913, 0.3680913, 0.3680913, 0.3680913,…
$ SD_Nstart_anc <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ SD_N1_anc <dbl> 0.025064967, 0.025064967, 0.025064967, 0.025064967, 0.…2 Effect sizes
We will use two types of fitness measures (for now without explicit errors): proportion difference and (natural log) of proportion ratio.
Code
alldata_2 <- alldata_2 %>%
mutate(fitness_ee_1 = mean_N1 - mean_Nstart,
fitness_ee_2 = mean_N1/mean_Nstart,
fitness_ee_3 = log(mean_N1/mean_Nstart),
fitness_anc_1 = mean_N1_anc - mean_Nstart_anc,
fitness_anc_2 = mean_N1_anc/mean_Nstart_anc,
fitness_anc_3 = log(mean_N1_anc/mean_Nstart_anc),
repr.type = relevel(repr.type, ref = "WT"))Let’s see the properties of these measures:
Code
ggplot(alldata_2) +
geom_density(mapping = aes(x = fitness_ee_1), colour = "brown") +
geom_density(mapping = aes(x = fitness_ee_2), colour = "red") +
geom_density(mapping = aes(x = fitness_ee_3), colour = "orange") +
theme_classic()
Code
ggplot(alldata_2) +
geom_density(mapping = aes(x = fitness_anc_1), colour = "darkblue") +
geom_density(mapping = aes(x = fitness_anc_2), colour = "blue") +
geom_density(mapping = aes(x = fitness_anc_3), colour = "turquoise") +
theme_classic()
As expected - the log(proportion_1/proportion_0) has much better distributional properties.
Let’s turn it into a simple effect size (fitness difference, for now without any sampling variance calculation):
Code
alldata_2 <- alldata_2 %>%
mutate(d_1 = fitness_ee_1 - fitness_anc_1,
d_2 = fitness_ee_2 - fitness_anc_2,
d_3 = fitness_ee_3 - fitness_anc_3)3 Analysis using lmer()
Code
model1 <- lmer(d_1 ~ isoline * temperature * repr.type + (1|population), data = alldata_2)boundary (singular) fit: see help('isSingular')Code
model2 <- lmer(d_2 ~ isoline * temperature * repr.type + (1|population), data = alldata_2)boundary (singular) fit: see help('isSingular')Code
model3 <- lmer(d_3 ~ isoline * temperature * repr.type + (1|population), data = alldata_2)boundary (singular) fit: see help('isSingular')Code
summary(model1)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: d_1 ~ isoline * temperature * repr.type + (1 | population)
Data: alldata_2
REML criterion at convergence: -144.9
Scaled residuals:
Min 1Q Median 3Q Max
-2.00316 -0.59848 0.05253 0.58230 2.86831
Random effects:
Groups Name Variance Std.Dev.
population (Intercept) 0.000000 0.00000
Residual 0.005085 0.07131
Number of obs: 80, groups: population, 35
Fixed effects:
Estimate Std. Error df t value
(Intercept) -0.04438 0.02521 68.00000 -1.760
isolineIz8 0.03813 0.03851 68.00000 0.990
isolineIz9 0.18667 0.04367 68.00000 4.275
temperature24C 0.03371 0.03382 68.00000 0.997
repr.typeFOG 0.06055 0.03851 68.00000 1.572
isolineIz8:temperature24C -0.05396 0.05712 68.00000 -0.945
isolineIz9:temperature24C -0.14368 0.05712 68.00000 -2.515
isolineIz8:repr.typeFOG -0.08151 0.05446 68.00000 -1.497
isolineIz9:repr.typeFOG -0.06854 0.06345 68.00000 -1.080
temperature24C:repr.typeFOG -0.05268 0.04806 68.00000 -1.096
isolineIz8:temperature24C:repr.typeFOG 0.16047 0.07688 68.00000 2.087
isolineIz9:temperature24C:repr.typeFOG 0.22186 0.09078 68.00000 2.444
Pr(>|t|)
(Intercept) 0.0828 .
isolineIz8 0.3256
isolineIz9 6.11e-05 ***
temperature24C 0.3225
repr.typeFOG 0.1205
isolineIz8:temperature24C 0.3482
isolineIz9:temperature24C 0.0143 *
isolineIz8:repr.typeFOG 0.1391
isolineIz9:repr.typeFOG 0.2838
temperature24C:repr.typeFOG 0.2769
isolineIz8:temperature24C:repr.typeFOG 0.0406 *
isolineIz9:temperature24C:repr.typeFOG 0.0171 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) islnI8 islnI9 tmp24C rp.FOG isI8:24C isI9:24C iI8:.F iI9:.F
isolineIz8 -0.655
isolineIz9 -0.577 0.378
tempertr24C -0.745 0.488 0.430
repr.typFOG -0.655 0.429 0.378 0.488
islnIz8:24C 0.441 -0.674 -0.255 -0.592 -0.289
islnIz9:24C 0.441 -0.289 -0.764 -0.592 -0.289 0.351
islnI8:.FOG 0.463 -0.707 -0.267 -0.345 -0.707 0.477 0.204
islnI9:.FOG 0.397 -0.260 -0.688 -0.296 -0.607 0.175 0.526 0.429
tmp24C:.FOG 0.525 -0.343 -0.303 -0.704 -0.801 0.417 0.417 0.567 0.486
iI8:24C:.FO -0.328 0.501 0.189 0.440 0.501 -0.743 -0.261 -0.708 -0.304
iI9:24C:.FO -0.278 0.182 0.481 0.373 0.424 -0.221 -0.629 -0.300 -0.699
t24C:. iI8:24C:
isolineIz8
isolineIz9
tempertr24C
repr.typFOG
islnIz8:24C
islnIz9:24C
islnI8:.FOG
islnI9:.FOG
tmp24C:.FOG
iI8:24C:.FO -0.625
iI9:24C:.FO -0.529 0.331
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')Code
summary(model2)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: d_2 ~ isoline * temperature * repr.type + (1 | population)
Data: alldata_2
REML criterion at convergence: -55.5
Scaled residuals:
Min 1Q Median 3Q Max
-2.3456 -0.5324 -0.0754 0.4152 3.5639
Random effects:
Groups Name Variance Std.Dev.
population (Intercept) 0.00000 0.0000
Residual 0.01893 0.1376
Number of obs: 80, groups: population, 35
Fixed effects:
Estimate Std. Error df t value
(Intercept) -0.07398 0.04865 68.00000 -1.521
isolineIz8 0.05132 0.07431 68.00000 0.691
isolineIz9 0.50076 0.08426 68.00000 5.943
temperature24C 0.05925 0.06527 68.00000 0.908
repr.typeFOG 0.09333 0.07431 68.00000 1.256
isolineIz8:temperature24C -0.06884 0.11022 68.00000 -0.625
isolineIz9:temperature24C -0.39185 0.11022 68.00000 -3.555
isolineIz8:repr.typeFOG -0.14814 0.10509 68.00000 -1.410
isolineIz9:repr.typeFOG -0.24933 0.12243 68.00000 -2.037
temperature24C:repr.typeFOG -0.06493 0.09273 68.00000 -0.700
isolineIz8:temperature24C:repr.typeFOG 0.25625 0.14835 68.00000 1.727
isolineIz9:temperature24C:repr.typeFOG 0.51977 0.17517 68.00000 2.967
Pr(>|t|)
(Intercept) 0.132962
isolineIz8 0.492187
isolineIz9 1.07e-07 ***
temperature24C 0.367198
repr.typeFOG 0.213405
isolineIz8:temperature24C 0.534366
isolineIz9:temperature24C 0.000693 ***
isolineIz8:repr.typeFOG 0.163193
isolineIz9:repr.typeFOG 0.045586 *
temperature24C:repr.typeFOG 0.486186
isolineIz8:temperature24C:repr.typeFOG 0.088653 .
isolineIz9:temperature24C:repr.typeFOG 0.004146 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) islnI8 islnI9 tmp24C rp.FOG isI8:24C isI9:24C iI8:.F iI9:.F
isolineIz8 -0.655
isolineIz9 -0.577 0.378
tempertr24C -0.745 0.488 0.430
repr.typFOG -0.655 0.429 0.378 0.488
islnIz8:24C 0.441 -0.674 -0.255 -0.592 -0.289
islnIz9:24C 0.441 -0.289 -0.764 -0.592 -0.289 0.351
islnI8:.FOG 0.463 -0.707 -0.267 -0.345 -0.707 0.477 0.204
islnI9:.FOG 0.397 -0.260 -0.688 -0.296 -0.607 0.175 0.526 0.429
tmp24C:.FOG 0.525 -0.343 -0.303 -0.704 -0.801 0.417 0.417 0.567 0.486
iI8:24C:.FO -0.328 0.501 0.189 0.440 0.501 -0.743 -0.261 -0.708 -0.304
iI9:24C:.FO -0.278 0.182 0.481 0.373 0.424 -0.221 -0.629 -0.300 -0.699
t24C:. iI8:24C:
isolineIz8
isolineIz9
tempertr24C
repr.typFOG
islnIz8:24C
islnIz9:24C
islnI8:.FOG
islnI9:.FOG
tmp24C:.FOG
iI8:24C:.FO -0.625
iI9:24C:.FO -0.529 0.331
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')Code
summary(model3)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: d_3 ~ isoline * temperature * repr.type + (1 | population)
Data: alldata_2
REML criterion at convergence: -31
Scaled residuals:
Min 1Q Median 3Q Max
-2.4977 -0.5537 -0.1069 0.5102 3.3330
Random effects:
Groups Name Variance Std.Dev.
population (Intercept) 0.00000 0.0000
Residual 0.02717 0.1648
Number of obs: 80, groups: population, 35
Fixed effects:
Estimate Std. Error df t value
(Intercept) -0.09001 0.05827 68.00000 -1.545
isolineIz8 0.05572 0.08901 68.00000 0.626
isolineIz9 0.47604 0.10093 68.00000 4.716
temperature24C 0.07288 0.07818 68.00000 0.932
repr.typeFOG 0.12127 0.08901 68.00000 1.362
isolineIz8:temperature24C -0.07232 0.13203 68.00000 -0.548
isolineIz9:temperature24C -0.35980 0.13203 68.00000 -2.725
isolineIz8:repr.typeFOG -0.21857 0.12588 68.00000 -1.736
isolineIz9:repr.typeFOG -0.16758 0.14665 68.00000 -1.143
temperature24C:repr.typeFOG -0.07670 0.11108 68.00000 -0.690
isolineIz8:temperature24C:repr.typeFOG 0.34893 0.17771 68.00000 1.963
isolineIz9:temperature24C:repr.typeFOG 0.46898 0.20984 68.00000 2.235
Pr(>|t|)
(Intercept) 0.12708
isolineIz8 0.53342
isolineIz9 1.24e-05 ***
temperature24C 0.35452
repr.typeFOG 0.17758
isolineIz8:temperature24C 0.58566
isolineIz9:temperature24C 0.00816 **
isolineIz8:repr.typeFOG 0.08704 .
isolineIz9:repr.typeFOG 0.25715
temperature24C:repr.typeFOG 0.49225
isolineIz8:temperature24C:repr.typeFOG 0.05368 .
isolineIz9:temperature24C:repr.typeFOG 0.02871 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) islnI8 islnI9 tmp24C rp.FOG isI8:24C isI9:24C iI8:.F iI9:.F
isolineIz8 -0.655
isolineIz9 -0.577 0.378
tempertr24C -0.745 0.488 0.430
repr.typFOG -0.655 0.429 0.378 0.488
islnIz8:24C 0.441 -0.674 -0.255 -0.592 -0.289
islnIz9:24C 0.441 -0.289 -0.764 -0.592 -0.289 0.351
islnI8:.FOG 0.463 -0.707 -0.267 -0.345 -0.707 0.477 0.204
islnI9:.FOG 0.397 -0.260 -0.688 -0.296 -0.607 0.175 0.526 0.429
tmp24C:.FOG 0.525 -0.343 -0.303 -0.704 -0.801 0.417 0.417 0.567 0.486
iI8:24C:.FO -0.328 0.501 0.189 0.440 0.501 -0.743 -0.261 -0.708 -0.304
iI9:24C:.FO -0.278 0.182 0.481 0.373 0.424 -0.221 -0.629 -0.300 -0.699
t24C:. iI8:24C:
isolineIz8
isolineIz9
tempertr24C
repr.typFOG
islnIz8:24C
islnIz9:24C
islnI8:.FOG
islnI9:.FOG
tmp24C:.FOG
iI8:24C:.FO -0.625
iI9:24C:.FO -0.529 0.331
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')Code
emmip(model3, ~temperature+repr.type|isoline, CIs = T,
CIarg = list(col = "blue")) +
theme_bw()
Code
emmip(model3, ~temperature+repr.type, CIs = T,
CIarg = list(col = "blue")) +
theme_bw()NOTE: Results may be misleading due to involvement in interactions
The patterns largely agree with what we can see in the model.
4 Sampling variances and meta-analysis
The sampling variances of subsequent effect sizes are based on several relationships. The sampling variances of raw proportions are (because we do not have \(N\) for all cases - we will have to make some assumptions, probably not too bold anyways)
\[ var(p)=\frac{p(1-p)}{N}, \]
and variance of the average of 4 proportions is (from stats theory)
\[var(\hat{p}) = \frac{\sum_{i} var(p_{i})}{n^{2}} = \frac{1}{16}\sum_{i} var(p_{i}).\] Variance of a difference of two proportions \(W_{1}=\hat{p}_{2}-\hat{p}_{1}\) is (from the delta method) \(var(W_{1})=var(\hat{p}_{2})+var(\hat{p}_{1})\). We are ignoring here covariance of two proportions (which is reasonable, they should not covary in general ’casue are the result of (unpredictable) competition) - but above we also ignored covariance of 4 replicated proportion which may not be reasonable, in such case we would have to increase the estimate (maybe worth trying) assuming some correlation \(r\):
\[var(\hat{p}) = \frac{\sum_{i} var(p_{i}) + 2\sum_{i}\sum_{j<1}cov(p_{i}p_{j})}{n^{2}} = \\ \frac{1}{16}\sum_{i} var(p_{i})+2\sum_{i}\sum_{j>i}r \sqrt{var(p_{i})}\sqrt{var(p_{j})}.\] This correction can be seen as an alternative to correcting for multiple comparisons (and it actually achieves what we need - decreases Type I error in zero-ES cases to roughly 5%).
Sampling variance of a log ratio is (from delta method):
\[var(W_{2}) = var[ln(\frac{\hat{p_{2}}}{\hat{p_{1}}})] = \\ [\frac{\partial}{\partial \hat{p_{2}}}(ln\:\hat{p_{2}}-ln\:\hat{p_{1}})]^{2}\:var(p_{2}) + \\ [\frac{\partial}{\partial \hat{p_{1}}}(ln\:\hat{p}_{2}-ln\:\hat{p}_{1})]^{2}\:var(\hat{p}_{1}) = \\ \frac{1}{\hat{p_{2}}^{2}}\:var(\hat{p_{2}})+\frac{1}{\hat{p_{1}}^{2}}\:var(\hat{p_{1}}).\]
And finally - using those variances and fitness estimates - we can calculate \(d\) and it’s sampling variance (in the usual way):
\[d = \frac {{W_{x,ee}} - {W_{x,anc}}} {s_{pooled}}J,\] \[s_{pooled} = \sqrt{\frac{(n_{1}-1)\:var(W_{x,ee})+(n_{2}-1)\:var(W_{x,anc})}{n_{ee}+n_{anc}-2}},\]
\[J = 1-\frac{3}{4(n_{ee}+n_{anc}-2)-1}.\]
First - let’s calculate relevant effect sizes and sampling variances. We need to repeat the above calculations adding proportion sampling variance to original data:
Code
correction <- 2*3*0.8*0.0005
# this correction assumes that on average variance of proportion estimation here is 0.0005
# (close to actual averages of 0.0003-0.0006) and that (due to coming from the same
# realisation of breeding the 4 replicates are strongly correlated (r = 0.8)
# number 3 comes from the fact that in a 4x4 covariance matrix there are 3 correlations
alldata_2 <- alldata %>%
filter(anc_pop != "X") %>%
mutate(varNstart = (propN_start*(1-propN_start))/sum,
varN1 = (propN_F1*(1-propN_F1))/(sum)) %>%
group_by(block_unique) %>%
summarise(mean_Nstart = mean(propN_start),
mean_N1 = mean(propN_F1),
var_Nstart = (sum(varNstart)/16)+correction,
var_N1 = (sum(varN1)/16)+correction,
population = unique(population),
isoline = unique(isoline),
temperature = unique(temperature),
repr.type = unique(repr.type),
anc_pop = unique(anc_pop),
block = unique(block),
.groups = "drop")Code
ancdata_2 <- ancdata %>%
mutate(varNstart = (propN_start*(1-propN_start))/sum,
varN1 = (propN_F1*(1-propN_F1))/(sum)) %>%
group_by(block_unique) %>%
summarise(mean_Nstart_anc = mean(propN_start),
mean_N1_anc = mean(propN_F1),
var_Nstart_anc = (sum(varNstart)/16)+correction,
var_N1_anc = (sum(varN1)/16)+correction)Merge actual data and ancestral data:
Code
alldata_2 <- alldata_2 %>%
left_join(ancdata_2, by = c("anc_pop" = "block_unique")) %>%
select(-anc_pop)
glimpse(alldata_2)Rows: 80
Columns: 14
$ block_unique <fct> I.I_E02_Iz6, I.I_E05_Iz6, I.I_E06_Iz6, I.I_E08_Iz6, I.…
$ mean_Nstart <dbl> 0.5818966, 0.5750916, 0.7437037, 0.6898734, 0.6875000,…
$ mean_N1 <dbl> 0.3776727, 0.4294794, 0.5099959, 0.3424359, 0.4647014,…
$ var_Nstart <dbl> 0.002519465, 0.002639630, 0.002484285, 0.002524753, 0.…
$ var_N1 <dbl> 0.002514981, 0.002638783, 0.002510318, 0.002531016, 0.…
$ population <fct> E02, E05, E06, E08, E09, E73, E74, E75, K11, K32, K60,…
$ isoline <fct> Iz6, Iz6, Iz6, Iz6, Iz6, Iz6, Iz6, Iz6, Iz6, Iz6, Iz6,…
$ temperature <fct> 24C, 24C, 24C, 24C, 24C, 24C, 24C, 24C, 20C, 20C, 20C,…
$ repr.type <fct> FOG, FOG, FOG, FOG, FOG, WT, WT, WT, WT, WT, FOG, FOG,…
$ block <fct> I.I, I.I, I.I, I.I, I.I, I.I, I.I, I.I, I.I, I.I, I.I,…
$ mean_Nstart_anc <dbl> 0.5939742, 0.5939742, 0.5939742, 0.5939742, 0.5939742,…
$ mean_N1_anc <dbl> 0.3680913, 0.3680913, 0.3680913, 0.3680913, 0.3680913,…
$ var_Nstart_anc <dbl> 0.002550122, 0.002550122, 0.002550122, 0.002550122, 0.…
$ var_N1_anc <dbl> 0.002543473, 0.002543473, 0.002543473, 0.002543473, 0.…Code
alldata_2 <- alldata_2 %>%
mutate(fitness_ee_1 = mean_N1 - mean_Nstart,
fitness_ee_2 = log(mean_N1/mean_Nstart),
var_f_ee_1 = var_Nstart + var_N1,
var_f_ee_2 = (1/mean_N1^2)*var_N1 + (1/mean_Nstart^2)*var_Nstart,
fitness_anc_1 = mean_N1_anc - mean_Nstart_anc,
fitness_anc_2 = log(mean_N1_anc/mean_Nstart_anc),
var_f_anc_1 = var_Nstart_anc + var_N1_anc,
var_f_anc_2 = (1/mean_N1_anc^2)*var_N1_anc + (1/mean_Nstart_anc^2)*var_Nstart_anc,
repr.type = relevel(repr.type, ref = "WT"))Finally let’s calculate Hedge’s \(g\):
Code
J <- 1 - (3/(4*(4+4-2) - 1))
alldata_2 <- alldata_2 %>%
mutate(d_1 = J*(fitness_ee_1 - fitness_anc_1)/sqrt((3*var_f_ee_1 + 3*var_f_anc_1)/(4+4-2)),
d_2 = J*(fitness_ee_2 - fitness_anc_2)/sqrt((3*var_f_ee_2 + 3*var_f_anc_2)/(4+4-2)),
var_d_1 = (8/16)+(d_1^2/16),
var_d_2 = (8/16)+(d_2^2/16))Add generation numbers for each population:
Code
alldata_2 <- alldata_2 %>%
left_join(gens, by = c("population" = "population"))Check if generation time is consistently related to the effect size.
Code
ggplot(alldata_2) +
geom_point(aes(x = gen_nr, y = d_2)) +
geom_smooth(aes(x = gen_nr, y = d_2), method = "lm") +
theme_classic() + theme(text = element_text(size = 15)) +
labs(x = "Number of generations", y = "Effect size")`geom_smooth()` using formula = 'y ~ x'
5 Meta-analysis
Code
alldata_2 <- mutate(alldata_2, weight_1 = 1/var_d_1, weight_2 = 1/var_d_2)
alldata_2$esid <- as.factor(1:nrow(alldata_2))
# ES = difference in proportions
model1 <- rma.mv(d_1 ~ isoline * temperature * repr.type,
random = list(~ 1|population,
~ 1|esid),
data = alldata_2,
V = var_d_1)
summary(model1)
Multivariate Meta-Analysis Model (k = 80; method: REML)
logLik Deviance AIC BIC AICc
-81.3927 162.7855 190.7855 221.8586 198.7100
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 35 no population
sigma^2.2 0.0000 0.0001 80 no esid
Test for Residual Heterogeneity:
QE(df = 68) = 77.4421, p-val = 0.2028
Test of Moderators (coefficients 2:12):
QM(df = 11) = 42.9071, p-val < .0001
Model Results:
estimate se zval pval
intrcpt -0.5032 0.2597 -1.9371 0.0527
isolineIz8 0.4688 0.3985 1.1763 0.2395
isolineIz9 2.1836 0.4902 4.4542 <.0001
temperature24C 0.3600 0.3488 1.0320 0.3021
repr.typeFOG 0.6847 0.3924 1.7449 0.0810
isolineIz8:temperature24C -0.6400 0.5851 -1.0938 0.2740
isolineIz9:temperature24C -1.7659 0.6200 -2.8482 0.0044
isolineIz8:repr.typeFOG -0.9321 0.5593 -1.6666 0.0956
isolineIz9:repr.typeFOG -1.3909 0.7060 -1.9701 0.0488
temperature24C:repr.typeFOG -0.5422 0.4911 -1.1040 0.2696
isolineIz8:temperature24C:repr.typeFOG 1.7765 0.7880 2.2543 0.0242
isolineIz9:temperature24C:repr.typeFOG 3.0495 1.0413 2.9286 0.0034
ci.lb ci.ub
intrcpt -1.0122 0.0059 .
isolineIz8 -0.3123 1.2498
isolineIz9 1.2228 3.1445 ***
temperature24C -0.3237 1.0437
repr.typeFOG -0.0844 1.4538 .
isolineIz8:temperature24C -1.7869 0.5068
isolineIz9:temperature24C -2.9811 -0.5507 **
isolineIz8:repr.typeFOG -2.0283 0.1641 .
isolineIz9:repr.typeFOG -2.7747 -0.0072 *
temperature24C:repr.typeFOG -1.5048 0.4204
isolineIz8:temperature24C:repr.typeFOG 0.2319 3.3210 *
isolineIz9:temperature24C:repr.typeFOG 1.0086 5.0904 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
model2 <- rma.mv(d_2 ~ isoline * temperature * repr.type,
random = list(~ 1|population,
~ 1|esid),
data = alldata_2,
V = var_d_2)
summary(model2)
Multivariate Meta-Analysis Model (k = 80; method: REML)
logLik Deviance AIC BIC AICc
-82.2477 164.4954 192.4954 223.5685 200.4199
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 35 no population
sigma^2.2 0.0252 0.1587 80 no esid
Test for Residual Heterogeneity:
QE(df = 68) = 78.9367, p-val = 0.1715
Test of Moderators (coefficients 2:12):
QM(df = 11) = 46.1937, p-val < .0001
Model Results:
estimate se zval pval
intrcpt -0.5066 0.2656 -1.9076 0.0564
isolineIz8 0.4350 0.4063 1.0706 0.2843
isolineIz9 2.2553 0.5021 4.4921 <.0001
temperature24C 0.3917 0.3566 1.0983 0.2721
repr.typeFOG 0.6391 0.4025 1.5877 0.1123
isolineIz8:temperature24C -0.5159 0.5969 -0.8643 0.3874
isolineIz9:temperature24C -1.7290 0.6347 -2.7241 0.0064
isolineIz8:repr.typeFOG -1.1218 0.5737 -1.9552 0.0506
isolineIz9:repr.typeFOG -1.4045 0.7205 -1.9494 0.0512
temperature24C:repr.typeFOG -0.4300 0.5040 -0.8532 0.3935
isolineIz8:temperature24C:repr.typeFOG 1.8047 0.8071 2.2360 0.0254
isolineIz9:temperature24C:repr.typeFOG 2.8437 1.0596 2.6838 0.0073
ci.lb ci.ub
intrcpt -1.0270 0.0139 .
isolineIz8 -0.3614 1.2314
isolineIz9 1.2713 3.2394 ***
temperature24C -0.3073 1.0906
repr.typeFOG -0.1498 1.4281
isolineIz8:temperature24C -1.6857 0.6540
isolineIz9:temperature24C -2.9730 -0.4850 **
isolineIz8:repr.typeFOG -2.2463 0.0027 .
isolineIz9:repr.typeFOG -2.8166 0.0076 .
temperature24C:repr.typeFOG -1.4177 0.5578
isolineIz8:temperature24C:repr.typeFOG 0.2228 3.3865 *
isolineIz9:temperature24C:repr.typeFOG 0.7669 4.9204 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Additional model with block effect - this will become our main final model.
Code
model3 <- rma.mv(d_2 ~ isoline * temperature * repr.type + gen_nr,
random = list(~ 1|population,
~ 1|block,
~ 1|esid),
data = alldata_2,
V = var_d_2)
summary(model3)
Multivariate Meta-Analysis Model (k = 80; method: REML)
logLik Deviance AIC BIC AICc
-80.8442 161.6884 193.6884 228.9635 204.5684
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 35 no population
sigma^2.2 0.0066 0.0812 8 no block
sigma^2.3 0.0208 0.1441 80 no esid
Test for Residual Heterogeneity:
QE(df = 67) = 77.4754, p-val = 0.1792
Test of Moderators (coefficients 2:13):
QM(df = 12) = 46.0337, p-val < .0001
Model Results:
estimate se zval pval
intrcpt -1.3425 0.7554 -1.7773 0.0755
isolineIz8 0.4329 0.4108 1.0537 0.2920
isolineIz9 2.2548 0.5053 4.4623 <.0001
temperature24C 0.2378 0.3775 0.6299 0.5288
repr.typeFOG 0.1707 0.5651 0.3021 0.7626
gen_nr 0.0079 0.0067 1.1837 0.2365
isolineIz8:temperature24C -0.4080 0.6010 -0.6789 0.4972
isolineIz9:temperature24C -1.6230 0.6385 -2.5417 0.0110
isolineIz8:repr.typeFOG -1.0735 0.5732 -1.8728 0.0611
isolineIz9:repr.typeFOG -1.1036 0.7631 -1.4463 0.1481
temperature24C:repr.typeFOG -0.0178 0.6120 -0.0290 0.9768
isolineIz8:temperature24C:repr.typeFOG 1.7267 0.8069 2.1400 0.0324
isolineIz9:temperature24C:repr.typeFOG 2.5976 1.0771 2.4117 0.0159
ci.lb ci.ub
intrcpt -2.8230 0.1380 .
isolineIz8 -0.3723 1.2381
isolineIz9 1.2644 3.2452 ***
temperature24C -0.5021 0.9777
repr.typeFOG -0.9369 1.2783
gen_nr -0.0052 0.0209
isolineIz8:temperature24C -1.5859 0.7699
isolineIz9:temperature24C -2.8745 -0.3714 *
isolineIz8:repr.typeFOG -2.1970 0.0500 .
isolineIz9:repr.typeFOG -2.5992 0.3920
temperature24C:repr.typeFOG -1.2173 1.1817
isolineIz8:temperature24C:repr.typeFOG 0.1452 3.3081 *
isolineIz9:temperature24C:repr.typeFOG 0.4866 4.7087 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Forest plot of all effect sizes
Code
alldata_2$repro.temp <- interaction(alldata_2$repr.type, alldata_2$temperature)
densscale <- 0.025
iso_labs <- c("Line 6", "Line 8", "Line 9")
names(iso_labs) <- c("Iz6", "Iz8", "Iz9")
ggplot(data = arrange(alldata_2, d_2)) +
geom_density(aes(x = d_2, y = (..count..)*densscale*nrow(alldata_2)),
colour = "gray", fill = "gray95", trim = F) +
geom_vline(xintercept = 0, col = 'gray50', lty = 2) +
geom_point(size = 0.5, aes(x = d_2, y = 1:nrow(alldata_2),
colour = repro.temp)) +
geom_errorbarh(aes(y = 1:nrow(alldata_2),
xmin = d_2-1.96*sqrt(var_d_2), xmax = d_2+1.96*sqrt(var_d_2),
colour = repro.temp)) +
theme_classic() +
scale_colour_manual(values = c("blue3", "purple3", "magenta", "firebrick1")) +
labs(x = 'Cohen\'s d', y = 'Replicate ID', ) +
facet_grid(~ isoline,
labeller = labeller(isoline = iso_labs)) +
theme(strip.background = element_rect(fill = "oldlace", linewidth = 0), text = element_text(size = 15))Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(count)` instead.
Divided by populations
Code
densscale <- 0.025
iso_labs <- c("Line 6", "Line 8", "Line 9")
names(iso_labs) <- c("Iz6", "Iz8", "Iz9")
alldata_2$population <- as.factor(alldata_2$population)
# create mock variables to form a good grid (6 x 6 populations)
pops <- levels(alldata_2$population)[-(1:6)]
mockids <- as.data.frame(cbind(
pops,
rep(1:6, times = 6),
rep(1:6, each = 6)
)[1:length(pops), ])Warning in cbind(pops, rep(1:6, times = 6), rep(1:6, each = 6)): number of rows
of result is not a multiple of vector length (arg 1)Code
alldata_2 <- alldata_2 %>%
left_join(mockids, by = c("population" = "pops")) %>%
arrange(isoline, repr.type, temperature)
# plots <- list()
# i <- 1
# for (popul in alldata_2$population) {
# alldata_plot <- arrange(filter(alldata_2, population == popul), d_2)
# alldata_plot$ys <- 1:nrow(alldata_plot)
# cat(alldata_plot$ys)
# cat("\n")
# plots[[i]] <- ggplot(data = alldata_plot) +
# # geom_density(aes(x = d_2, y = (..count..)*densscale*nrow(alldata_2)),
# # colour = "gray", fill = "gray95", trim = F) +
# geom_vline(xintercept = 0, col = "gray50", lty = 2) +
# geom_point(size = 0.5, aes(
# x = d_2, y = ys
# )) +
# geom_errorbarh(aes(
# y = ys,
# xmin = d_2 - 1.96 * sqrt(var_d_2), xmax = d_2 + 1.96 * sqrt(var_d_2)
# )) +
# geom_text(label = popul, x = -3, y = 4.25, size = 20, colour = "gray") +
# theme_classic() +
# # scale_colour_manual(values = c("blue3", "purple3", "magenta", "firebrick1")) +
# labs(x = "Cohen's d", y = "Replicate", ) +
# # facet_grid(V2 ~ V3) +
# xlim(min(alldata_2$d_2 - 2 * sqrt(alldata_2$var_d_2)), max(alldata_2$d_2 + 2 * sqrt(alldata_2$var_d_2))) +
# ylim(0.5, 4.5) +
# theme(
# strip.background = element_blank(),
# strip.text = element_blank(),
# panel.background = element_rect(fill = "gray98"),
# text = element_text(size = 15)
# )
# names(plots)[i] <- popul
# i <- i + 1
# }
# nCol <- ceiling(sqrt(length(plots)))
# grid.arrange(grobs = plots, ncol = nCol)
alldata_2 <- alldata_2 %>%
group_by(population) %>%
mutate(ys = 1:n())
alldata_2 %>%
ggplot() +
geom_text(aes(label = population), x = 3, y = 4.25, size = 5, colour = "gray", hjust = 0, vjust = 1) +
geom_vline(xintercept = 0, col = "gray50", lty = 2) +
geom_point(size = 0.9, aes(
x = d_2, y = ys
)) +
geom_errorbarh(aes(
y = ys,
xmin = d_2 - 1.96 * sqrt(var_d_2), xmax = d_2 + 1.96 * sqrt(var_d_2)
), height = 0.35) +
theme_classic() +
# scale_colour_manual(values = c("blue3", "purple3", "magenta", "firebrick1")) +
labs(x = "Cohen's d", y = "Replicate", ) +
facet_grid(V2 ~ V3) +
xlim(min(alldata_2$d_2 - 2 * sqrt(alldata_2$var_d_2)), max(alldata_2$d_2 + 2 * sqrt(alldata_2$var_d_2))) +
ylim(0.5, 4.5) +
theme(
strip.background = element_blank(),
strip.text = element_blank(),
panel.background = element_rect(fill = "gray98"),
text = element_text(size = 15)
)
Predictions from metafor model.
Code
gen_av <- mean(alldata_2$gen_nr)
predictions <- predict(model3, newmods = rbind(
# i8 i9 t24 rFOG gen_av it it ir ir tr itr itr
c(0, 0, 0, 0, gen_av, 0, 0, 0, 0, 0, 0, 0), #iz6 t20 rOUT
c(1, 0, 0, 0, gen_av, 0, 0, 0, 0, 0, 0, 0), #iz8 t20 rOUT
c(0, 1, 0, 0, gen_av, 0, 0, 0, 0, 0, 0, 0), #iz9 t20 rOUT
c(0, 0, 1, 0, gen_av, 0, 0, 0, 0, 0, 0, 0), #iz6 t24 rOUT
c(0, 0, 0, 1, gen_av, 0, 0, 0, 0, 0, 0, 0), #iz6 t20 rFOG
c(1, 0, 1, 0, gen_av, 1, 0, 0, 0, 0, 0, 0), #iz8 t24 rOUT
c(0, 1, 1, 0, gen_av, 0, 1, 0, 0, 0, 0, 0), #iz9 t24 rOUT
c(1, 0, 0, 1, gen_av, 0, 0, 1, 0, 0, 0, 0), #iz8 t20 rFOG
c(0, 1, 0, 1, gen_av, 0, 0, 0, 1, 0, 0, 0), #iz9 t20 rFOG
c(0, 0, 1, 1, gen_av, 0, 0, 0, 0, 1, 0, 0), #iz6 t24 rFOG
c(1, 0, 1, 1, gen_av, 1, 0, 1, 0, 1, 1, 0), #iz8 t24 rFOG
c(0, 1, 1, 1, gen_av, 0, 1, 0, 1, 1, 0, 1) #iz9 t24 rFOG
))
predictions <- as.data.frame(predictions)
predictions$izo <- c("i6", "i8", "i9", "i6", "i6", "i8", "i9", "i8", "i9", "i6", "i8", "i9")
predictions$temp <- c("t20", "t20", "t20", "t24", "t20", "t24", "t24", "t20", "t20", "t24", "t24", "t24")
predictions$repro <- c("wt", "wt", "wt", "wt", "fog", "wt", "wt", "fog", "fog", "fog", "fog", "fog")
predictions$temprepr <- paste0(predictions$temp, " x ", predictions$repro)
iso_labs <- c("Line 6", "Line 8", "Line 9")
names(iso_labs) <- c("i6", "i8", "i9")
ggplot(
data = as.data.frame(predictions),
mapping = aes(x = temprepr, y = pred)
) +
geom_line(aes(group = izo)) +
geom_errorbar(aes(ymin = ci.lb, ymax = ci.ub), colour = "blue", width = 0.25) +
geom_point(size = 2) +
scale_x_discrete(limits = c("t20 x wt", "t24 x wt", "t20 x fog", "t24 x fog")) +
facet_wrap(~izo, scales = "free_x", labeller = labeller(izo = iso_labs)) +
theme_classic() +
theme(
text = element_text(size = 15),
axis.text.x = element_text(angle = 45, hjust = 1),
strip.background = element_rect(fill = "oldlace", linewidth = 0)
) +
labs(x = "Temperature x Reproduction type", y = "Predicted effect size")
Code
# emmip(model3, ~temperature+repr.type, CIs = T,
# CIarg = list(col = "blue")) +
# theme_bw()