PCA and Plasticity
Anoob Prakash
02/06/2021
Methods
Adaptive plasticity model
- Height was selected as proxy for fitness
- Absolute plasticity or Plasticity of the trait correlated with fitness to check whether the plasticity is adaptive, maladaptive or neutral
Absolute_plasticity or Plasticity <- Height_Growth ~ Abs_P + Region + Abs_P * Region
Estimating the adaptive or maladaptive nature of trait plasticity
notes
gw - growth
bb - budbreak
bs - budset
B - lmer family blups (means)
P - plasticity value for the trait
packages
suppressPackageStartupMessages({
require(ggplot2)
require(ggpubr)
require(wesanderson)
require(factoextra)
require(dplyr)
require(lmerTest)
require(tidyverse)
require(viridis)
require(cowplot)
})# setwd("Z:/Anoob/MCMCglmm")
# data
budset_2019 <- read.csv("./trait_data/BudSet_2019.csv")
budset_2020 <- read.csv("./trait_data/BudSet_2020.csv")
budbreak_2020 <- read.csv("./trait_data/BudBreak_2020_cGDD.csv")
heightGrowth_2019 <- read.csv("./trait_data/Growth_2019.csv")
heightGrowth_2020 <- read.csv("./trait_data/Growth_2020.csv")
exome_meta <- read.table("./Exome/RS_Exome_metadata.txt", header = T, sep="\t")
Plasticity <- read.csv("./trait_data/Plasticity.csv", stringsAsFactors = T); Plasticity <- Plasticity[,-1]
Plasticity$Region <- factor(Plasticity$Region, levels=c("Core","Margin","Edge"))
# Region palette
reg_palette <- c("goldenrod2","green2","steelblue")# exome data
exome_meta <- read.table("./Exome/RS_Exome_metadata.txt", sep="\t", header = T)
# climateNA
Population_climateNA <- read.csv("./ClimateNA/Family/PCA_sprucePops.txt", header=T, sep="\t")
colnames(Population_climateNA)[colnames(Population_climateNA)=="PC1_sel"] <- "PC1"
Family_climateNA <- read.csv("./ClimateNA/Family/PCA_spruceFams.txt", header=T, sep="\t")
colnames(Family_climateNA)[colnames(Family_climateNA)=="PC1_sel"] <- "PC1"Climate data
Provenance climate
Provenance_MAT <- read.table("./ClimateNA/Family/PCA_spruceFams.txt", sep="\t", header = T, stringsAsFactors = T)
selVars<-c('DD_0','DD18','MAR','PAS','MSP','RH','EXT','CMD','TD','eFFP','current_30arcsec_annualPET') Garden climate
Garden_clim <- read.table("./ClimateNA/climNA_gardens_2019-20.txt", header=T)
Garden_clim## garden year MSP DD_0 DD5 DD18 PAS eFFP RH TD CMD MAR EXT PET
## 1 MD 2019 469 341 2649 368 94 297 66 25.4 212 15.2 36.4 866.3545
## 2 MD 2020 492 224 2548 362 44 299 67 23.7 131 14.8 36.4 872.7068
## 3 NC 2019 478 55 3873 803 5 324 64 20.4 277 17.1 37.7 1147.0480
## 4 NC 2020 723 52 3747 698 4 321 66 20.9 166 16.2 37.8 1094.4220
## 5 VT 2019 471 678 2451 407 172 293 67 31.9 147 13.9 39.6 746.0178
## 6 VT 2020 403 385 2633 483 74 294 67 28.8 226 14.2 38.9 793.8243selVars1 <- c('DD_0','DD18','MAR','PAS','MSP','RH','EXT','CMD','TD','eFFP','PET')PCA
# create ID for provenance (family level)
identity <- Provenance_MAT[,(1:2)]
colnames(exome_meta)[1] <- "Population"
fam_id <- right_join(identity,exome_meta)## Joining, by = c("Population", "Family")fam_id$Region <- factor(fam_id$Region,levels=c("C","M","E"))
fam_id$Region <- plyr::revalue(fam_id$Region, c("C"="Core","M"="Margin","E"="Edge"))
# selected clim of provenance
data_p2 <-Provenance_MAT[,selVars]
colnames(data_p2)[11] <- "PET"# create ID for garden means
# garden_identity <- Garden_climateNA1[,1]
# create ID for garden
garden_identity <- Garden_clim[,1:2]
# selected clim of gardens
data_p1 <- Garden_clim[,selVars1]# combined ID
head(fam_id)## Population Family Tree State Region Date Latitude Longitude Elevation
## 1 AB AB_05 5 TN Edge 9002016 35.55297 83.49438 1812
## 2 AB AB_08 8 TN Edge 9002016 35.55212 83.49259 1785
## 3 AB AB_12 12 TN Edge 9002016 35.53890 83.49463 1750
## 4 AB AB_16 16 TN Edge 9002016 35.53882 83.49534 1738
## 5 AB AB_18 18 TN Edge 9002016 35.53945 83.49375 1769
## 6 ALB ALB_01 1 ME Core 11112017 44.30670 70.84019 455
## DBH
## 1 57.7
## 2 44.3
## 3 36.8
## 4 44.0
## 5 68.2
## 6 28.1head(garden_identity)## garden year
## 1 MD 2019
## 2 MD 2020
## 3 NC 2019
## 4 NC 2020
## 5 VT 2019
## 6 VT 2020fam_id <- fam_id[,c("Population","Family","Region")]
fam_id$Garden <- NA
fam_id$Year <- NA
garden_id <- garden_identity
garden_id$Population <- NA
garden_id$Family <- NA
garden_id$Region <- NA
colnames(garden_id)[1] <- "Garden"
colnames(garden_id)[2] <- "Year"
garden_id <- garden_id[,c("Population","Family","Region","Garden","Year")]
id_c <- rbind(garden_id,fam_id)
# combined data
data_c <- rbind(data_p1,data_p2)
# combined PCA
data_pca <- prcomp(data_c,scale. = T)
# Extract PC axes for plotting
PCAvalues <- data.frame(Garden = id_c$Garden, Year = id_c$Year,
Population = id_c$Population, Family = id_c$Family, Region = id_c$Region,
data_pca$x)
PCAvalues$Year <- as.factor(PCAvalues$Year)
PCAvalues[!is.na(PCAvalues$Garden),"Group"] <- PCAvalues[!is.na(PCAvalues$Garden),"Garden"]
PCAvalues[!is.na(PCAvalues$Region),"Group"] <- PCAvalues[!is.na(PCAvalues$Region),"Region"]
PCAvalues$Group <- as.factor(PCAvalues$Group)
# Extract loadings of the variables
PCAloadings <- data.frame(Variables = rownames(data_pca$rotation), data_pca$rotation)
# write.csv(PCAvalues, row.names = F, "./trait_data/PCA/PCA_score.csv")
# write.csv(PCAloadings, row.names = F, "./trait_data/PCA/PCA_loadings.csv")# full plot
full_p <- ggplot(PCAvalues[!is.na(PCAvalues$Region),], aes(x = PC1, y = PC2)) +
# Provenance data
geom_point(aes(color=Group), size=4, pch = 20) +
# add loadings
geom_segment(data = PCAloadings, aes(x = 0, y = 0, xend = (PC1*7),
yend = (PC2*7)), arrow = arrow(length = unit(1/2, "picas")),
color = "black") +
annotate("text", x = (PCAloadings$PC1*7.5), y = (PCAloadings$PC2*7.5),
label = PCAloadings$Variables, size =5) +
scale_color_manual(values = c(
"Core" = "goldenrod2",
"Margin" = "green2",
"Edge" = "steelblue",
"VT" = "#F1BB7B",
"NC"="#5B1A18",
"MD"="#FD6467"),
breaks = c("Core","Margin","Edge","VT","MD","NC"),
labels = c("Core","Margin","Edge","Vermont", "Maryland", "North Carolina")) +
# Garden data
geom_point(data = PCAvalues[!is.na(PCAvalues$Garden),],
aes(x = PC1, y = PC2,
shape=Year,
color=Group), size=4) +
scale_shape_manual(values = c('2019' = 17, '2020' = 16)) +
labs(x = "PC1 (48%)", y = "PC2 (30.7%)") +
theme_bw() +
theme(legend.position='none',
axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) +
geom_vline(xintercept = 0, col="darkred", alpha=0.3)+
geom_hline(yintercept = 0, col="darkred", alpha=0.3) +
stat_ellipse(
mapping = aes(PC1, PC2, color = Region),
data = PCAvalues[!is.na(PCAvalues$Region),],
geom = "path",
position = "identity",
type = "t",
level = 0.95,
segments = 51,
na.rm = FALSE,
show.legend = NA,
inherit.aes = TRUE)
# create legend for garden
garden_legend_1 <- ggplot(data = PCAvalues[!is.na(PCAvalues$Garden),]) +
geom_point(aes(x = PC1, y = PC2,
shape=Year), size=4) +
scale_shape_manual(values = c('2019' = 17, '2020' = 16)) +
theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13))
garden_legend_2 <- ggplot(data = PCAvalues[!is.na(PCAvalues$Garden),]) +
geom_point(aes(x = PC1, y = PC2,
color=Group), size=4) +
scale_color_manual(values = c( "VT" = "#F1BB7B",
"NC"="#5B1A18",
"MD"="#FD6467"),
name = "Garden",
breaks = c("VT","MD","NC"),
labels = c("Vermont", "Maryland", "North Carolina")) +
theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13))
# create legend for Region
region_legend <- ggplot(data = PCAvalues[!is.na(PCAvalues$Region),]) +
geom_point(aes(x = PC1, y = PC2,
color=Group), size=4) +
scale_shape_manual(values = c('2019' = 17, '2020' = 16)) +
scale_color_manual(values = c(
"Core" = "goldenrod2",
"Margin" = "green2",
"Edge" = "steelblue"),
name = "Region",
breaks = c("Core","Margin","Edge"),
labels = c("Core","Margin","Edge")) +
theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13))
cowplot::plot_grid(
full_p
, cowplot::plot_grid(
get_legend(garden_legend_1),
get_legend(garden_legend_2)
, get_legend(region_legend)
, nrow = 1
)
, nrow = 2
, rel_heights = c(8,2)
)
Plasticity vs PC1
Plasticity vs PC1 Models
Budset 2019
# budset 2019
budset2019_PC1_mod1 <- lmer(BudSet_2019 ~ PC1 + (1|Population), data=Plasticity[!is.na(Plasticity$BudSet_2019),])
anova(budset2019_PC1_mod1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## PC1 0.03176 0.03176 1 61.243 0.0543 0.8165summary(budset2019_PC1_mod1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: BudSet_2019 ~ PC1 + (1 | Population)
## Data: Plasticity[!is.na(Plasticity$BudSet_2019), ]
##
## REML criterion at convergence: 807.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.50530 -0.69399 0.02162 0.64177 3.09669
##
## Random effects:
## Groups Name Variance Std.Dev.
## Population (Intercept) 0.07226 0.2688
## Residual 0.58483 0.7647
## Number of obs: 334, groups: Population, 65
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.16682 0.05403 58.98503 58.615 <2e-16 ***
## PC1 -0.00589 0.02528 61.24270 -0.233 0.817
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## PC1 -0.025Budset 2020
# budset 2020
budset2020_PC1_mod1 <- lmer(BudSet_2020 ~ PC1 + (1|Population), data=Plasticity[!is.na(Plasticity$BudSet_2020),])
anova(budset2020_PC1_mod1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## PC1 0.012373 0.012373 1 60.461 0.0968 0.7568summary(budset2020_PC1_mod1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: BudSet_2020 ~ PC1 + (1 | Population)
## Data: Plasticity[!is.na(Plasticity$BudSet_2020), ]
##
## REML criterion at convergence: 282.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9419 -0.6652 0.0111 0.7081 2.8398
##
## Random effects:
## Groups Name Variance Std.Dev.
## Population (Intercept) 0.004167 0.06455
## Residual 0.127860 0.35758
## Number of obs: 334, groups: Population, 65
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.832691 0.021250 57.985915 274.485 <2e-16 ***
## PC1 -0.003103 0.009974 60.461382 -0.311 0.757
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## PC1 -0.028Budbreak 2020
# budbreak 2020
budbreak2020_PC1_mod1 <- lmer(BudBreak_cGDD ~ PC1 + (1|Population), data=Plasticity[!is.na(Plasticity$BudBreak_cGDD),])
anova(budbreak2020_PC1_mod1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## PC1 560.27 560.27 1 66.249 33.366 2.226e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(budbreak2020_PC1_mod1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: BudBreak_cGDD ~ PC1 + (1 | Population)
## Data: Plasticity[!is.na(Plasticity$BudBreak_cGDD), ]
##
## REML criterion at convergence: 1954.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6209 -0.5602 0.0713 0.6796 1.8507
##
## Random effects:
## Groups Name Variance Std.Dev.
## Population (Intercept) 5.757 2.399
## Residual 16.792 4.098
## Number of obs: 334, groups: Population, 65
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 28.0225 0.3764 62.9877 74.456 < 2e-16 ***
## PC1 -1.0115 0.1751 66.2493 -5.776 2.23e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## PC1 -0.021Growth 2019
# growth 2019
growth2019_PC1_mod1 <- lmer(Growth_2019 ~ PC1 + (1|Population), data=Plasticity[!is.na(Plasticity$Growth_2019),])
anova(growth2019_PC1_mod1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## PC1 0.44243 0.44243 1 68.98 60.936 4.529e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(growth2019_PC1_mod1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Growth_2019 ~ PC1 + (1 | Population)
## Data: Plasticity[!is.na(Plasticity$Growth_2019), ]
##
## REML criterion at convergence: -628.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8264 -0.6755 -0.1184 0.6271 3.4848
##
## Random effects:
## Groups Name Variance Std.Dev.
## Population (Intercept) 0.001847 0.04298
## Residual 0.007261 0.08521
## Number of obs: 334, groups: Population, 65
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.164158 0.007156 66.049367 22.939 < 2e-16 ***
## PC1 -0.026047 0.003337 68.980217 -7.806 4.53e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## PC1 -0.022Growth 2020
# growth 2020
growth2020_PC1_mod1 <- lmer(Growth_2020 ~ PC1 + (1|Population), data=Plasticity[!is.na(Plasticity$Growth_2020),])
anova(growth2020_PC1_mod1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## PC1 0.16999 0.16999 1 67.73 3.5533 0.06372 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(growth2020_PC1_mod1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Growth_2020 ~ PC1 + (1 | Population)
## Data: Plasticity[!is.na(Plasticity$Growth_2020), ]
##
## REML criterion at convergence: 19
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.89344 -0.47660 0.07766 0.61271 2.32820
##
## Random effects:
## Groups Name Variance Std.Dev.
## Population (Intercept) 0.02092 0.1446
## Residual 0.04784 0.2187
## Number of obs: 334, groups: Population, 65
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.84749 0.02177 63.92545 38.934 <2e-16 ***
## PC1 0.01905 0.01011 67.73023 1.885 0.0637 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## PC1 -0.019Plasticity vs PC1 Viz
Phenology
pheno_plastic <- subset(Plasticity, select=-c(BudBreak_DOY,Growth_2019,Growth_2020)) # growth_2019,growth_2020
colnames(pheno_plastic)[colnames(pheno_plastic)=="BudBreak_cGDD"] <- "BudBreak_2020"
Plasticity_pheno <- pheno_plastic %>% tidyr::pivot_longer(cols=c("BudBreak_2020":"BudSet_2020"),
names_to = "Traits",
values_to = "value")
Plasticity_pheno$TraitType <- ifelse(startsWith(Plasticity_pheno$Traits, "BudSet"), "BudSet", "BudBreak")
Plasticity_pheno$Year <- ifelse(endsWith(Plasticity_pheno$Traits, "2020"), "2020", "2019")
# Plasticity_pheno <- dplyr::rename(Plasticity_pheno, PC1 = PC1_sel)
# Population means
Provenance_pop <- read.table("./ClimateNA/Family/PCA_sprucePops.txt", sep="\t", header = T, stringsAsFactors = T)
Provenance_pop <- subset(Provenance_pop, select=c(Population,PC1_sel,MAT))
Provenance_pop <- dplyr::rename(Provenance_pop,c(PC1_pop=PC1_sel,MAT_pop=MAT))
# Add population means to the main dataframe
Plasticity_pheno <- left_join(Plasticity_pheno,Provenance_pop)## Joining, by = "Population"# combined
Phenology <- ggplot(data=Plasticity_pheno,
mapping = aes(x=PC1,y=abs(value),group=Traits)) +
geom_smooth(method='lm', formula= y~x, show.legend = FALSE, col="darkred") + theme_minimal() +
geom_point(aes(colour = Region, shape=TraitType), size = 4, alpha = .4) +
scale_color_manual(values = c("Core" = "goldenrod2","Margin"="green2","Edge"="steelblue")) +
labs(x = "PC1", y = "Plasticity (CV)") +
theme_bw(base_size = 11, base_family = "Times") +
theme(axis.text=element_text(size=14),
axis.title=element_text(size=18),
panel.background = element_blank(),
legend.background = element_blank(),
panel.grid = element_blank(),
plot.background = element_blank(),
legend.text=element_text(size=rel(.8)),
strip.text = element_text(size=30)) +
facet_grid(.~Year) +
theme(strip.text.x = element_text(size = 30))
# create legend
test1 <- ggplot(data=Plasticity_pheno,
mapping = aes(x=PC1,y=abs(value),group=Traits)) +
geom_smooth(method='lm', formula= y~x, show.legend = FALSE, col="darkred") + theme_minimal() +
geom_point(aes(colour = Region, shape=TraitType), size = 3, alpha = 1) +
scale_color_manual(values = c("Core" = "goldenrod2","Margin"="green2","Edge"="steelblue")) +
labs(x = "PC1", y = "Absolute Plasticity") +
scale_color_manual(values = c("Core" = "goldenrod2","Margin"="green2","Edge"="steelblue")) +
theme_bw() + theme(axis.text=element_text(size=10),
axis.title=element_text(size=12),
legend.position = "right",
legend.title=element_text(size=21),
legend.text=element_text(size=20),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
plot.title = element_blank()) +
guides(color = guide_legend(override.aes = list(size = 8)))## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.legend_ba <- get_legend(test1 + theme(legend.position="right"))
# final
cowplot::plot_grid(Phenology+ theme(legend.position="none"), legend_ba, ncol = 2, rel_widths = c(1,.15))## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
Growth
growth_plastic <- subset(Plasticity, select=-c(BudBreak_DOY,BudBreak_cGDD,BudSet_2020,BudSet_2019));
growth_plastic <- growth_plastic %>% tidyr::pivot_longer(cols=c("Growth_2019":"Growth_2020"),
names_to = "Traits",
values_to = "value")
growth_plastic$TraitType <- "Growth"
growth_plastic$Year <- ifelse(endsWith(growth_plastic$Traits, "2020"), "2020", "2019")
# growth_plastic <- dplyr::rename(growth_plastic, PC1 = PC1_sel)
# Add population means to the main dataframe
growth_plastic <- left_join(growth_plastic,Provenance_pop)## Joining, by = "Population"# growth
Growth <- ggplot(data=growth_plastic,
mapping = aes(x=PC1,y=abs(value),group=Traits)) +
geom_smooth(method='lm', formula= y~x, show.legend = FALSE, col="darkred") + theme_minimal() +
geom_point(aes(colour = Region), size = 4, alpha = .4) +
scale_color_manual(values = c("Core" = "goldenrod2","Margin"="green2","Edge"="steelblue",
"growth2019_Plasticity"="darkred", "growth2020_Plasticity"="limegreen")) +
labs(x = "PC1", y = "Plasticity (CV)") +
theme_bw(base_size = 11, base_family = "Times") +
theme(axis.text=element_text(size=14),
axis.title=element_text(size=18),
panel.background = element_blank(),
legend.background = element_blank(),
panel.grid = element_blank(),
plot.background = element_blank(),
legend.text=element_text(size=rel(.8)),
strip.text = element_text(size=30)) +
facet_grid(.~Year)
# create legend
test2 <- ggplot(data=growth_plastic,
mapping = aes(x=PC1,y=abs(value),group=Traits)) +
geom_smooth(method='lm', formula= y~x, show.legend = FALSE, col="darkred") + theme_minimal() +
geom_point(aes(colour = Region), size = 3, alpha = 1) +
scale_color_manual(values = c("Core" = "goldenrod2","Margin"="green2","Edge"="steelblue")) +
theme_bw() + theme(axis.text=element_text(size=10),
axis.title=element_text(size=12),
legend.position = "right",
legend.title=element_text(size=21),
legend.text=element_text(size=20),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
plot.title = element_blank()) +
guides(color = guide_legend(override.aes = list(size = 8)))
legend_bb <- get_legend(test2 + theme(legend.position="right"))
# final
cowplot::plot_grid(Growth+ theme(legend.position="none"), legend_bb, ncol = 2, rel_widths = c(1,.09))## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
ANOVA
sjPlot::tab_model(budset2019_PC1_mod1, budset2020_PC1_mod1, budbreak2020_PC1_mod1,growth2019_PC1_mod1, growth2020_PC1_mod1,
dv.labels = c("Budset 2019", "Budset 2020", "Budbreak 2020", "Growth 2019", "Growth 2020"),
string.pred = "Coeffcient",
string.ci = "CI (95%)",
string.p = "P-Value",
p.style = "stars")| Â | Budset 2019 | Budset 2020 | Budbreak 2020 | Growth 2019 | Growth 2020 | |||||
|---|---|---|---|---|---|---|---|---|---|---|
| Coeffcient | Estimates | CI (95%) | Estimates | CI (95%) | Estimates | CI (95%) | Estimates | CI (95%) | Estimates | CI (95%) |
| (Intercept) | 3.17 *** | 3.06 – 3.27 | 5.83 *** | 5.79 – 5.87 | 28.02 *** | 27.28 – 28.76 | 0.16 *** | 0.15 – 0.18 | 0.85 *** | 0.80 – 0.89 |
| PC1 | -0.01 | -0.06 – 0.04 | -0.00 | -0.02 – 0.02 | -1.01 *** | -1.35 – -0.67 | -0.03 *** | -0.03 – -0.02 | 0.02 | -0.00 – 0.04 |
| Random Effects | ||||||||||
| σ2 | 0.58 | 0.13 | 16.79 | 0.01 | 0.05 | |||||
| Ï„00 | 0.07 Population | 0.00 Population | 5.76 Population | 0.00 Population | 0.02 Population | |||||
| ICC | 0.11 | 0.03 | 0.26 | 0.20 | 0.30 | |||||
| N | 65 Population | 65 Population | 65 Population | 65 Population | 65 Population | |||||
| Observations | 334 | 334 | 334 | 334 | 334 | |||||
| Marginal R2 / Conditional R2 | 0.000 / 0.110 | 0.000 / 0.032 | 0.171 / 0.382 | 0.252 / 0.404 | 0.023 / 0.320 | |||||
|
||||||||||
Figure for paper
## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family not
## found in Windows font database
## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family not
## found in Windows font database## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family not
## found in Windows font database
## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family not
## found in Windows font database
## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family not
## found in Windows font database## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
Plasticity vs Height Growth
Height growth of 2019
- Means vs Plasticity
Family data
# blups for height growth
heightGrowth_2019_2 <- heightGrowth_2019
heightGrowth_2019_2 <- heightGrowth_2019_2%>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)
HeightGrowth_2019_mod <- lmer(Growth ~ Garden + (1|mBed) + (1|Family), heightGrowth_2019_2)
HeightGrowth_2019_blup <- ranef(HeightGrowth_2019_mod)
HeightGrowth_2019_blup <- HeightGrowth_2019_blup$Family
HeightGrowth_2019_blup <- cbind(Family=rownames(HeightGrowth_2019_blup),HeightGrowth_2019_blup)
rownames(HeightGrowth_2019_blup) <- NULL
names(HeightGrowth_2019_blup)[2] <- "Height_Growth"
# plasticity based on height growth
HeightGrowth_2019_blup <- merge(x=HeightGrowth_2019_blup,
y=Plasticity[,c("Family","Growth_2019")],
by="Family")
# rename plasticity column
colnames(HeightGrowth_2019_blup)[colnames(HeightGrowth_2019_blup)=="Growth_2019"] <- "Plasticity"
# remove NAs
HeightGrowth_2019_blup <- HeightGrowth_2019_blup[!is.na(HeightGrowth_2019_blup$Plasticity),]
# Population column creation
HeightGrowth_2019_blup$Population <- gsub("\\_.*", "", HeightGrowth_2019_blup$Family)
# merge climateNA data
HeightGrowth_2019_blup <- merge(x=HeightGrowth_2019_blup,
y=Family_climateNA[,c("Population","PC1","MAT")],
by="Population")
# merge Region data
HeightGrowth_2019_blup <- merge(x=HeightGrowth_2019_blup,
y=heightGrowth_2019[,c("Population","Region")],
by="Population")
# reorder regions
HeightGrowth_2019_blup$Region <- factor(HeightGrowth_2019_blup$Region,levels = c("Core", "Margin", "Edge"))
HeightGrowth_2019_blup <- distinct(HeightGrowth_2019_blup)
HeightGrowth_2019_blup$Abs_P <- abs(HeightGrowth_2019_blup$Plasticity)Population data
HeightGrowth_2019_blup$Population <- as.factor(HeightGrowth_2019_blup$Population)
HeightGrowth_2019_pop <- HeightGrowth_2019_blup %>%
select(Population, Height_Growth, Plasticity, Region) %>%
group_by(Population) %>%
dplyr::summarise(Height_Growth = mean(Height_Growth), Plasticity = mean(Plasticity))
HeightGrowth_2019_pop <- merge(x=HeightGrowth_2019_pop,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
HeightGrowth_2019_pop <- merge(x=HeightGrowth_2019_pop,
y=heightGrowth_2019[,c("Population","Region")],
by="Population")
# reorder regions
HeightGrowth_2019_pop$Region <- factor(HeightGrowth_2019_pop$Region,levels = c("Core", "Margin", "Edge"))
HeightGrowth_2019_pop <- distinct(HeightGrowth_2019_pop)
HeightGrowth_2019_pop$Abs_P <- abs(HeightGrowth_2019_pop$Plasticity)Family plot
HeightGrowth_2019_fam_abs_mod <- lm(Height_Growth ~ Abs_P + Region + Abs_P*Region,
data=HeightGrowth_2019_blup)
summary(HeightGrowth_2019_fam_abs_mod)##
## Call:
## lm(formula = Height_Growth ~ Abs_P + Region + Abs_P * Region,
## data = HeightGrowth_2019_blup)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5173 -0.0751 0.0890 0.2742 0.7124
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.44017 0.03362 -102.320 < 2e-16 ***
## Abs_P 20.84251 0.25121 82.970 < 2e-16 ***
## RegionMargin 0.26170 0.09085 2.881 0.00402 **
## RegionEdge 0.40777 0.06659 6.124 1.11e-09 ***
## Abs_P:RegionMargin -1.17380 0.37479 -3.132 0.00176 **
## Abs_P:RegionEdge -1.52333 0.35440 -4.298 1.81e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6059 on 1851 degrees of freedom
## Multiple R-squared: 0.9367, Adjusted R-squared: 0.9365
## F-statistic: 5474 on 5 and 1851 DF, p-value: < 2.2e-16- Absolute plasticity is not significant.
- Absolute plasticity’s interaction with Core and Edge region are significant***.
Population plot
HeightGrowth_2019_pop_abs_mod <- lm(Height_Growth ~ Abs_P + Region + Abs_P*Region,
data=HeightGrowth_2019_pop)
summary(HeightGrowth_2019_pop_abs_mod)##
## Call:
## lm(formula = Height_Growth ~ Abs_P + Region + Abs_P * Region,
## data = HeightGrowth_2019_pop)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.79780 -0.05818 0.00806 0.12916 0.61914
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.0053 0.1501 -26.691 < 2e-16 ***
## Abs_P 25.9086 1.3238 19.571 < 2e-16 ***
## RegionMargin 0.8896 0.3405 2.613 0.011382 *
## RegionEdge 1.0275 0.3312 3.103 0.002942 **
## Abs_P:RegionMargin -6.3872 1.7397 -3.671 0.000521 ***
## Abs_P:RegionEdge -6.8532 1.9572 -3.502 0.000888 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2934 on 59 degrees of freedom
## Multiple R-squared: 0.9721, Adjusted R-squared: 0.9697
## F-statistic: 410.8 on 5 and 59 DF, p-value: < 2.2e-16- Absolute plasticity is not significant.
- Absolute plasticity’s interaction with Core, Margin and Edge region are significant*.
Height growth of 2020
- Means vs Plasticity
Family data
# blups for height growth
heightGrowth_2020_2 <- heightGrowth_2020
heightGrowth_2020_2 <- heightGrowth_2020_2%>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)
HeightGrowth_2020_mod <- lmer(Growth ~ Garden + (1|mBed) + (1|Family), heightGrowth_2020_2)
HeightGrowth_2020_blup <- ranef(HeightGrowth_2020_mod)
HeightGrowth_2020_blup <- HeightGrowth_2020_blup$Family
HeightGrowth_2020_blup <- cbind(Family=rownames(HeightGrowth_2020_blup),HeightGrowth_2020_blup)
rownames(HeightGrowth_2020_blup) <- NULL
names(HeightGrowth_2020_blup)[2] <- "Height_Growth"
# plasticity based on height growth
HeightGrowth_2020_blup <- merge(x=HeightGrowth_2020_blup,
y=Plasticity[,c("Family","Growth_2020")],
by="Family")
# rename plasticity column
colnames(HeightGrowth_2020_blup)[colnames(HeightGrowth_2020_blup)=="Growth_2020"] <- "Plasticity"
# remove NAs
HeightGrowth_2020_blup <- HeightGrowth_2020_blup[!is.na(HeightGrowth_2020_blup$Plasticity),]
# Population column creation
HeightGrowth_2020_blup$Population <- gsub("\\_.*", "", HeightGrowth_2020_blup$Family)
# merge climateNA data
HeightGrowth_2020_blup <- merge(x=HeightGrowth_2020_blup,
y=Family_climateNA[,c("Population","PC1","MAT")],
by="Population")
# merge Region data
HeightGrowth_2020_blup <- merge(x=HeightGrowth_2020_blup,
y=heightGrowth_2020[,c("Population","Region")],
by="Population")
# reorder regions
HeightGrowth_2020_blup$Region <- factor(HeightGrowth_2020_blup$Region,levels = c("Core", "Margin", "Edge"))
HeightGrowth_2020_blup <- distinct(HeightGrowth_2020_blup)
HeightGrowth_2020_blup$Abs_P <- abs(HeightGrowth_2020_blup$Plasticity)Population data
HeightGrowth_2020_blup$Population <- as.factor(HeightGrowth_2020_blup$Population)
HeightGrowth_2020_pop <- HeightGrowth_2020_blup %>%
select(Population, Height_Growth, Plasticity, Region) %>%
group_by(Population) %>%
dplyr::summarise(Height_Growth = mean(Height_Growth), Plasticity = mean(Plasticity))
HeightGrowth_2020_pop <- merge(x=HeightGrowth_2020_pop,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
HeightGrowth_2020_pop <- merge(x=HeightGrowth_2020_pop,
y=heightGrowth_2019[,c("Population","Region")],
by="Population")
# reorder regions
HeightGrowth_2020_pop$Region <- factor(HeightGrowth_2020_pop$Region,levels = c("Core", "Margin", "Edge"))
HeightGrowth_2020_pop <- distinct(HeightGrowth_2020_pop)
HeightGrowth_2020_pop$Abs_P <- abs(HeightGrowth_2020_pop$Plasticity)Family plot
HeightGrowth_2020_fam_abs_mod <- lm(Height_Growth ~ Abs_P + Region + Abs_P*Region,
data=HeightGrowth_2020_blup)
summary(HeightGrowth_2020_fam_abs_mod)##
## Call:
## lm(formula = Height_Growth ~ Abs_P + Region + Abs_P * Region,
## data = HeightGrowth_2020_blup)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.6286 -0.4908 -0.0495 0.3532 15.5621
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.40454 0.22541 28.413 < 2e-16 ***
## Abs_P -7.89036 0.23892 -33.024 < 2e-16 ***
## RegionMargin -1.06800 0.31047 -3.440 0.000595 ***
## RegionEdge -0.06272 0.41411 -0.151 0.879626
## Abs_P:RegionMargin 5.05131 0.41562 12.154 < 2e-16 ***
## Abs_P:RegionEdge 0.03699 0.45338 0.082 0.934984
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.791 on 1851 degrees of freedom
## Multiple R-squared: 0.6359, Adjusted R-squared: 0.635
## F-statistic: 646.7 on 5 and 1851 DF, p-value: < 2.2e-16- Absolute plasticity is significant***.
- Absolute plasticity’s interaction with Core and Margin region is significant***.
Population plot
HeightGrowth_2020_pop_abs_mod <- lm(Height_Growth ~ Abs_P + Region + Abs_P*Region,
data=HeightGrowth_2020_pop)
summary(HeightGrowth_2020_pop_abs_mod)##
## Call:
## lm(formula = Height_Growth ~ Abs_P + Region + Abs_P * Region,
## data = HeightGrowth_2020_pop)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.8747 -0.1743 -0.0479 0.1071 9.1229
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.4983 1.6483 3.942 0.000217 ***
## Abs_P -7.9608 1.7941 -4.437 4.04e-05 ***
## RegionMargin -3.3675 2.4740 -1.361 0.178649
## RegionEdge -0.3484 2.8208 -0.124 0.902120
## Abs_P:RegionMargin 8.6277 3.6101 2.390 0.020066 *
## Abs_P:RegionEdge 0.3684 3.1435 0.117 0.907113
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.393 on 59 degrees of freedom
## Multiple R-squared: 0.6384, Adjusted R-squared: 0.6078
## F-statistic: 20.83 on 5 and 59 DF, p-value: 6.245e-12- Absolute plasticity is significant*.
- Absolute plasticity’s interaction with Margin region is significant***.
Budbreak 2020 Plasticity
Family data
# blups for height growth
heightGrowth_2020_2 <- heightGrowth_2020
heightGrowth_2020_2 <- heightGrowth_2020_2 %>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)
HeightGrowth_2020_mod <- lmer(Growth ~ Garden + (1|mBed) + (1|Family), heightGrowth_2020_2)
HeightGrowth_2020_blup <- ranef(HeightGrowth_2020_mod)
HeightGrowth_2020_blup <- HeightGrowth_2020_blup$Family
HeightGrowth_2020_blup <- cbind(Family=rownames(HeightGrowth_2020_blup),HeightGrowth_2020_blup)
rownames(HeightGrowth_2020_blup) <- NULL
names(HeightGrowth_2020_blup)[2] <- "Height_Growth"
# plasticity based on budbreak 2020
gw_20_B_bb_20_P <- merge(x=HeightGrowth_2020_blup,
y=Plasticity[,c("Family","BudBreak_cGDD")],
by="Family")
# rename plasticity column if needed
colnames(gw_20_B_bb_20_P)[colnames(gw_20_B_bb_20_P) == "BudBreak_cGDD"] <- "BudBreak2020_Plasticity"
# remove NAs
gw_20_B_bb_20_P <- gw_20_B_bb_20_P[!is.na(gw_20_B_bb_20_P$BudBreak2020_Plasticity),]
# Population column creation
gw_20_B_bb_20_P$Population <- gsub("\\_.*", "", gw_20_B_bb_20_P$Family)
# merge climateNA data
gw_20_B_bb_20_P <- merge(x=gw_20_B_bb_20_P,
y=Family_climateNA[,c("Population","PC1","MAT")],
by="Population")
# merge Region data
gw_20_B_bb_20_P <- merge(x=gw_20_B_bb_20_P,
y=heightGrowth_2020[,c("Population","Region")],
by="Population")
# reorder regions
gw_20_B_bb_20_P$Region <- factor(gw_20_B_bb_20_P$Region,levels = c("Core", "Margin", "Edge"))
gw_20_B_bb_20_P <- distinct(gw_20_B_bb_20_P)
gw_20_B_bb_20_P$Abs_P <- abs(gw_20_B_bb_20_P$BudBreak2020_Plasticity)Population data
gw_20_B_bb_20_P$Population <- as.factor(gw_20_B_bb_20_P$Population)
gw_20_B_bb_20_P_pop <- gw_20_B_bb_20_P %>%
select(Population, Height_Growth, BudBreak2020_Plasticity, Region) %>%
group_by(Population) %>%
dplyr::summarise(Height_Growth = mean(Height_Growth), BudBreak2020_Plasticity = mean(BudBreak2020_Plasticity))
gw_20_B_bb_20_P_pop <- merge(x=gw_20_B_bb_20_P_pop,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
gw_20_B_bb_20_P_pop <- merge(x=gw_20_B_bb_20_P_pop,
y=heightGrowth_2019[,c("Population","Region")],
by="Population")
# reorder regions
gw_20_B_bb_20_P_pop$Region <- factor(gw_20_B_bb_20_P_pop$Region,levels = c("Core", "Margin", "Edge"))
gw_20_B_bb_20_P_pop <- distinct(gw_20_B_bb_20_P_pop)
gw_20_B_bb_20_P_pop$Abs_P <- abs(gw_20_B_bb_20_P_pop$BudBreak2020_Plasticity)Family plot
gw_20_B_bb_20_P_U <- dplyr::distinct(gw_20_B_bb_20_P, Abs_P,Height_Growth, .keep_all=T)
gw_20_B_bb_20_P_fam_abs_mod <- lm(Height_Growth ~ Abs_P + Region + Abs_P*Region,
data=gw_20_B_bb_20_P_U)
summary(gw_20_B_bb_20_P_fam_abs_mod)##
## Call:
## lm(formula = Height_Growth ~ Abs_P + Region + Abs_P * Region,
## data = gw_20_B_bb_20_P_U)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.7710 -1.3773 -0.2445 0.9990 13.7460
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.52943 1.00317 -2.521 0.01216 *
## Abs_P 0.06750 0.03807 1.773 0.07710 .
## RegionMargin 10.11687 1.87943 5.383 1.4e-07 ***
## RegionEdge -3.27404 2.31683 -1.413 0.15856
## Abs_P:RegionMargin -0.21133 0.06983 -3.026 0.00267 **
## Abs_P:RegionEdge 0.09805 0.07583 1.293 0.19692
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.336 on 328 degrees of freedom
## Multiple R-squared: 0.3493, Adjusted R-squared: 0.3394
## F-statistic: 35.22 on 5 and 328 DF, p-value: < 2.2e-16anova(gw_20_B_bb_20_P_fam_abs_mod)## Analysis of Variance Table
##
## Response: Height_Growth
## Df Sum Sq Mean Sq F value Pr(>F)
## Abs_P 1 0.01 0.01 0.0010 0.974489
## Region 2 884.43 442.22 81.0307 < 2.2e-16 ***
## Abs_P:Region 2 76.52 38.26 7.0109 0.001043 **
## Residuals 328 1790.02 5.46
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Absolute plasticity is not significant.
- Absolute plasticity’s interaction with Margin region is significant***.
bb_A1 <- sjPlot::plot_model(gw_20_B_bb_20_P_fam_abs_mod, type = "int", colors=reg_palette) +
labs(x = "Plasticity of bud break in 2020", y = "Height Growth (BLUP)") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=18),
legend.position = "bottom",
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
plot.title = element_blank())
bb_A <- bb_A1 + geom_point(data=gw_20_B_bb_20_P_U,
aes(x = abs(BudBreak2020_Plasticity), y = Height_Growth, color = Region,
shape = Region), size=4, alpha = 0.25,
inherit.aes = FALSE)
bb_A
Population plot
gw_20_B_bb_20_P_pop_U <- dplyr::distinct(gw_20_B_bb_20_P_pop, Abs_P,Height_Growth, .keep_all=T)
gw_20_B_bb_20_P_pop_abs_mod <- lm(Height_Growth ~ Abs_P + Region + Abs_P*Region,
data=gw_20_B_bb_20_P_pop_U)
summary(gw_20_B_bb_20_P_pop_abs_mod)##
## Call:
## lm(formula = Height_Growth ~ Abs_P + Region + Abs_P * Region,
## data = gw_20_B_bb_20_P_pop_U)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5590 -0.8155 -0.0406 0.5900 6.0057
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.42312 2.62312 -0.543 0.589498
## Abs_P 0.02664 0.10087 0.264 0.792659
## RegionMargin 20.87003 4.26543 4.893 8.04e-06 ***
## RegionEdge -7.13249 5.39241 -1.323 0.191042
## Abs_P:RegionMargin -0.62286 0.16042 -3.883 0.000263 ***
## Abs_P:RegionEdge 0.22737 0.17977 1.265 0.210921
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.414 on 59 degrees of freedom
## Multiple R-squared: 0.6275, Adjusted R-squared: 0.5959
## F-statistic: 19.88 on 5 and 59 DF, p-value: 1.462e-11anova(gw_20_B_bb_20_P_pop_abs_mod)## Analysis of Variance Table
##
## Response: Height_Growth
## Df Sum Sq Mean Sq F value Pr(>F)
## Abs_P 1 8.740 8.740 4.3731 0.04083 *
## Region 2 143.922 71.961 36.0075 6.014e-11 ***
## Abs_P:Region 2 45.975 22.987 11.5023 6.052e-05 ***
## Residuals 59 117.911 1.998
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Absolute plasticity is not significant.
- Absolute plasticity’s interaction with Margin region is significant***.
# updated fig for the paper
bb_B1 <- sjPlot::plot_model(gw_20_B_bb_20_P_pop_abs_mod, type = "int", colors=reg_palette) +
labs(x = "Absolute plasticity of bud break in 2020 (Population)", y = "") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=18),
legend.position = "bottom",
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
plot.title = element_blank())
bb_B <- bb_B1 + geom_point(data=gw_20_B_bb_20_P_pop_U,
aes(x = abs(BudBreak2020_Plasticity), y = Height_Growth, color = Region,
shape = Region), size=4, alpha = 0.25,
inherit.aes = FALSE)
bb_B
# combined for the paper
# create legend
test <- ggplot() + geom_point(data=gw_20_B_bb_20_P_U,
aes(x = abs(BudBreak2020_Plasticity), y = Height_Growth, color = Region,
shape = Region), size=3, alpha = 1,
inherit.aes = FALSE) +
scale_color_manual(values = c("goldenrod2", "green2", "steelblue")) +
theme_bw() + theme(axis.text=element_text(size=10),
axis.title=element_text(size=12),
legend.position = "bottom",
legend.title=element_text(size=21),
legend.text=element_text(size=20),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
plot.title = element_blank()) +
guides(color = guide_legend(override.aes = list(size = 8)))
legend_bb <- get_legend(test + theme(legend.position="bottom"))
BB_grid <- cowplot::plot_grid( bb_A + theme(legend.position="none"),
bb_B + theme(legend.position="none"),
align = 'vh',
labels = c("A", "B"),
hjust = -1,
nrow = 1
)
# bud break final
bb_comb <- cowplot::plot_grid(BB_grid, legend_bb, ncol = 1, rel_heights = c(1,.2))Budset 2019 Plasticity
Family data
# blups for height growth
heightGrowth_2019_2 <- heightGrowth_2019
heightGrowth_2019_2 <- heightGrowth_2019_2%>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)
HeightGrowth_2019_mod <- lmer(Growth ~ Garden + (1|mBed) + (1|Family), heightGrowth_2019_2)
# HeightGrowth_2019_blup <- coef(HeightGrowth_2019_mod)
HeightGrowth_2019_blup <- ranef(HeightGrowth_2019_mod)
HeightGrowth_2019_blup <- HeightGrowth_2019_blup$Family
HeightGrowth_2019_blup <- cbind(Family=rownames(HeightGrowth_2019_blup),HeightGrowth_2019_blup)
rownames(HeightGrowth_2019_blup) <- NULL
names(HeightGrowth_2019_blup)[2] <- "Height_Growth"
# plasticity based on budset 2019
gw_19_B_bs_19_P <- merge(x=HeightGrowth_2019_blup,
y=Plasticity[,c("Family","BudSet_2019")],
by="Family")
# rename plasticity column if needed
colnames(gw_19_B_bs_19_P)[colnames(gw_19_B_bs_19_P)=="BudSet_2019"] <- "BudSet2019_Plasticity"
# remove NAs
gw_19_B_bs_19_P <- gw_19_B_bs_19_P[!is.na(gw_19_B_bs_19_P$BudSet2019_Plasticity),]
# Population column creation
gw_19_B_bs_19_P$Population <- gsub("\\_.*", "", gw_19_B_bs_19_P$Family)
# merge climateNA data
gw_19_B_bs_19_P <- merge(x=gw_19_B_bs_19_P,
y=Family_climateNA[,c("Population","PC1","MAT")],
by="Population")
# merge Region data
gw_19_B_bs_19_P <- merge(x=gw_19_B_bs_19_P,
y=heightGrowth_2019[,c("Population","Region")],
by="Population")
# reorder regions
gw_19_B_bs_19_P$Region <- factor(gw_19_B_bs_19_P$Region,levels = c("Core", "Margin", "Edge"))
gw_19_B_bs_19_P <- distinct(gw_19_B_bs_19_P)
gw_19_B_bs_19_P$Abs_P <- abs(gw_19_B_bs_19_P$BudSet2019_Plasticity)Population data
gw_19_B_bs_19_P$Population <- as.factor(gw_19_B_bs_19_P$Population)
gw_19_B_bs_19_P_pop <- gw_19_B_bs_19_P %>%
select(Population, Height_Growth, BudSet2019_Plasticity, Region) %>%
group_by(Population) %>%
dplyr::summarise(Height_Growth = mean(Height_Growth), BudSet2019_Plasticity = mean(BudSet2019_Plasticity))
gw_19_B_bs_19_P_pop <- merge(x=gw_19_B_bs_19_P_pop,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
gw_19_B_bs_19_P_pop <- merge(x=gw_19_B_bs_19_P_pop,
y=heightGrowth_2019[,c("Population","Region")],
by="Population")
# reorder regions
gw_19_B_bs_19_P_pop$Region <- factor(gw_19_B_bs_19_P_pop$Region,levels = c("Core", "Margin", "Edge"))
gw_19_B_bs_19_P_pop <- distinct(gw_19_B_bs_19_P_pop)
gw_19_B_bs_19_P_pop$Abs_P <- abs(gw_19_B_bs_19_P_pop$BudSet2019_Plasticity)Family plot
gw_19_B_bs_19_P_U <- dplyr::distinct(gw_19_B_bs_19_P, Abs_P,Height_Growth, .keep_all=T)
gw_19_B_bs_19_P_fam_abs_mod <- lm(Height_Growth ~ Abs_P + Region + Abs_P*Region,
data=gw_19_B_bs_19_P_U)
summary(gw_19_B_bs_19_P_fam_abs_mod)##
## Call:
## lm(formula = Height_Growth ~ Abs_P + Region + Abs_P * Region,
## data = gw_19_B_bs_19_P_U)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.8943 -1.3822 -0.0872 1.2888 6.4847
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.9894 0.5542 -3.590 0.000381 ***
## Abs_P 0.2514 0.1709 1.471 0.142265
## RegionMargin 3.8057 1.2747 2.986 0.003044 **
## RegionEdge 4.5872 1.0435 4.396 1.49e-05 ***
## Abs_P:RegionMargin -0.1640 0.3677 -0.446 0.655832
## Abs_P:RegionEdge -0.7991 0.3247 -2.461 0.014363 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.938 on 328 degrees of freedom
## Multiple R-squared: 0.3241, Adjusted R-squared: 0.3138
## F-statistic: 31.45 on 5 and 328 DF, p-value: < 2.2e-16anova(gw_19_B_bs_19_P_fam_abs_mod)## Analysis of Variance Table
##
## Response: Height_Growth
## Df Sum Sq Mean Sq F value Pr(>F)
## Abs_P 1 8.31 8.313 2.2125 0.13786
## Region 2 559.67 279.836 74.4775 < 2e-16 ***
## Abs_P:Region 2 22.86 11.428 3.0416 0.04911 *
## Residuals 328 1232.40 3.757
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Absolute plasticity is significant***.
- Absolute plasticity’s interaction with Core and Edge region is significant***.
# updated fig for the paper
bs_A1 <- sjPlot::plot_model(gw_19_B_bs_19_P_fam_abs_mod, type = "int", colors=reg_palette) +
labs(x = "Plasticity of bud set in 2019 (Family)", y = "Height Growth (BLUP)") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=18),
legend.position = "bottom",
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
plot.title = element_blank())
bs_A <- bs_A1 + geom_point(data=gw_19_B_bs_19_P_U,
aes(x = abs(BudSet2019_Plasticity), y = Height_Growth, color = Region,
shape = Region), size=4, alpha = 0.25,
inherit.aes = FALSE)
bs_A
Population plot
gw_19_B_bs_19_P_pop_U <- dplyr::distinct(gw_19_B_bs_19_P_pop, Abs_P,Height_Growth, .keep_all=T)
gw_19_B_bs_19_P_pop_abs_mod <- lm(Height_Growth ~ Abs_P + Region + Abs_P*Region,
data=gw_19_B_bs_19_P_pop_U)
summary(gw_19_B_bs_19_P_pop_abs_mod)##
## Call:
## lm(formula = Height_Growth ~ Abs_P + Region + Abs_P * Region,
## data = gw_19_B_bs_19_P_pop_U)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5522 -0.5923 -0.2067 0.6528 3.7638
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.9029 1.2791 -2.269 0.0269 *
## Abs_P 0.5352 0.4089 1.309 0.1956
## RegionMargin 4.6610 2.7031 1.724 0.0899 .
## RegionEdge 3.6782 2.9341 1.254 0.2149
## Abs_P:RegionMargin -0.4869 0.8089 -0.602 0.5496
## Abs_P:RegionEdge -0.5147 0.9299 -0.554 0.5820
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.139 on 59 degrees of freedom
## Multiple R-squared: 0.5788, Adjusted R-squared: 0.5431
## F-statistic: 16.22 on 5 and 59 DF, p-value: 4.884e-10anova(gw_19_B_bs_19_P_pop_abs_mod)## Analysis of Variance Table
##
## Response: Height_Growth
## Df Sum Sq Mean Sq F value Pr(>F)
## Abs_P 1 9.366 9.366 7.2142 0.009379 **
## Region 2 95.195 47.598 36.6626 4.484e-11 ***
## Abs_P:Region 2 0.711 0.355 0.2738 0.761435
## Residuals 59 76.598 1.298
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Absolute plasticity is significant*.
- Absolute plasticity’s interaction with Core*** and Edge* region is significant***.
# updated fig for the paper
bs_B1 <- sjPlot::plot_model(gw_19_B_bs_19_P_pop_abs_mod, type = "int", colors=reg_palette) +
labs(x = "Plasticity of bud set in 2019 (Population)", y = "") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=18),
legend.position = "bottom",
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
plot.title = element_blank())
bs_B <- bs_B1 + geom_point(data=gw_19_B_bs_19_P_pop_U,
aes(x = abs(BudSet2019_Plasticity), y = Height_Growth, color = Region,
shape = Region), size=4, alpha = 0.25,
inherit.aes = FALSE)
bs_B
Budset 2020 Plasticity
Family data
# blups for height growth
heightGrowth_2020_2 <- heightGrowth_2020
heightGrowth_2020_2 <- heightGrowth_2020_2 %>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)
HeightGrowth_2020_mod <- lmer(Growth ~ Garden + (1|mBed) + (1|Family), heightGrowth_2020_2)
HeightGrowth_2020_blup <- ranef(HeightGrowth_2020_mod)
HeightGrowth_2020_blup <- HeightGrowth_2020_blup$Family
HeightGrowth_2020_blup <- cbind(Family=rownames(HeightGrowth_2020_blup),HeightGrowth_2020_blup)
rownames(HeightGrowth_2020_blup) <- NULL
names(HeightGrowth_2020_blup)[2] <- "Height_Growth"
# plasticity based on budset 2020
gw_20_B_bs_20_P <- merge(x=HeightGrowth_2020_blup,
y=Plasticity[,c("Family","BudSet_2020")],
by="Family")
# rename plasticity column if needed
colnames(gw_20_B_bs_20_P)[colnames(gw_20_B_bs_20_P)=="BudSet_2020"] <- "BudSet2020_Plasticity"
# remove NAs
gw_20_B_bs_20_P <- gw_20_B_bs_20_P[!is.na(gw_20_B_bs_20_P$BudSet2020_Plasticity),]
# Population column creation
gw_20_B_bs_20_P$Population <- gsub("\\_.*", "", gw_20_B_bs_20_P$Family)
# merge climateNA data
gw_20_B_bs_20_P <- merge(x=gw_20_B_bs_20_P,
y=Family_climateNA[,c("Population","PC1","MAT")],
by="Population")
# merge Region data
gw_20_B_bs_20_P <- merge(x=gw_20_B_bs_20_P,
y=heightGrowth_2020[,c("Population","Region")],
by="Population")
# reorder regions
gw_20_B_bs_20_P$Region <- factor(gw_20_B_bs_20_P$Region,levels = c("Core", "Margin", "Edge"))
gw_20_B_bs_20_P <- distinct(gw_20_B_bs_20_P)
gw_20_B_bs_20_P$Abs_P <- abs(gw_20_B_bs_20_P$BudSet2020_Plasticity)Population data
gw_20_B_bs_20_P$Population <- as.factor(gw_20_B_bs_20_P$Population)
gw_20_B_bs_20_P_pop <- gw_20_B_bs_20_P %>%
select(Population, Height_Growth, BudSet2020_Plasticity, Region) %>%
group_by(Population) %>%
dplyr::summarise(Height_Growth = mean(Height_Growth), BudSet2020_Plasticity = mean(BudSet2020_Plasticity))
gw_20_B_bs_20_P_pop <- merge(x=gw_20_B_bs_20_P_pop,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
gw_20_B_bs_20_P_pop <- merge(x=gw_20_B_bs_20_P_pop,
y=heightGrowth_2019[,c("Population","Region")],
by="Population")
# reorder regions
gw_20_B_bs_20_P_pop$Region <- factor(gw_20_B_bs_20_P_pop$Region,levels = c("Core", "Margin", "Edge"))
gw_20_B_bs_20_P_pop <- distinct(gw_20_B_bs_20_P_pop)
gw_20_B_bs_20_P_pop$Abs_P <- abs(gw_20_B_bs_20_P_pop$BudSet2020_Plasticity)Family plot
gw_20_B_bs_20_P_U <- dplyr::distinct(gw_20_B_bs_20_P, Abs_P,Height_Growth, .keep_all=T)
gw_20_B_bs_20_P_fam_abs_mod <- lm(Height_Growth ~ Abs_P + Region + Abs_P*Region,
data=gw_20_B_bs_20_P_U)
summary(gw_20_B_bs_20_P_fam_abs_mod)##
## Call:
## lm(formula = Height_Growth ~ Abs_P + Region + Abs_P * Region,
## data = gw_20_B_bs_20_P_U)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.6269 -1.4662 -0.1439 0.9603 13.2596
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.52439 2.91739 0.523 0.602
## Abs_P -0.39603 0.50069 -0.791 0.430
## RegionMargin 7.87580 6.99085 1.127 0.261
## RegionEdge 0.26813 4.67887 0.057 0.954
## Abs_P:RegionMargin -0.54337 1.17012 -0.464 0.643
## Abs_P:RegionEdge -0.01143 0.80611 -0.014 0.989
##
## Residual standard error: 2.384 on 328 degrees of freedom
## Multiple R-squared: 0.3222, Adjusted R-squared: 0.3119
## F-statistic: 31.19 on 5 and 328 DF, p-value: < 2.2e-16anova(gw_20_B_bs_20_P_fam_abs_mod)## Analysis of Variance Table
##
## Response: Height_Growth
## Df Sum Sq Mean Sq F value Pr(>F)
## Abs_P 1 7.92 7.92 1.3927 0.2388
## Region 2 877.28 438.64 77.1650 <2e-16 ***
## Abs_P:Region 2 1.30 0.65 0.1142 0.8921
## Residuals 328 1864.49 5.68
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Absolute plasticity is significant***.
- Absolute plasticity’s interaction with Core, Margin and Edge region is significant***.
# updated fig for the paper
bs_C1 <- sjPlot::plot_model(gw_20_B_bs_20_P_fam_abs_mod, type = "int", colors=reg_palette) +
labs(x = "Plasticity of bud set in 2020 (Family)", y = "Height Growth (BLUP)") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=18),
legend.position = "bottom",
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
plot.title = element_blank())
bs_C <- bs_C1 + geom_point(data=gw_20_B_bs_20_P_U,
aes(x = abs(BudSet2020_Plasticity), y = Height_Growth, color = Region,
shape = Region), size=4, alpha = 0.25,
inherit.aes = FALSE)
bs_C
Population plot
gw_20_B_bs_20_P_pop_U <- dplyr::distinct(gw_20_B_bs_20_P_pop, Abs_P,Height_Growth, .keep_all=T)
gw_20_B_bs_20_P_pop_abs_mod <- lm(Height_Growth ~ Abs_P + Region + Abs_P*Region,
data=gw_20_B_bs_20_P_pop_U)
summary(gw_20_B_bs_20_P_pop_abs_mod)##
## Call:
## lm(formula = Height_Growth ~ Abs_P + Region + Abs_P * Region,
## data = gw_20_B_bs_20_P_pop_U)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0045 -0.8303 -0.1386 0.6193 9.0547
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.791 11.531 1.283 0.205
## Abs_P -2.669 1.982 -1.347 0.183
## RegionMargin 15.868 22.808 0.696 0.489
## RegionEdge -14.380 15.570 -0.924 0.359
## Abs_P:RegionMargin -1.844 3.824 -0.482 0.631
## Abs_P:RegionEdge 2.507 2.680 0.936 0.353
##
## Residual standard error: 1.644 on 59 degrees of freedom
## Multiple R-squared: 0.4962, Adjusted R-squared: 0.4535
## F-statistic: 11.62 on 5 and 59 DF, p-value: 7.74e-08anova(gw_20_B_bs_20_P_pop_abs_mod)## Analysis of Variance Table
##
## Response: Height_Growth
## Df Sum Sq Mean Sq F value Pr(>F)
## Abs_P 1 7.889 7.889 2.9188 0.09281 .
## Region 2 144.588 72.294 26.7479 5.39e-09 ***
## Abs_P:Region 2 4.607 2.303 0.8523 0.43163
## Residuals 59 159.464 2.703
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Absolute plasticity is significant**.
- Absolute plasticity’s interaction with Core, Margin and Edge region is significant**.
# updated fig for the paper
bs_D1 <- sjPlot::plot_model(gw_20_B_bs_20_P_pop_abs_mod, type = "int", colors=reg_palette) +
labs(x = "Plasticity of bud set in 2020 (Population)", y = "") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=18),
legend.position = "bottom",
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
plot.title = element_blank())
bs_D <- bs_D1 + geom_point(data=gw_20_B_bs_20_P_pop_U,
aes(x = abs(BudSet2020_Plasticity), y = Height_Growth, color = Region,
shape = Region), size=4, alpha = 0.25,
inherit.aes = FALSE)
bs_D
# combined
# create legend
legend_bb <- get_legend(test + theme(legend.position="bottom"))
BS_grid <- cowplot::plot_grid( bs_A + theme(legend.position="none"),
bs_B + theme(legend.position="none"),
bs_C + theme(legend.position="none"),
bs_D + theme(legend.position="none"),
align = 'vh',
labels = c("A", "B", "C", "D"),
hjust = -1,
nrow = 2
)
# bud set final
bs_plot <- cowplot::plot_grid(BS_grid, legend_bb, ncol = 1, rel_heights = c(1,.2))Figure for paper
In text
- Adaptive plasticity within population (Family level)
In supplementary
- Adaptive plasticity among population (Population level)
