NSF Red Spruce common garden reaction norms (2019-20)
Anoob Prakash
16/11/2020
suppressPackageStartupMessages({
require(tidyr)
require(lme4) # for running models
require(sjPlot) # for visualizing model interactions
require(sjmisc)
require(sjlabelled)
require(tidyverse)
require(viridis)
require(lmerTest)
require(cowplot)
})Climate data
Population_climateNA <- read.csv("C:/Dropbox/UVM/Research/Code/ClimateData/ClimateNA/Family/PCA_sprucePops.txt",
header=T, sep="\t")
colnames(Population_climateNA)[colnames(Population_climateNA)=="PC1_sel"] <- "PC1"
Family_climateNA <- read.csv("C:/Dropbox/UVM/Research/Code/ClimateData/ClimateNA/Family/PCA_spruceFams.txt",
header=T, sep="\t")
colnames(Family_climateNA)[colnames(Family_climateNA)=="PC1_sel"] <- "PC1"Trait data
budset_2019 <- read.csv("./trait_data/BudSet_2019.csv", stringsAsFactors = T)
budset_2020 <- read.csv("./trait_data/BudSet_2020.csv", stringsAsFactors = T)
budbreak_2020 <- read.csv("./trait_data/BudBreak_2020_cGDD.csv", stringsAsFactors = T)
growth_2019 <- read.csv("./trait_data/Growth_2019.csv", stringsAsFactors = T)
growth_2020 <- read.csv("./trait_data/Growth_2020.csv", stringsAsFactors = T)
Plasticity <- read.csv("./trait_data/Plasticity.csv", stringsAsFactors = T)- Reorder levels
# Garden
budset_2019$Garden <- factor(budset_2019$Garden,levels = c("Vermont", "Maryland", "North_Carolina"))
budset_2020$Garden <- factor(budset_2020$Garden,levels = c("Vermont", "Maryland", "North_Carolina"))
budbreak_2020$Garden <- factor(budbreak_2020$Garden,levels = c("Vermont", "Maryland", "North_Carolina"))
growth_2019$Garden <- factor(growth_2019$Garden,levels = c("Vermont", "Maryland", "North_Carolina"))
growth_2020$Garden <- factor(growth_2020$Garden,levels = c("Vermont", "Maryland", "North_Carolina"))
# Region
budset_2019$Region <- factor(budset_2019$Region,levels = c("Core", "Margin", "Edge"))
budset_2020$Region <- factor(budset_2020$Region,levels = c("Core", "Margin", "Edge"))
budbreak_2020$Region <- factor(budbreak_2020$Region,levels = c("Core", "Margin", "Edge"))
growth_2019$Region <- factor(growth_2019$Region,levels = c("Core", "Margin", "Edge"))
growth_2020$Region <- factor(growth_2020$Region,levels = c("Core", "Margin", "Edge"))Models
Growth 2019
## working model for growth - ****
Growth19_mod1 <- lmer(data = growth_2019, na.action = na.omit,
Growth ~ Garden * Region +
(1|Population) + (1|mBed) +
(Garden|Family)) # G x E## boundary (singular) fit: see ?isSingularGrowth19_mod1a <- lmer(data = growth_2019, na.action = na.omit,
Growth ~ Garden * Region +
(1|Population) + (1|mBed) + (1|Family))
Growth19_mod1b <- lmer(data = growth_2019, na.action = na.omit,
Growth ~ Garden * Region +
(1|Population) + (1|mBed))
# model 1 anova
anova(Growth19_mod1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Garden 75.68 37.84 2 12.32 1.0818 0.369101
## Region 2991.90 1495.95 2 64.21 42.7664 1.562e-12 ***
## Garden:Region 532.33 133.08 4 722.10 3.8046 0.004536 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# impact of GxE
anova(Growth19_mod1, Growth19_mod1a) ## refitting model(s) with ML (instead of REML)## Data: growth_2019
## Models:
## Growth19_mod1a: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | Family)
## Growth19_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Growth19_mod1a 13 30358 30442 -15166 30332
## Growth19_mod1 18 30351 30467 -15158 30315 16.499 5 0.005555 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# effect of genotype
anova(Growth19_mod1a, Growth19_mod1b)## refitting model(s) with ML (instead of REML)## Data: growth_2019
## Models:
## Growth19_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth19_mod1a: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Growth19_mod1b 12 30469 30547 -15223 30445
## Growth19_mod1a 13 30358 30442 -15166 30332 113.39 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# model 1 vs 3
anova(Growth19_mod1, Growth19_mod1b)## refitting model(s) with ML (instead of REML)## Data: growth_2019
## Models:
## Growth19_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth19_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Growth19_mod1b 12 30469 30547 -15223 30445
## Growth19_mod1 18 30351 30467 -15158 30315 129.89 6 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# the effect of genotype has higher magnitude of significance compared to the
# effect of GxE
# GxE is of higher magnitude effect when not accounting for family
anova(Growth19_mod1, Growth19_mod1b)## refitting model(s) with ML (instead of REML)## Data: growth_2019
## Models:
## Growth19_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth19_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Growth19_mod1b 12 30469 30547 -15223 30445
## Growth19_mod1 18 30351 30467 -15158 30315 129.89 6 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# LRT
anova(Growth19_mod1,Growth19_mod1a,Growth19_mod1b)## refitting model(s) with ML (instead of REML)## Data: growth_2019
## Models:
## Growth19_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth19_mod1a: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | Family)
## Growth19_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Growth19_mod1b 12 30469 30547 -15223 30445
## Growth19_mod1a 13 30358 30442 -15166 30332 113.392 1 < 2.2e-16 ***
## Growth19_mod1 18 30351 30467 -15158 30315 16.499 5 0.005555 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# AIC
AIC(Growth19_mod1,Growth19_mod1a,Growth19_mod1b)## df AIC
## Growth19_mod1 18 30344.87
## Growth19_mod1a 13 30351.38
## Growth19_mod1b 12 30462.51# effect size estimates for growth trait
plot_model(Growth19_mod1, vline.color = "red", sort.est=TRUE, show.values = T, value.offset = .3) +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
# diagnostic tests for the model
# sjPlot::plot_model(Growth19_mod1, type = "diag")
tab_model(Growth19_mod1)| Growth | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 9.53 | 7.86 – 11.20 | <0.001 |
| Garden [Maryland] | 0.03 | -2.24 – 2.31 | 0.976 |
| Garden [North_Carolina] | -0.41 | -2.69 – 1.86 | 0.721 |
| Region [Margin] | 4.74 | 3.49 – 6.00 | <0.001 |
| Region [Edge] | 3.02 | 2.05 – 3.98 | <0.001 |
|
Garden [Maryland] * Region [Margin] |
0.41 | -0.79 – 1.61 | 0.506 |
|
Garden [North_Carolina] * Region [Margin] |
-1.73 | -2.95 – -0.51 | 0.005 |
|
Garden [Maryland] * Region [Edge] |
0.29 | -0.65 – 1.23 | 0.548 |
|
Garden [North_Carolina] * Region [Edge] |
-0.96 | -1.92 – -0.00 | 0.050 |
| Random Effects | |||
| σ2 | 34.98 | ||
| τ00 Family | 3.55 | ||
| τ00 Population | 1.03 | ||
| τ00 mBed | 3.13 | ||
| τ11 Family.GardenMaryland | 0.89 | ||
| τ11 Family.GardenNorth_Carolina | 0.47 | ||
| ρ01 Family.GardenMaryland | 0.35 | ||
| ρ01 Family.GardenNorth_Carolina | -0.95 | ||
| N Population | 65 | ||
| N mBed | 15 | ||
| N Family | 334 | ||
| Observations | 4681 | ||
| Marginal R2 / Conditional R2 | 0.089 / NA | ||
Growth 2020
## working model for growth - ****
Growth_mod1 <- lmer(data = growth_2020, na.action = na.omit,
Growth ~ Garden * Region +
(1|Population) + (1|mBed) +
(Garden|Family)) # G x E## boundary (singular) fit: see ?isSingular## Warning: Model failed to converge with 1 negative eigenvalue: -2.3e+01Growth_mod1a <- lmer(data = growth_2020, na.action = na.omit,
Growth ~ Garden * Region +
(1|Population) + (1|mBed) + (1|Family))
Growth_mod1b <- lmer(data = growth_2020, na.action = na.omit,
Growth ~ Garden * Region +
(1|Population) + (1|mBed))
# model 1 anova
anova(Growth_mod1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Garden 1602.55 801.27 2 13.77 17.6768 0.0001573 ***
## Region 2830.07 1415.03 2 67.79 31.2169 2.454e-10 ***
## Garden:Region 815.96 203.99 4 368.64 4.5002 0.0014624 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# impact of GxE
anova(Growth_mod1, Growth_mod1a) ## refitting model(s) with ML (instead of REML)## Data: growth_2020
## Models:
## Growth_mod1a: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | Family)
## Growth_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Growth_mod1a 13 30161 30244 -15067 30135
## Growth_mod1 18 30164 30280 -15064 30128 6.1533 5 0.2916# effect of genotype
anova(Growth_mod1a, Growth_mod1b)## refitting model(s) with ML (instead of REML)## Data: growth_2020
## Models:
## Growth_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth_mod1a: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Growth_mod1b 12 30216 30293 -15096 30192
## Growth_mod1a 13 30161 30244 -15067 30135 57.164 1 4.009e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# model 1 vs 3
anova(Growth_mod1, Growth_mod1b)## refitting model(s) with ML (instead of REML)## Data: growth_2020
## Models:
## Growth_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Growth_mod1b 12 30216 30293 -15096 30192
## Growth_mod1 18 30164 30280 -15064 30128 63.317 6 9.51e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# the effect of genotype has higher magnitude of significance compared to the
# effect of GxE
# GxE is of higher magnitude effect when not accounting for family
anova(Growth_mod1, Growth_mod1b)## refitting model(s) with ML (instead of REML)## Data: growth_2020
## Models:
## Growth_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Growth_mod1b 12 30216 30293 -15096 30192
## Growth_mod1 18 30164 30280 -15064 30128 63.317 6 9.51e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# LRT
anova(Growth_mod1,Growth_mod1a,Growth_mod1b)## refitting model(s) with ML (instead of REML)## Data: growth_2020
## Models:
## Growth_mod1b: Growth ~ Garden * Region + (1 | Population) + (1 | mBed)
## Growth_mod1a: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | Family)
## Growth_mod1: Growth ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Growth_mod1b 12 30216 30293 -15096 30192
## Growth_mod1a 13 30161 30244 -15067 30135 57.1641 1 4.009e-14 ***
## Growth_mod1 18 30164 30280 -15064 30128 6.1533 5 0.2916
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# AIC
AIC(Growth_mod1,Growth_mod1a,Growth_mod1b)## df AIC
## Growth_mod1 18 30157.95
## Growth_mod1a 13 30154.50
## Growth_mod1b 12 30209.53# effect size estimates for growth trait
plot_model(Growth_mod1, vline.color = "red", sort.est=TRUE, show.values = T, value.offset = .3) +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
# diagnostic tests for the model
# sjPlot::plot_model(Growth19_mod1, type = "diag")
tab_model(Growth_mod1)| Growth | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 16.04 | 14.91 – 17.17 | <0.001 |
| Garden [Maryland] | 0.87 | -0.50 – 2.25 | 0.212 |
| Garden [North_Carolina] | -2.94 | -4.33 – -1.56 | <0.001 |
| Region [Margin] | 3.76 | 2.16 – 5.37 | <0.001 |
| Region [Edge] | -0.16 | -1.39 – 1.07 | 0.796 |
|
Garden [Maryland] * Region [Margin] |
3.08 | 1.56 – 4.59 | <0.001 |
|
Garden [North_Carolina] * Region [Margin] |
2.69 | 1.12 – 4.26 | 0.001 |
|
Garden [Maryland] * Region [Edge] |
0.49 | -0.71 – 1.69 | 0.420 |
|
Garden [North_Carolina] * Region [Edge] |
0.67 | -0.57 – 1.92 | 0.289 |
| Random Effects | |||
| σ2 | 45.33 | ||
| τ00 Family | 0.00 | ||
| τ00 Population | 3.10 | ||
| τ00 mBed | 0.86 | ||
| τ11 Family.GardenMaryland | 4.07 | ||
| τ11 Family.GardenNorth_Carolina | 6.30 | ||
| ρ01 Family.GardenMaryland | |||
| ρ01 Family.GardenNorth_Carolina | |||
| N Population | 65 | ||
| N mBed | 15 | ||
| N Family | 334 | ||
| Observations | 4475 | ||
| Marginal R2 / Conditional R2 | 0.136 / NA | ||
Bud Set 2019
budset_2019 <- budset_2019%>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)
## working model for growth - ****
budset19_mod1 <- lmer(data = budset_2019, na.action = na.omit,
BudSet ~ Garden * Region +
(1|Population) + (1|mBed) +
(Garden|Family)) # G x E## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues## Warning: Model failed to converge with 1 negative eigenvalue: -3.2e+01budset19_mod1a <- lmer(data = budset_2019, na.action = na.omit,
BudSet ~ Garden * Region +
(1|Population) + (1|mBed) + (1|Family))
budset19_mod1b <- lmer(data = budset_2019, na.action = na.omit,
BudSet ~ Garden * Region +
(1|Population) + (1|mBed))
# model 1 anova
anova(budset19_mod1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Garden 39470 19735.0 2 13.60 31.1200 8.408e-06 ***
## Region 39318 19658.9 2 69.07 30.9999 2.465e-10 ***
## Garden:Region 12902 3225.6 4 604.80 5.0865 0.0004895 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# impact of GxE
anova(budset19_mod1, budset19_mod1a) ## refitting model(s) with ML (instead of REML)## Data: budset_2019
## Models:
## budset19_mod1a: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | Family)
## budset19_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budset19_mod1a 13 44060 44144 -22017 44034
## budset19_mod1 18 44019 44135 -21991 43983 50.972 5 8.765e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# effect of genotype
anova(budset19_mod1a, budset19_mod1b)## refitting model(s) with ML (instead of REML)## Data: budset_2019
## Models:
## budset19_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset19_mod1a: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budset19_mod1b 12 44135 44212 -22055 44111
## budset19_mod1a 13 44060 44144 -22017 44034 77.107 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# model 1 vs 3
anova(budset19_mod1, budset19_mod1b)## refitting model(s) with ML (instead of REML)## Data: budset_2019
## Models:
## budset19_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset19_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budset19_mod1b 12 44135 44212 -22055 44111
## budset19_mod1 18 44019 44135 -21991 43983 128.08 6 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# the effect of genotype has higher magnitude of significance compared to the
# effect of GxE
# GxE is of higher magnitude effect when not accounting for family
anova(budset19_mod1, budset19_mod1b)## refitting model(s) with ML (instead of REML)## Data: budset_2019
## Models:
## budset19_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset19_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budset19_mod1b 12 44135 44212 -22055 44111
## budset19_mod1 18 44019 44135 -21991 43983 128.08 6 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# LRT
anova(budset19_mod1,budset19_mod1a,budset19_mod1b)## refitting model(s) with ML (instead of REML)## Data: budset_2019
## Models:
## budset19_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset19_mod1a: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | Family)
## budset19_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budset19_mod1b 12 44135 44212 -22055 44111
## budset19_mod1a 13 44060 44144 -22017 44034 77.107 1 < 2.2e-16 ***
## budset19_mod1 18 44019 44135 -21991 43983 50.972 5 8.765e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# AIC
AIC(budset19_mod1,budset19_mod1a,budset19_mod1b)## df AIC
## budset19_mod1 18 43989.61
## budset19_mod1a 13 44029.72
## budset19_mod1b 12 44104.62# effect size estimates for budset trait
plot_model(budset19_mod1, vline.color = "red", sort.est=TRUE, show.values = T, value.offset = .3) +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
tab_model(budset19_mod1)| BudSet | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 260.61 | 255.82 – 265.40 | <0.001 |
| Garden [Maryland] | -8.61 | -14.69 – -2.52 | 0.006 |
| Garden [North_Carolina] | -20.62 | -26.83 – -14.40 | <0.001 |
| Region [Margin] | 8.89 | 3.12 – 14.66 | 0.003 |
| Region [Edge] | 16.94 | 12.49 – 21.39 | <0.001 |
|
Garden [Maryland] * Region [Margin] |
-4.30 | -9.64 – 1.05 | 0.115 |
|
Garden [North_Carolina] * Region [Margin] |
-8.88 | -15.02 – -2.74 | 0.005 |
|
Garden [Maryland] * Region [Edge] |
-7.34 | -11.55 – -3.13 | 0.001 |
|
Garden [North_Carolina] * Region [Edge] |
-1.85 | -6.66 – 2.96 | 0.450 |
| Random Effects | |||
| σ2 | 634.16 | ||
| τ00 Family | 103.15 | ||
| τ00 Population | 20.11 | ||
| τ00 mBed | 19.75 | ||
| τ11 Family.GardenMaryland | 42.42 | ||
| τ11 Family.GardenNorth_Carolina | 100.27 | ||
| ρ01 Family.GardenMaryland | -1.00 | ||
| ρ01 Family.GardenNorth_Carolina | -0.39 | ||
| ICC | 0.16 | ||
| N Population | 65 | ||
| N mBed | 15 | ||
| N Family | 334 | ||
| Observations | 4685 | ||
| Marginal R2 / Conditional R2 | 0.143 / 0.277 | ||
Bud Set 2020
budset_2020 <- budset_2020%>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)
## working model for growth - ****
budset20_mod1 <- lmer(data = budset_2020, na.action = na.omit,
BudSet ~ Garden * Region +
(1|Population) + (1|mBed) +
(Garden|Family)) # G x E## boundary (singular) fit: see ?isSingularbudset20_mod1a <- lmer(data = budset_2020, na.action = na.omit,
BudSet ~ Garden * Region +
(1|Population) + (1|mBed) + (1|Family))
budset20_mod1b <- lmer(data = budset_2020, na.action = na.omit,
BudSet ~ Garden * Region +
(1|Population) + (1|mBed))
# model 1 anova
anova(budset20_mod1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Garden 66924 33462 2 13.29 107.4274 6.252e-09 ***
## Region 6295 3148 2 62.25 10.1050 0.0001583 ***
## Garden:Region 3817 954 4 518.61 3.0635 0.0163764 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# impact of GxE
anova(budset20_mod1, budset20_mod1a) ## refitting model(s) with ML (instead of REML)## Data: budset_2020
## Models:
## budset20_mod1a: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | Family)
## budset20_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budset20_mod1a 13 37141 37224 -18558 37115
## budset20_mod1 18 37116 37231 -18540 37080 35.178 5 1.386e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# effect of genotype
anova(budset20_mod1a, budset20_mod1b)## refitting model(s) with ML (instead of REML)## Data: budset_2020
## Models:
## budset20_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset20_mod1a: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budset20_mod1b 12 37185 37261 -18580 37161
## budset20_mod1a 13 37141 37224 -18558 37115 45.36 1 1.64e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# model 1 vs 3
anova(budset20_mod1, budset20_mod1b)## refitting model(s) with ML (instead of REML)## Data: budset_2020
## Models:
## budset20_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset20_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budset20_mod1b 12 37185 37261 -18580 37161
## budset20_mod1 18 37116 37231 -18540 37080 80.538 6 2.766e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# the effect of genotype has higher magnitude of significance compared to the
# effect of GxE
# GxE is of higher magnitude effect when not accounting for family
anova(budset20_mod1, budset20_mod1b)## refitting model(s) with ML (instead of REML)## Data: budset_2020
## Models:
## budset20_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset20_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budset20_mod1b 12 37185 37261 -18580 37161
## budset20_mod1 18 37116 37231 -18540 37080 80.538 6 2.766e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# LRT
anova(budset20_mod1,budset20_mod1a,budset20_mod1b)## refitting model(s) with ML (instead of REML)## Data: budset_2020
## Models:
## budset20_mod1b: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed)
## budset20_mod1a: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | Family)
## budset20_mod1: BudSet ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budset20_mod1b 12 37185 37261 -18580 37161
## budset20_mod1a 13 37141 37224 -18558 37115 45.360 1 1.640e-11 ***
## budset20_mod1 18 37116 37231 -18540 37080 35.178 5 1.386e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# AIC
AIC(budset20_mod1,budset20_mod1a,budset20_mod1b)## df AIC
## budset20_mod1 18 37092.72
## budset20_mod1a 13 37118.40
## budset20_mod1b 12 37161.59# effect size estimates for budset trait
plot_model(budset20_mod1, vline.color = "red", sort.est=TRUE, show.values = T, value.offset = .3) +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
tab_model(budset20_mod1)| BudSet | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 235.37 | 232.04 – 238.70 | <0.001 |
| Garden [Maryland] | -25.29 | -29.69 – -20.89 | <0.001 |
| Garden [North_Carolina] | -28.81 | -33.25 – -24.37 | <0.001 |
| Region [Margin] | 6.71 | 2.89 – 10.53 | 0.001 |
| Region [Edge] | 2.02 | -0.94 – 4.97 | 0.181 |
|
Garden [Maryland] * Region [Margin] |
0.02 | -4.04 – 4.08 | 0.991 |
|
Garden [North_Carolina] * Region [Margin] |
-2.36 | -6.50 – 1.79 | 0.266 |
|
Garden [Maryland] * Region [Edge] |
-2.15 | -5.34 – 1.04 | 0.187 |
|
Garden [North_Carolina] * Region [Edge] |
2.00 | -1.33 – 5.33 | 0.240 |
| Random Effects | |||
| σ2 | 311.48 | ||
| τ00 Family | 57.16 | ||
| τ00 Population | 4.93 | ||
| τ00 mBed | 10.02 | ||
| τ11 Family.GardenMaryland | 37.27 | ||
| τ11 Family.GardenNorth_Carolina | 38.27 | ||
| ρ01 Family.GardenMaryland | -0.97 | ||
| ρ01 Family.GardenNorth_Carolina | -0.78 | ||
| N Population | 65 | ||
| N mBed | 15 | ||
| N Family | 334 | ||
| Observations | 4282 | ||
| Marginal R2 / Conditional R2 | 0.361 / NA | ||
Bud Break 2020
- cGDD
budbreak_2020 <- budbreak_2020%>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)
## working model for bud break - ****
budbreak_mod1 <- lmer(data = budbreak_2020, na.action = na.omit,
cGDD ~ Garden * Region +
(1|Population) + (1|mBed) +
(Garden|Family)) # G x E## boundary (singular) fit: see ?isSingular## Warning: Model failed to converge with 2 negative eigenvalues: -6.2e-02 -6.7e+01budbreak_mod1a <- lmer(data = budbreak_2020, na.action = na.omit,
cGDD ~ Garden * Region +
(1|Population) + (1|mBed) + (1|Family))
budbreak_mod1b <- lmer(data = budbreak_2020, na.action = na.omit,
cGDD ~ Garden * Region +
(1|Population) + (1|mBed))
# model 1 anova
anova(budbreak_mod1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Garden 2521553 1260777 2 17.95 1104.830 < 2.2e-16 ***
## Region 122334 61167 2 90.15 53.601 4.592e-16 ***
## Garden:Region 45697 11424 4 370.29 10.011 1.056e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# impact of GxE
anova(budbreak_mod1, budbreak_mod1a) ## refitting model(s) with ML (instead of REML)## Data: budbreak_2020
## Models:
## budbreak_mod1a: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | Family)
## budbreak_mod1: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budbreak_mod1a 13 47668 47752 -23821 47642
## budbreak_mod1 18 47524 47641 -23744 47488 153.5 5 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# effect of genotype
anova(budbreak_mod1a, budbreak_mod1b)## refitting model(s) with ML (instead of REML)## Data: budbreak_2020
## Models:
## budbreak_mod1b: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed)
## budbreak_mod1a: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budbreak_mod1b 12 47825 47903 -23901 47801
## budbreak_mod1a 13 47668 47752 -23821 47642 159.29 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# model 1 vs 3
anova(budbreak_mod1, budbreak_mod1b)## refitting model(s) with ML (instead of REML)## Data: budbreak_2020
## Models:
## budbreak_mod1b: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed)
## budbreak_mod1: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budbreak_mod1b 12 47825 47903 -23901 47801
## budbreak_mod1 18 47524 47641 -23744 47488 312.79 6 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# the effect of genotype has higher magnitude of significance compared to the
# effect of GxE
# GxE is of higher magnitude effect when not accounting for family
anova(budbreak_mod1, budbreak_mod1b)## refitting model(s) with ML (instead of REML)## Data: budbreak_2020
## Models:
## budbreak_mod1b: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed)
## budbreak_mod1: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budbreak_mod1b 12 47825 47903 -23901 47801
## budbreak_mod1 18 47524 47641 -23744 47488 312.79 6 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# LRT
anova(budbreak_mod1,budbreak_mod1a,budbreak_mod1b)## refitting model(s) with ML (instead of REML)## Data: budbreak_2020
## Models:
## budbreak_mod1b: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed)
## budbreak_mod1a: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (1 | Family)
## budbreak_mod1: cGDD ~ Garden * Region + (1 | Population) + (1 | mBed) + (Garden | Family)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## budbreak_mod1b 12 47825 47903 -23901 47801
## budbreak_mod1a 13 47668 47752 -23821 47642 159.29 1 < 2.2e-16 ***
## budbreak_mod1 18 47524 47641 -23744 47488 153.50 5 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# AIC
AIC(budbreak_mod1,budbreak_mod1a,budbreak_mod1b)## df AIC
## budbreak_mod1 18 47491.26
## budbreak_mod1a 13 47635.36
## budbreak_mod1b 12 47792.33# effect size estimates for budbreak trait
plot_model(budbreak_mod1, vline.color = "red", sort.est=TRUE, show.values = T, value.offset = .3) +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
tab_model(budbreak_mod1)| cGDD | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 562.01 | 557.99 – 566.04 | <0.001 |
| Garden [Maryland] | 90.76 | 85.39 – 96.13 | <0.001 |
| Garden [North_Carolina] | 120.18 | 113.90 – 126.45 | <0.001 |
| Region [Margin] | 7.47 | 1.55 – 13.39 | 0.013 |
| Region [Edge] | 12.03 | 7.48 – 16.59 | <0.001 |
|
Garden [Maryland] * Region [Margin] |
-1.43 | -8.77 – 5.90 | 0.702 |
|
Garden [North_Carolina] * Region [Margin] |
-6.34 | -16.28 – 3.60 | 0.211 |
|
Garden [Maryland] * Region [Edge] |
10.39 | 4.63 – 16.15 | <0.001 |
|
Garden [North_Carolina] * Region [Edge] |
21.99 | 14.18 – 29.80 | <0.001 |
| Random Effects | |||
| σ2 | 1141.15 | ||
| τ00 Family | 0.00 | ||
| τ00 Population | 24.59 | ||
| τ00 mBed | 10.48 | ||
| τ11 Family.GardenMaryland | 102.53 | ||
| τ11 Family.GardenNorth_Carolina | 560.79 | ||
| ρ01 Family.GardenMaryland | |||
| ρ01 Family.GardenNorth_Carolina | |||
| N Population | 65 | ||
| N mBed | 15 | ||
| N Family | 334 | ||
| Observations | 4757 | ||
| Marginal R2 / Conditional R2 | 0.725 / NA | ||
Visualization
Growth
2019
Family plot
### Reaction norm for growth
# remove mbed to get proper reaction norms
# the presence of mbeds complicates the reaction norm otherwise
Growth19_mod_nMb <- lmer(data = growth_2019,
Growth ~ Garden * Region +
(1|Population) +
(Garden|Family))## boundary (singular) fit: see ?isSingularheightData19 <- growth_2019
heightData19$pred <- predict(Growth19_mod_nMb)
# merge climateNA data
heightData19 <- merge(x=heightData19,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
# faceting by regions
# the geographic region labels are based on the genotypic differences as well
G19 <- ggplot(heightData19,aes(Garden,Growth))+
facet_grid(.~Region)
# zero_margin
# Rxn - family
G19_A <- G19 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
xlab("") + ylab("Growth in 2019 (cm)") +
scale_color_viridis(option = "A", direction = -1) +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_x_discrete(labels=c("Vermont" = "VT", "Maryland" = "MD", "North_Carolina" = "NC")) +
theme(strip.text.x = element_text(size = 30))Region plot
# Rxn - population
growth19_pop_rxn <- G19 + geom_line(aes(y=pred,group=Family,colour=Population)) +
xlab("Garden site") + ylab("Vertical height growth (in cm)") +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
## Reaction norm for geo genetic groups/regions
# growth reaction norm
G19_B <- ggplot(heightData19,
aes(x=Garden, y=Growth, color=Region)) +
stat_summary(aes(group = Region), fun.y = mean, geom = "path", size = 1.5) +
stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1, size=1.5) +
stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 5) +
# Lines by species using grouping
stat_summary(aes(group = Region), geom = "line", fun = mean) +
ylab("") + xlab("") +
theme(legend.position="right") +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_color_manual(values = c("Core" = "goldenrod2","Margin"="green2", "Edge"="steelblue"))## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.y` is deprecated. Use `fun` instead.2020
Family plot
### Reaction norm for growth
# remove mbed to get proper reaction norms
# the presence of mbeds complicates the reaction norm otherwise
Growth20_mod_nMb <- lmer(data = growth_2020,
Growth ~ Garden * Region +
(1|Population) +
(Garden|Family))## boundary (singular) fit: see ?isSingular## Warning: Model failed to converge with 1 negative eigenvalue: -1.2e+02heightData20 <- growth_2020
heightData20$pred <- predict(Growth20_mod_nMb)
# merge climateNA data
heightData20 <- merge(x=heightData20,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
# faceting by regions
# the geographic region labels are based on the genotypic differences as well
G20 <- ggplot(heightData20,aes(Garden,Growth))+
facet_grid(.~Region)
# zero_margin
# Rxn - family
G20_C <- G20 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
xlab("") + ylab("Growth in 2020 (cm)") +
scale_color_viridis(option = "A", direction = -1) +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_x_discrete(labels=c("Vermont" = "VT", "Maryland" = "MD", "North_Carolina" = "NC")) +
theme(strip.text.x = element_text(size = 30))Region plot
# Rxn - population
growth20_pop_rxn <- G20 + geom_line(aes(y=pred,group=Family,colour=Population)) +
xlab("Garden site") + ylab("Vertical height growth (in cm)") +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
## Reaction norm for geo genetic groups/regions
# growth reaction norm
G20_D <- ggplot(heightData20,
aes(x=Garden, y=Growth, color=Region)) +
stat_summary(aes(group = Region), fun.y = mean, geom = "path", size = 1.5) +
stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1, size=1.5) +
stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 5) +
# Lines by species using grouping
stat_summary(aes(group = Region), geom = "line", fun = mean) +
ylab("") + xlab("") +
theme(legend.position="right") +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_color_manual(values = c("Core" = "goldenrod2","Margin"="green2", "Edge"="steelblue"))## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.y` is deprecated. Use `fun` instead.Combined plot for growth
legend_PC <- get_legend(G20_C)
legend_RE <- get_legend(G20_D)
leg_full <- plot_grid(legend_PC, legend_RE,
nrow = 2)
G_grid <- plot_grid(G19_A + theme(legend.position="none"),
G19_B + theme(legend.position="none"),
G20_C + theme(legend.position="none"),
G20_D + theme(legend.position="none"),
align = 'vh',
labels = c("A", "B","C", "D"),
hjust = -1,
nrow = 2,
rel_widths = c(1,.7))## Warning: Graphs cannot be vertically aligned unless the axis parameter is set.
## Placing graphs unaligned.## Warning: Graphs cannot be horizontally aligned unless the axis parameter is set.
## Placing graphs unaligned.# bud set final
plot_grid(G_grid, leg_full, ncol = 2, rel_widths = c(1,.05))
Phenology
Bud set 2019 (DOY)
Family plot
### Reaction norm for growth
# remove mbed to get proper reaction norms
# the presence of mbeds complicates the reaction norm otherwise
budset19_mod_nMb <- lmer(data = budset_2019,
BudSet ~ Garden * Region +
(1|Population) +
(Garden|Family))## boundary (singular) fit: see ?isSingular## Warning: Model failed to converge with 1 negative eigenvalue: -8.5e+00budset_2019_2 <- budset_2019
budset_2019_2$pred <- predict(budset19_mod_nMb)
# merge climateNA data
budset_2019_2 <- merge(x=budset_2019_2,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
# faceting by regions
# the geographic region labels are based on the genotypic differences as well
BS19 <- ggplot(budset_2019_2,aes(Garden,BudSet))+
facet_grid(.~Region)
# zero_margin
# Rxn - family
BS19_A <- BS19 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
xlab("") + ylab("Bud set in 2019 (DOY)") +
scale_color_viridis(option = "A", direction = -1) +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_x_discrete(labels=c("Vermont" = "VT", "Maryland" = "MD", "North_Carolina" = "NC")) +
theme(strip.text.x = element_text(size = 30))Region plot
# Rxn - population
budset19_pop_rxn <- BS19 + geom_line(aes(y=pred,group=Family,colour=Population)) +
xlab("Garden site") + ylab("Bud set in 2019 (DOY)") +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
## Reaction norm for geo genetic groups/regions
# growth reaction norm
BS19_B <- ggplot(budset_2019_2,
aes(x=Garden, y=BudSet, color=Region)) +
stat_summary(aes(group = Region), fun.y = mean, geom = "path", size = 1.5) +
stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1, size=1.5) +
stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 5) +
# Lines by species using grouping
stat_summary(aes(group = Region), geom = "line", fun = mean) +
ylab("") + xlab("") +
theme(legend.position="right") +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_color_manual(values = c("Core" = "goldenrod2","Margin"="green2", "Edge"="steelblue"))## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.y` is deprecated. Use `fun` instead.Bud set 2020 (DOY)
Family plot
### Reaction norm for growth
# remove mbed to get proper reaction norms
# the presence of mbeds complicates the reaction norm otherwise
budset20_mod_nMb <- lmer(data = budset_2020,
BudSet ~ Garden * Region +
(1|Population) +
(Garden|Family))## boundary (singular) fit: see ?isSingular## Warning: Model failed to converge with 1 negative eigenvalue: -6.8e+01budset_2020_2 <- budset_2020
budset_2020_2$pred <- predict(budset20_mod_nMb)
# merge climateNA data
budset_2020_2 <- merge(x=budset_2020_2,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
# faceting by regions
# the geographic region labels are based on the genotypic differences as well
BS20 <- ggplot(budset_2020_2,aes(Garden,BudSet))+
facet_grid(.~Region)
# zero_margin
# Rxn - family
BS20_C <- BS20 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
xlab("") + ylab("Bud set in 2020 (DOY)") +
scale_color_viridis(option = "A", direction = -1) +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_x_discrete(labels=c("Vermont" = "VT", "Maryland" = "MD", "North_Carolina" = "NC")) +
theme(strip.text.x = element_text(size = 30))Region plot
# Rxn - population
budset20_pop_rxn <- BS20 + geom_line(aes(y=pred,group=Family,colour=Population)) +
xlab("Garden site") + ylab("Bud set in 2020 (DOY)") +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
## Reaction norm for geo genetic groups/regions
# growth reaction norm
BS20_D <- ggplot(budset_2020_2,
aes(x=Garden, y=BudSet, color=Region)) +
stat_summary(aes(group = Region), fun.y = mean, geom = "path", size = 1.5) +
stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1, size=1.5) +
stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 5) +
# Lines by species using grouping
stat_summary(aes(group = Region), geom = "line", fun = mean) +
ylab("") + xlab("") +
theme(legend.position="right") +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_color_manual(values = c("Core" = "goldenrod2","Margin"="green2", "Edge"="steelblue"))## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.y` is deprecated. Use `fun` instead.Combined plot for bud set
legend_PC <- get_legend(BS20_C)
legend_RE <- get_legend(BS20_D)
leg_full <- plot_grid(legend_PC, legend_RE,
nrow = 2)
BS_grid <- plot_grid(BS19_A + theme(legend.position="none"),
BS19_B + theme(legend.position="none"),
BS20_C + theme(legend.position="none"),
BS20_D + theme(legend.position="none"),
align = 'vh',
labels = c("A", "B","C", "D"),
hjust = -1,
nrow = 2,
rel_widths = c(1,.7))## Warning: Graphs cannot be vertically aligned unless the axis parameter is set.
## Placing graphs unaligned.## Warning: Graphs cannot be horizontally aligned unless the axis parameter is set.
## Placing graphs unaligned.# bud set final
plot_grid(BS_grid, leg_full, ncol = 2, rel_widths = c(1,.05))
Bud Break 2020 (cGDD)
Family plot
### Reaction norm for growth
# remove mbed to get proper reaction norms
# the presence of mbeds complicates the reaction norm otherwise
budbreak_mod_nMb <- lmer(data = budbreak_2020,
cGDD ~ Garden * Region +
(1|Population) +
(Garden|Family))## boundary (singular) fit: see ?isSingular## Warning: Model failed to converge with 1 negative eigenvalue: -1.3e+02budbreak <- budbreak_2020[!is.na(budbreak_2020$cGDD),]
budbreak$pred <- predict(budbreak_mod_nMb)
# merge climateNA data
budbreak <- merge(x=budbreak,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
# faceting by regions
# the geographic region labels are based on the genotypic differences as well
BB20 <- ggplot(budbreak,aes(Garden,cGDD))+
facet_grid(.~Region)
# zero_margin
# Rxn - family
BB20_A <- BB20 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
xlab("") + ylab("Bud break in 2020 (cGDD)") +
scale_color_viridis(option = "A", direction = -1) +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_x_discrete(labels=c("Vermont" = "VT", "Maryland" = "MD", "North_Carolina" = "NC")) +
theme(strip.text.x = element_text(size = 30))Region plot
# Rxn - population
budbreak_pop_rxn <- BB20 + geom_line(aes(y=pred,group=Family,colour=Population)) +
xlab("Garden site") + ylab("cGDD") +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
## Reaction norm for geo genetic groups/regions
BB20_B <- ggplot(budbreak,
aes(x=Garden, y=cGDD, color=Region)) +
stat_summary(aes(group = Region), fun.y = mean, geom = "path", size = 1.5) +
stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1, size=1.5) +
stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 5) +
# Lines by species using grouping
stat_summary(aes(group = Region), geom = "line", fun = mean) +
ylab("") + xlab("") +
theme(legend.position="right") +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_color_manual(values = c("Core" = "goldenrod2","Margin"="green2", "Edge"="steelblue"))## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.y` is deprecated. Use `fun` instead.Combined plot for bud break - cGDD
legend_PC <- get_legend(BB20_A)
legend_RE <- get_legend(BB20_B)
leg_full <- plot_grid(legend_PC, legend_RE,
nrow = 2)
BB_grid <- plot_grid(BB20_A + theme(legend.position="none"),
BB20_B + theme(legend.position="none"),
align = 'vh',
labels = c("A", "B"),
hjust = -1,
nrow = 1,
rel_widths = c(1,.7))## Warning: Graphs cannot be vertically aligned unless the axis parameter is set.
## Placing graphs unaligned.## Warning: Graphs cannot be horizontally aligned unless the axis parameter is set.
## Placing graphs unaligned.# bud set final
plot_grid(BB_grid, leg_full, ncol = 2, rel_widths = c(1,.08))
Figure in text
A <- ggplot(budset_2019_2,
aes(x=Garden, y=BudSet, color=Region)) +
stat_summary(aes(group = Region), fun.y = mean, geom = "path", size = 1.5) +
stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1, size=1.5) +
stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 5) +
# Lines by species using grouping
stat_summary(aes(group = Region), geom = "line", fun = mean) +
ylab("Bud set (DOY)") + xlab("") +
theme(legend.position="right") +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_color_manual(values = c("Core" = "goldenrod2","Margin"="green2", "Edge"="steelblue"))## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.y` is deprecated. Use `fun` instead.B <- ggplot(budset_2020_2,
aes(x=Garden, y=BudSet, color=Region)) +
stat_summary(aes(group = Region), fun.y = mean, geom = "path", size = 1.5) +
stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1, size=1.5) +
stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 5) +
# Lines by species using grouping
stat_summary(aes(group = Region), geom = "line", fun = mean) +
ylab("") + xlab("") +
theme(legend.position="right") +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_color_manual(values = c("Core" = "goldenrod2","Margin"="green2", "Edge"="steelblue"))## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.y` is deprecated. Use `fun` instead.C <- ggplot(budbreak,
aes(x=Garden, y=cGDD, color=Region)) +
stat_summary(aes(group = Region), fun.y = mean, geom = "path", size = 1.5) +
stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1, size=1.5) +
stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 5) +
# Lines by species using grouping
stat_summary(aes(group = Region), geom = "line", fun = mean) +
ylab("Bud break (cGDD)") + xlab("") +
theme(legend.position="right") +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_color_manual(values = c("Core" = "goldenrod2","Margin"="green2", "Edge"="steelblue"))## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.y` is deprecated. Use `fun` instead.D <- ggplot(heightData19,
aes(x=Garden, y=Growth, color=Region)) +
stat_summary(aes(group = Region), fun.y = mean, geom = "path", size = 1.5) +
stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1, size=1.5) +
stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 5) +
# Lines by species using grouping
stat_summary(aes(group = Region), geom = "line", fun = mean) +
ylab("Growth (cm)") + xlab("2019") +
theme(legend.position="right") +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_color_manual(values = c("Core" = "goldenrod2","Margin"="green2", "Edge"="steelblue"))## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.y` is deprecated. Use `fun` instead.E <- ggplot(heightData20,
aes(x=Garden, y=Growth, color=Region)) +
stat_summary(aes(group = Region), fun.y = mean, geom = "path", size = 1.5) +
stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1, size=1.5) +
stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 5) +
# Lines by species using grouping
stat_summary(aes(group = Region), geom = "line", fun = mean) +
ylab("") + xlab("2020") +
theme(legend.position="right") +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 20),
legend.text = element_text(color = "black", size = 18)) +
scale_color_manual(values = c("Core" = "goldenrod2","Margin"="green2", "Edge"="steelblue"))## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.y` is deprecated. Use `fun` instead.legend <- get_legend(E + theme(legend.position="right"))
HA_grid <- cowplot::plot_grid( A + theme(legend.position="none"),
B + theme(legend.position="none"),
legend,
C + theme(legend.position="none"),
D + theme(legend.position="none"),
E + theme(legend.position="none"),
align = 'vh',
labels = c("A", "B","" ,"C", "D","E"),
hjust = -1,
nrow = 3
)## Warning: Graphs cannot be vertically aligned unless the axis parameter is set.
## Placing graphs unaligned.## Warning: Graphs cannot be horizontally aligned unless the axis parameter is set.
## Placing graphs unaligned.HA_grid
Supplementary
GxE at family level
S.A <- BS19 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
xlab("") + ylab("Bud set (DOY)") +
scale_color_viridis(option = "A", direction = -1) +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_x_discrete(labels=c("Vermont" = "VT", "Maryland" = "MD", "North_Carolina" = "NC")) +
theme(strip.text.x = element_text(size = 30))
S.B <- BS20 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
xlab("") + ylab("") +
scale_color_viridis(option = "A", direction = -1) +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_x_discrete(labels=c("Vermont" = "VT", "Maryland" = "MD", "North_Carolina" = "NC")) +
theme(strip.text.x = element_text(size = 30))
S.C <- BB20 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
xlab("") + ylab("Bud break (cGDD)") +
scale_color_viridis(option = "A", direction = -1) +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_x_discrete(labels=c("Vermont" = "VT", "Maryland" = "MD", "North_Carolina" = "NC")) +
theme(strip.text.x = element_text(size = 30))
S.D <- G19 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
xlab("2019") + ylab("Growth (cm)") +
scale_color_viridis(option = "A", direction = -1) +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_x_discrete(labels=c("Vermont" = "VT", "Maryland" = "MD", "North_Carolina" = "NC")) +
theme(strip.text.x = element_text(size = 30))
S.E <- G20 + geom_line(aes(y=pred,group=Family,colour=PC1)) +
xlab("2020") + ylab("") +
scale_color_viridis(option = "A", direction = -1) +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_x_discrete(labels=c("Vermont" = "VT", "Maryland" = "MD", "North_Carolina" = "NC")) +
theme(strip.text.x = element_text(size = 30))
legend <- get_legend(S.E + theme(legend.position="right"))
SA_grid <- cowplot::plot_grid( S.A + theme(legend.position="none"),
S.B + theme(legend.position="none"),
legend,
S.C + theme(legend.position="none"),
S.D + theme(legend.position="none"),
S.E + theme(legend.position="none"),
align = 'vh',
labels = c("A", "B","" ,"C", "D","E"),
hjust = -1,
nrow = 3
)## Warning: Graphs cannot be vertically aligned unless the axis parameter is set.
## Placing graphs unaligned.## Warning: Graphs cannot be horizontally aligned unless the axis parameter is set.
## Placing graphs unaligned.SA_grid
Bud break 2020 (DOY)
Family plot
### Reaction norm for growth
# remove mbed to get proper reaction norms
# the presence of mbeds complicates the reaction norm otherwise
budbreak_DOY_mod_nMb <- lmer(data = budbreak_2020,
BudBreak ~ Garden * Region +
(1|Population) +
(Garden|Family))## boundary (singular) fit: see ?isSingular## Warning: Model failed to converge with 1 negative eigenvalue: -1.1e+01budbreak <- budbreak_2020[!is.na(budbreak_2020$BudBreak),]
budbreak$pred <- predict(budbreak_DOY_mod_nMb)
# merge climateNA data
budbreak <- merge(x=budbreak,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
# faceting by regions
# the geographic region labels are based on the genotypic differences as well
BB20_DOY <- ggplot(budbreak,aes(Garden,BudBreak))+
facet_grid(.~Region)
# zero_margin
# Rxn - family
BB20_DOY_A <- BB20_DOY + geom_line(aes(y=pred,group=Family,colour=PC1)) +
xlab("") + ylab("Bud break in 2020 (DOY)") +
scale_color_viridis(option = "A", direction = -1) +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_x_discrete(labels=c("Vermont" = "VT", "Maryland" = "MD", "North_Carolina" = "NC")) +
theme(strip.text.x = element_text(size = 30))Region plot
# Rxn - population
budbreak_DOY_pop_rxn <- BB20_DOY + geom_line(aes(y=pred,group=Family,colour=Population)) +
xlab("Garden site") + ylab("DOY") +
theme_minimal() +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=12,face="bold"),
strip.text.x = element_text(size = 12, colour = "black", angle = 0, face="bold"))
## Reaction norm for geo genetic groups/regions
BB20_DOY_B <- ggplot(budbreak,
aes(x=Garden, y=BudBreak, color=Region)) +
stat_summary(aes(group = Region), fun.y = mean, geom = "path", size = 1.5) +
stat_summary(aes(color = Region), fun.data = mean_cl_boot, geom = "errorbar", width = 0.1, size=1.5) +
stat_summary(aes(color = Region), fun.y = mean, geom = "point", size = 5) +
# Lines by species using grouping
stat_summary(aes(group = Region), geom = "line", fun = mean) +
ylab("") + xlab("") +
theme(legend.position="right") +
theme_bw() + theme(axis.text=element_text(size=18),
axis.title=element_text(size=20,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
scale_color_manual(values = c("Core" = "goldenrod2","Margin"="green2", "Edge"="steelblue"))## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.y` is deprecated. Use `fun` instead.Combined plot - Supplementary figure
legend_PC <- get_legend(BB20_DOY_A)
legend_RE <- get_legend(BB20_DOY_B)
leg_full <- cowplot::plot_grid(legend_PC, legend_RE,
nrow = 2)
BB_DOY_grid <- cowplot::plot_grid(BB20_DOY_A + theme(legend.position="none"),
BB20_DOY_B + theme(legend.position="none"),
align = 'vh',
labels = c("A", "B"),
hjust = -1,
nrow = 1,
rel_widths = c(1,.7))## Warning: Graphs cannot be vertically aligned unless the axis parameter is set.
## Placing graphs unaligned.## Warning: Graphs cannot be horizontally aligned unless the axis parameter is set.
## Placing graphs unaligned.# bud set DOY final
cowplot::plot_grid(BB_DOY_grid, leg_full, ncol = 2, rel_widths = c(1,.08))