Load packages:
library(dplyr)
library(tidyr)
library(lme4)
library(lmerTest)
library(car)
library(emmeans)
library(pbkrtest)
library(parallel)
library(ggplot2)
library(geodist)Read in data:
field<-read.csv('~/Desktop/BCMA archive.R1/BCMA-final-Field-2013-2016-archive.csv', header=T)
surv<-read.csv('~/Desktop/BCMA archive.R1/CFRs-2017surv.csv', header=T)
rp<-read.csv('~/Desktop/BCMA archive.R1/BCMA-final-RefPops-AllMerged.csv', header=T)
dry<-read.csv('~/Desktop/BCMA archive.R1/BCMA-final-DryDownCFR.csv', header=T)
nuc.div<-read.table("~/Desktop/BCMA archive.R1/Bs54.fold4.div.gene.chr.txt", header=T, na="-9")Temporary arrays:
library(pbkrtest)
# Supplementary Table 1A
## subset data
use.a<-subset(field, Expt=="arrays")
## fit model
d4<-lmer(logLAR2 ~ Site + (1|Blk) + Geno + Category + (1|Fam) + Geno*Site + Geno*Category, data=use.a,
control=lmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)), na.action = na.omit)
summary(d4)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logLAR2 ~ Site + (1 | Blk) + Geno + Category + (1 | Fam) + Geno *
## Site + Geno * Category
## Data: use.a
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## REML criterion at convergence: 7778.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5676 -0.6813 -0.1647 0.5183 7.1082
##
## Random effects:
## Groups Name Variance Std.Dev.
## Blk (Intercept) 0.0526425 0.22944
## Fam (Intercept) 0.0004151 0.02037
## Residual 0.2504217 0.50042
## Number of obs: 5193, groups: Blk, 60; Fam, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.82697 0.09276 60.09774 8.915 1.35e-12 ***
## SiteKC -0.19488 0.10978 56.98266 -1.775 0.08120 .
## SiteLA -0.15946 0.10939 56.17866 -1.458 0.15047
## SiteMM -0.20490 0.11009 57.62952 -1.861 0.06783 .
## SiteRM -0.18783 0.10932 56.04588 -1.718 0.09130 .
## SiteUA -0.27554 0.10932 56.03514 -2.521 0.01459 *
## GenoSS 0.15195 0.04653 284.11214 3.266 0.00122 **
## CategoryL 0.10875 0.08107 55.02170 1.341 0.18530
## CategoryM 0.02283 0.07657 57.15397 0.298 0.76662
## SiteKC:GenoSS -0.12412 0.05285 5156.81455 -2.349 0.01888 *
## SiteLA:GenoSS -0.02542 0.05189 5156.57491 -0.490 0.62420
## SiteMM:GenoSS -0.05487 0.05406 5155.93256 -1.015 0.31017
## SiteRM:GenoSS -0.13429 0.05166 5156.76223 -2.600 0.00936 **
## SiteUA:GenoSS -0.10267 0.05156 5157.47298 -1.991 0.04649 *
## GenoSS:CategoryL -0.04208 0.03755 5116.29865 -1.121 0.26249
## GenoSS:CategoryM -0.02329 0.03470 5119.50198 -0.671 0.50214
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## random effects tests
ranova(d4)
## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## logLAR2 ~ Site + Geno + Category + (1 | Blk) + (1 | Fam) + Site:Geno +
## Geno:Category
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 19 -3889.1 7816.2
## (1 | Blk) 18 -4154.2 8344.4 530.21 1 <2e-16 ***
## (1 | Fam) 18 -3890.1 7816.2 1.98 1 0.1597
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## ANOVA with KR ddf
d4.a2<-anova(d4, ddf="Kenward-Roger") # takes a few minutes to run
d4.a2
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Site 2.84488 0.56898 5 52.0 2.2721 0.06077 .
## Geno 2.10329 2.10329 1 9.0 8.3990 0.01767 *
## Category 0.37910 0.18955 2 52.3 0.7569 0.47418
## Site:Geno 2.95790 0.59158 5 5150.0 2.3623 0.03761 *
## Geno:Category 0.31615 0.15807 2 5121.5 0.6312 0.53198
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Figure 2A (left pane)
## extract least squares means
emmeans(d4, pairwise ~ Geno*Site)$emmeans
## Geno Site emmean SE df asymp.LCL asymp.UCL
## LL 401 0.871 0.0790 Inf 0.716 1.026
## SS 401 1.001 0.0778 Inf 0.848 1.154
## LL KC 0.676 0.0775 Inf 0.524 0.828
## SS KC 0.682 0.0771 Inf 0.531 0.833
## LL LA 0.711 0.0770 Inf 0.561 0.862
## SS LA 0.816 0.0770 Inf 0.665 0.967
## LL MM 0.666 0.0779 Inf 0.513 0.819
## SS MM 0.741 0.0775 Inf 0.589 0.893
## LL RM 0.683 0.0769 Inf 0.532 0.834
## SS RM 0.679 0.0769 Inf 0.528 0.830
## LL UA 0.595 0.0769 Inf 0.445 0.746
## SS UA 0.623 0.0769 Inf 0.472 0.773
##
## Results are averaged over the levels of: Category
## Degrees-of-freedom method: asymptotic
## Confidence level used: 0.95
## pairwise comparisons within environments
emmeans(d4, pairwise ~ Geno | Site)$contrasts
## Site = 401:
## contrast estimate SE df z.ratio p.value
## LL - SS -0.13016 0.0414 Inf -3.146 0.0017
##
## Site = KC:
## contrast estimate SE df z.ratio p.value
## LL - SS -0.00604 0.0377 Inf -0.160 0.8726
##
## Site = LA:
## contrast estimate SE df z.ratio p.value
## LL - SS -0.10473 0.0363 Inf -2.885 0.0039
##
## Site = MM:
## contrast estimate SE df z.ratio p.value
## LL - SS -0.07528 0.0393 Inf -1.914 0.0556
##
## Site = RM:
## contrast estimate SE df z.ratio p.value
## LL - SS 0.00414 0.0360 Inf 0.115 0.9084
##
## Site = UA:
## contrast estimate SE df z.ratio p.value
## LL - SS -0.02748 0.0358 Inf -0.767 0.4430
##
## Results are averaged over the levels of: Category
## Degrees-of-freedom method: asymptoticTransplant experiments:
library(pbkrtest)
# Supplementary Table 1B
## Subset data
use.t<-subset(field, Expt=="transplants")
## Fit model
d3<-lmer(logLAR2 ~ Env + (1|Blk) + Geno + (1|Fam) + Geno*Env, data=use.t,
control=lmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
summary(d3)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logLAR2 ~ Env + (1 | Blk) + Geno + (1 | Fam) + Geno * Env
## Data: use.t
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## REML criterion at convergence: 10267.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9669 -0.6982 -0.1579 0.6200 3.8882
##
## Random effects:
## Groups Name Variance Std.Dev.
## Blk (Intercept) 0.15584 0.3948
## Fam (Intercept) 0.01019 0.1009
## Residual 0.88776 0.9422
## Number of obs: 3674, groups: Blk, 149; Fam, 19
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.9044 0.1290 191.5422 14.764 < 2e-16 ***
## EnvALD-2014 -0.6209 0.2019 167.6795 -3.075 0.002460 **
## EnvALD-2015 0.3895 0.1740 227.7899 2.238 0.026206 *
## EnvGTH-2016 -0.8372 0.1806 358.4077 -4.637 4.97e-06 ***
## EnvMAH-2014 -0.7247 0.1951 201.2145 -3.715 0.000263 ***
## EnvMAH-2015 -0.9154 0.1877 200.7535 -4.876 2.20e-06 ***
## EnvSCH-2016 -0.2998 0.1516 201.2438 -1.978 0.049277 *
## EnvSIL-2014 -0.9076 0.1936 188.1999 -4.689 5.26e-06 ***
## EnvSIL-2015 -1.2731 0.1759 187.2644 -7.238 1.14e-11 ***
## GenoSS 0.4168 0.1136 132.8703 3.670 0.000350 ***
## EnvALD-2014:GenoSS -0.3611 0.1464 3538.4064 -2.466 0.013701 *
## EnvALD-2015:GenoSS -0.4212 0.1471 2728.9738 -2.863 0.004225 **
## EnvGTH-2016:GenoSS -0.6919 0.1610 3287.7743 -4.297 1.78e-05 ***
## EnvMAH-2014:GenoSS -0.2475 0.1548 3553.5025 -1.600 0.109797
## EnvMAH-2015:GenoSS -0.6139 0.1457 2600.6488 -4.212 2.61e-05 ***
## EnvSCH-2016:GenoSS -0.4330 0.1190 2314.0075 -3.637 0.000282 ***
## EnvSIL-2014:GenoSS -0.4764 0.1523 3551.2107 -3.129 0.001768 **
## EnvSIL-2015:GenoSS -0.6775 0.1408 2489.5631 -4.813 1.58e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## random effects tests
ranova(d3)
## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## logLAR2 ~ Env + Geno + (1 | Blk) + (1 | Fam) + Env:Geno
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 21 -5133.7 10309
## (1 | Blk) 20 -5312.2 10664 356.91 1 < 2.2e-16 ***
## (1 | Fam) 20 -5142.1 10324 16.80 1 4.148e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## ANOVA with KR ddf
d3.a2<-anova(d3, ddf="Kenward-Roger") # takes ~2.5 min to run
d3.a2
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Env 199.208 24.9010 8 135.15 28.0462 < 2.2e-16 ***
## Geno 0.083 0.0830 1 15.96 0.0935 0.7638
## Env:Geno 29.754 3.7193 8 2840.39 4.1884 5.333e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Figure 2A (center pane)
## extract least squares means
emmeans(d3, pairwise ~ Geno* Env)$emmeans
## Geno Env emmean SE df asymp.LCL asymp.UCL
## LL ALD-2013 1.904 0.1290 Inf 1.652 2.157
## SS ALD-2013 2.321 0.1260 Inf 2.074 2.568
## LL ALD-2014 1.284 0.1665 Inf 0.957 1.610
## SS ALD-2014 1.339 0.1630 Inf 1.020 1.659
## LL ALD-2015 2.294 0.1264 Inf 2.046 2.542
## SS ALD-2015 2.289 0.1221 Inf 2.050 2.529
## LL GTH-2016 1.067 0.1371 Inf 0.799 1.336
## SS GTH-2016 0.792 0.1117 Inf 0.573 1.011
## LL MAH-2014 1.180 0.1580 Inf 0.870 1.489
## SS MAH-2014 1.349 0.1543 Inf 1.047 1.652
## LL MAH-2015 0.989 0.1443 Inf 0.706 1.272
## SS MAH-2015 0.792 0.1395 Inf 0.519 1.065
## LL SCH-2016 1.605 0.0948 Inf 1.419 1.790
## SS SCH-2016 1.588 0.0905 Inf 1.411 1.766
## LL SIL-2014 0.997 0.1562 Inf 0.691 1.303
## SS SIL-2014 0.937 0.1542 Inf 0.635 1.239
## LL SIL-2015 0.631 0.1287 Inf 0.379 0.884
## SS SIL-2015 0.371 0.1312 Inf 0.113 0.628
##
## Degrees-of-freedom method: asymptotic
## Confidence level used: 0.95
## pairwise comparisons within environments
emmeans(d3, pairwise ~ Geno | Env)$contrasts
## Env = ALD-2013:
## contrast estimate SE df z.ratio p.value
## LL - SS -0.41682 0.114 Inf -3.670 0.0002
##
## Env = ALD-2014:
## contrast estimate SE df z.ratio p.value
## LL - SS -0.05573 0.123 Inf -0.452 0.6512
##
## Env = ALD-2015:
## contrast estimate SE df z.ratio p.value
## LL - SS 0.00442 0.115 Inf 0.039 0.9693
##
## Env = GTH-2016:
## contrast estimate SE df z.ratio p.value
## LL - SS 0.27507 0.136 Inf 2.016 0.0438
##
## Env = MAH-2014:
## contrast estimate SE df z.ratio p.value
## LL - SS -0.16928 0.133 Inf -1.273 0.2030
##
## Env = MAH-2015:
## contrast estimate SE df z.ratio p.value
## LL - SS 0.19711 0.113 Inf 1.749 0.0803
##
## Env = SCH-2016:
## contrast estimate SE df z.ratio p.value
## LL - SS 0.01618 0.082 Inf 0.197 0.8436
##
## Env = SIL-2014:
## contrast estimate SE df z.ratio p.value
## LL - SS 0.05963 0.130 Inf 0.458 0.6468
##
## Env = SIL-2015:
## contrast estimate SE df z.ratio p.value
## LL - SS 0.26066 0.106 Inf 2.455 0.0141
##
## Degrees-of-freedom method: asymptoticLaboratory assay:
library(pbkrtest)
# Supplementary Table 1C
## subset data
use.l<-subset(field, Expt=="Lab_Tni")
## fit model
d5<-lmer(logLAR2 ~ (1|Blk) + Geno + (1|Fam), data=use.l,
control=lmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
summary(d5)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logLAR2 ~ (1 | Blk) + Geno + (1 | Fam)
## Data: use.l
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## REML criterion at convergence: 681.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6669 -0.4081 0.2268 0.6597 1.6656
##
## Random effects:
## Groups Name Variance Std.Dev.
## Fam (Intercept) 0.0000 0.0000
## Blk (Intercept) 0.1957 0.4424
## Residual 1.3450 1.1598
## Number of obs: 214, groups: Fam, 9; Blk, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.5028 0.2296 5.3897 15.254 1.22e-05 ***
## GenoSS -0.9513 0.1590 208.0067 -5.982 9.48e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## GenoSS -0.372
## convergence code: 0
## boundary (singular) fit: see ?isSingular
## random effects tests
ranova(d5)
## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## logLAR2 ~ Geno + (1 | Blk) + (1 | Fam)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 5 -340.86 691.72
## (1 | Blk) 4 -348.70 705.40 15.683 1 7.49e-05 ***
## (1 | Fam) 4 -340.86 689.72 0.000 1 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## ANOVA with KR ddf
d5.a2<-anova(d5, ddf="Kenward-Roger")
d5.a2
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Geno 48.114 48.114 1 6.8801 35.772 0.00059 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Figure 2A (right pane)
## extract least squares means
emmeans(d5, pairwise ~ Geno)$emmeans
## Geno emmean SE df lower.CL upper.CL
## LL 3.50 0.230 5.10 2.92 4.09
## SS 2.55 0.226 4.84 1.97 3.14
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## pairwise comparisons within environments
emmeans(d5, pairwise ~ Geno)$contrasts
## contrast estimate SE df t.ratio p.value
## LL - SS 0.951 0.159 6.88 5.981 0.0006
##
## Degrees-of-freedom method: kenward-rogerEffect of damage on the probability of reproduction:
library(car)
# Supplementary Table 2A
## filter field data to vernalized adult transplants
adult<-filter(field, Cohort != "s2013" & Expt=="transplants" & CensusYr < 2017)
## fit model
prrep2<-glmer(didrep ~ Env*LAR2 + (1|Blk), data=adult, family=binomial(link="logit"),
control=glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
summary(prrep2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: didrep ~ Env * LAR2 + (1 | Blk)
## Data: adult
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## AIC BIC logLik deviance df.resid
## 3144.7 3247.3 -1555.3 3110.7 3077
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.6727 -0.7020 0.2494 0.5977 5.4780
##
## Random effects:
## Groups Name Variance Std.Dev.
## Blk (Intercept) 0.6283 0.7927
## Number of obs: 3094, groups: Blk, 132
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.739506 0.318324 2.323 0.02017 *
## EnvALD-2015 -0.305176 0.411330 -0.742 0.45813
## EnvGTH-2016 -0.527352 0.389808 -1.353 0.17610
## EnvMAH-2014 0.733268 0.448010 1.637 0.10169
## EnvMAH-2015 -0.378877 0.425719 -0.890 0.37348
## EnvSCH-2016 2.460741 0.397835 6.185 6.20e-10 ***
## EnvSIL-2014 -1.174447 0.436319 -2.692 0.00711 **
## EnvSIL-2015 -1.151725 0.400879 -2.873 0.00407 **
## LAR2 -0.033620 0.013593 -2.473 0.01339 *
## EnvALD-2015:LAR2 0.006939 0.016984 0.409 0.68284
## EnvGTH-2016:LAR2 0.004065 0.035227 0.115 0.90813
## EnvMAH-2014:LAR2 0.017932 0.018926 0.947 0.34341
## EnvMAH-2015:LAR2 -0.010229 0.021379 -0.478 0.63234
## EnvSCH-2016:LAR2 -0.026026 0.018573 -1.401 0.16112
## EnvSIL-2014:LAR2 -0.073384 0.043195 -1.699 0.08934 .
## EnvSIL-2015:LAR2 0.153389 0.034531 4.442 8.91e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(prrep2, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: didrep
## Chisq Df Pr(>Chisq)
## (Intercept) 5.3969 1 0.02017 *
## Env 157.1284 7 < 2.2e-16 ***
## LAR2 6.1172 1 0.01339 *
## Env:LAR2 33.2727 7 2.355e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(prrep2)
## [1] 3094
## random effect test
prrep2b<-glm(didrep ~ Env*LAR2, data=adult, family=binomial(link="logit"))
anova(prrep2, prrep2b)
## Data: adult
## Models:
## prrep2b: didrep ~ Env * LAR2
## prrep2: didrep ~ Env * LAR2 + (1 | Blk)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## prrep2b 16 3281.1 3377.7 -1624.5 3249.1
## prrep2 17 3144.7 3247.3 -1555.3 3110.7 138.37 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Effect of damage on total fruit size:
library(car)
# Supplementary Table 2B
## filter field data to plants in 2016 Schofield Garden (best replication) that did reproduce
use.sch<-filter(adult, Env=="SCH-2016" & didrep==1 & !is.na(TotalRepro) & TotalRepro > 0)
## fit model
repamt2<-lmer(log(TotalRepro) ~ LAR2 + (1|Blk), data=use.sch)
summary(repamt2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(TotalRepro) ~ LAR2 + (1 | Blk)
## Data: use.sch
##
## REML criterion at convergence: 1696.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9464 -0.5097 0.1398 0.6902 2.2962
##
## Random effects:
## Groups Name Variance Std.Dev.
## Blk (Intercept) 0.03459 0.1860
## Residual 0.40868 0.6393
## Number of obs: 848, groups: Blk, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.435306 0.044865 42.236673 54.281 <2e-16 ***
## LAR2 -0.007115 0.003025 845.345277 -2.352 0.0189 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## LAR2 -0.415
Anova(repamt2, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: log(TotalRepro)
## Chisq Df Pr(>Chisq)
## (Intercept) 2946.4000 1 < 2e-16 ***
## LAR2 5.5318 1 0.01867 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(repamt2)
## [1] 848
ranova(repamt2)
## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log(TotalRepro) ~ LAR2 + (1 | Blk)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 4 -848.04 1704.1
## (1 | Blk) 3 -864.20 1734.4 32.327 1 1.303e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Effect of damage on survival:
library(car)
# Supplementary Table 4
s1<-glmer(X2017_Surv0 ~ X2016_LAR2 + (1|Blk), data=surv, family=binomial(link="logit"))
summary(s1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: X2017_Surv0 ~ X2016_LAR2 + (1 | Blk)
## Data: surv
##
## AIC BIC logLik deviance df.resid
## 686.5 701.8 -340.2 680.5 1216
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.0877 0.1869 0.2483 0.2983 0.9724
##
## Random effects:
## Groups Name Variance Std.Dev.
## Blk (Intercept) 0.8843 0.9404
## Number of obs: 1219, groups: Blk, 30
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.76623 0.23575 11.734 <2e-16 ***
## X2016_LAR2 -0.01366 0.01199 -1.139 0.255
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## X2016_LAR2 -0.385
Anova(s1, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: X2017_Surv0
## Chisq Df Pr(>Chisq)
## (Intercept) 137.6811 1 <2e-16 ***
## X2016_LAR2 1.2975 1 0.2547
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(s1)
## [1] 1219
## random effect test
s1b<-glm(X2017_Surv0 ~ X2016_LAR2, data=surv, family=binomial(link="logit"))
anova(s1, s1b)
## Data: surv
## Models:
## s1b: X2017_Surv0 ~ X2016_LAR2
## s1: X2017_Surv0 ~ X2016_LAR2 + (1 | Blk)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## s1b 2 737.21 747.43 -366.61 733.21
## s1 3 686.49 701.81 -340.25 680.49 52.721 1 3.845e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Temporary arrays:
library(car)
# Supplementary Table 3A
## define array dataset
a<-use.a
## fit full model
s2<-glmer(Surv2 ~ Site + Geno + Category + Geno*Site + Geno*Category + (1|Blk) + (1|Fam), data=a,
family=binomial(link="logit"),
control=glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
summary(s2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Surv2 ~ Site + Geno + Category + Geno * Site + Geno * Category +
## (1 | Blk) + (1 | Fam)
## Data: a
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## AIC BIC logLik deviance df.resid
## 3246.8 3367.0 -1605.4 3210.8 5862
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -11.5696 0.0990 0.1702 0.3416 2.3645
##
## Random effects:
## Groups Name Variance Std.Dev.
## Blk (Intercept) 0.8219 0.9066
## Fam (Intercept) 0.1814 0.4259
## Number of obs: 5880, groups: Blk, 60; Fam, 10
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.11661 0.43191 2.585 0.009730 **
## SiteKC 1.10565 0.44723 2.472 0.013427 *
## SiteLA 2.62562 0.48901 5.369 7.91e-08 ***
## SiteMM 0.32161 0.43655 0.737 0.461299
## SiteRM 4.17569 0.64793 6.445 1.16e-10 ***
## SiteUA 3.72388 0.58413 6.375 1.83e-10 ***
## GenoSS 0.82496 0.35940 2.295 0.021709 *
## CategoryL -0.12367 0.37978 -0.326 0.744701
## CategoryM -0.43696 0.35515 -1.230 0.218566
## SiteKC:GenoSS -0.03145 0.28434 -0.111 0.911929
## SiteLA:GenoSS -1.12606 0.35326 -3.188 0.001434 **
## SiteMM:GenoSS -0.42347 0.24258 -1.746 0.080858 .
## SiteRM:GenoSS -2.18724 0.57477 -3.805 0.000142 ***
## SiteUA:GenoSS -1.19828 0.53969 -2.220 0.026398 *
## GenoSS:CategoryL 0.16940 0.26739 0.634 0.526381
## GenoSS:CategoryM 0.15592 0.23746 0.657 0.511431
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(s2)
## [1] 5880
## Wald test
Anova(s2, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: Surv2
## Chisq Df Pr(>Chisq)
## (Intercept) 6.6837 1 0.00973 **
## Site 89.2772 5 < 2.2e-16 ***
## Geno 5.2690 1 0.02171 *
## Category 1.7201 2 0.42313
## Site:Geno 26.6144 5 6.78e-05 ***
## Geno:Category 0.5318 2 0.76653
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## random effects tests
s2b<-glmer(Surv2 ~ Site + Geno + Category + Geno*Site + Geno*Category + (1|Fam), data=a,
family=binomial(link="logit"),
control=glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
s2c<-glmer(Surv2 ~ Site + Geno + Category + Geno*Site + Geno*Category + (1|Blk), data=a,
family=binomial(link="logit"),
control=glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
anova(s2, s2b)
## Data: a
## Models:
## s2b: Surv2 ~ Site + Geno + Category + Geno * Site + Geno * Category +
## s2b: (1 | Fam)
## s2: Surv2 ~ Site + Geno + Category + Geno * Site + Geno * Category +
## s2: (1 | Blk) + (1 | Fam)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## s2b 17 3510.2 3623.7 -1738.1 3476.2
## s2 18 3246.8 3367.0 -1605.4 3210.8 265.42 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(s2, s2c)
## Data: a
## Models:
## s2c: Surv2 ~ Site + Geno + Category + Geno * Site + Geno * Category +
## s2c: (1 | Blk)
## s2: Surv2 ~ Site + Geno + Category + Geno * Site + Geno * Category +
## s2: (1 | Blk) + (1 | Fam)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## s2c 17 3308.5 3422 -1637.2 3274.5
## s2 18 3246.8 3367 -1605.4 3210.8 63.704 1 1.446e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## parametric bootstrap
# parametric bootstrap#
# method is very slow,#
# i.e. hours, so they #
# are commented out. ##
# delete #s to run it #
### define reduced models eliminating effects
#### drop geno*site interaction
s2.gs<-glmer(Surv2 ~ Site + Geno + Category + Geno*Category + (1|Blk) + (1|Fam), data=a,
family=binomial(link="logit"),
control=glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
#### drop genotype (main effect + interactions)
s2.g<-glmer(Surv2 ~ Site + Category + (1|Blk) + (1|Fam), data=a,
family=binomial(link="logit"),
control=glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
#### drop site (main effect + interactions)
s2.s<-glmer(Surv2 ~ Geno + Category + Geno*Category + (1|Blk) + (1|Fam), data=a,
family=binomial(link="logit"),
control=glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
#### drop geno*category interaction
s2.gc<-glmer(Surv2 ~ Site + Geno + Category + Geno*Site + (1|Blk) + (1|Fam), data=a,
family=binomial(link="logit"),
control=glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
#### drop category (main effect + interactions)
s2.c<-glmer(Surv2 ~ Site + Geno + Geno*Site +(1|Blk) + (1|Fam), data=a,
family=binomial(link="logit"),
control=glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
### run bootstrapping in parallel
library(parallel)
nc<-detectCores()
cl<-makeCluster(rep("localhost", nc))
#### test geno*site
# s2.pb.gs<-PBmodcomp(s2, s2.gs, cl=cl)
# s2.pb.gs
#### test geno
# s2.pb.g<-PBmodcomp(s2, s2.g, cl=cl)
# s2.pb.g
#### test site
# s2.pb.s<-PBmodcomp(s2, s2.s, cl=cl)
# s2.pb.s
#### test category
# s2.pb.c<-PBmodcomp(s2, s2.c, cl=cl)
# s2.pb.c
#### test geno*category
# s2.pb.gc<-PBmodcomp(s2, s2.gc, cl=cl)
# s2.pb.gc
### turn off cluster
stopCluster(cl)
# Figure 2C (left pane)
## extract least squares means
emmeans(s2, pairwise ~ Geno* Site, type="response")$emmeans
## Geno Site prob SE df asymp.LCL asymp.UCL
## LL 401 0.717 0.07394 Inf 0.554 0.838
## SS 401 0.866 0.04337 Inf 0.756 0.930
## LL KC 0.884 0.03848 Inf 0.785 0.941
## SS KC 0.950 0.01891 Inf 0.897 0.976
## LL LA 0.972 0.01152 Inf 0.938 0.988
## SS LA 0.967 0.01343 Inf 0.928 0.985
## LL MM 0.778 0.06282 Inf 0.632 0.877
## SS MM 0.853 0.04616 Inf 0.739 0.923
## LL RM 0.994 0.00361 Inf 0.981 0.998
## SS RM 0.979 0.00905 Inf 0.952 0.991
## LL UA 0.991 0.00497 Inf 0.974 0.997
## SS UA 0.988 0.00590 Inf 0.969 0.995
##
## Results are averaged over the levels of: Category
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
## pairwise comparisons within environments
emmeans(s2, pairwise ~ Geno | Site, type="response")$contrasts
## Site = 401:
## contrast odds.ratio SE df z.ratio p.value
## LL / SS 0.393 0.127 Inf -2.896 0.0038
##
## Site = KC:
## contrast odds.ratio SE df z.ratio p.value
## LL / SS 0.406 0.143 Inf -2.567 0.0103
##
## Site = LA:
## contrast odds.ratio SE df z.ratio p.value
## LL / SS 1.212 0.497 Inf 0.470 0.6384
##
## Site = MM:
## contrast odds.ratio SE df z.ratio p.value
## LL / SS 0.601 0.192 Inf -1.592 0.1115
##
## Site = RM:
## contrast odds.ratio SE df z.ratio p.value
## LL / SS 3.504 2.147 Inf 2.046 0.0407
##
## Site = UA:
## contrast odds.ratio SE df z.ratio p.value
## LL / SS 1.303 0.755 Inf 0.457 0.6474
##
## Results are averaged over the levels of: Category
## Tests are performed on the log odds ratio scaleTransplant experiments:
library(car)
# Supplementary Table 3B
## define transplant dataset
t<-use.t
## fit full model
s3<-glmer(Surv2 ~ Env + Geno + Geno*Env + (1|Blk) + (1|Fam), data=use.t, family=binomial(link="logit"), control=glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
summary(s3)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Surv2 ~ Env + Geno + Geno * Env + (1 | Blk) + (1 | Fam)
## Data: use.t
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## AIC BIC logLik deviance df.resid
## 7134.1 7270.7 -3547.0 7094.1 6840
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.7314 -0.6663 0.2340 0.6305 4.9096
##
## Random effects:
## Groups Name Variance Std.Dev.
## Blk (Intercept) 1.0591 1.0291
## Fam (Intercept) 0.0546 0.2337
## Number of obs: 6860, groups: Blk, 150; Fam, 19
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.63717 0.30195 2.110 0.034847 *
## EnvALD-2014 1.60151 0.53629 2.986 0.002824 **
## EnvALD-2015 -1.16282 0.39404 -2.951 0.003168 **
## EnvGTH-2016 -2.70750 0.37217 -7.275 3.47e-13 ***
## EnvMAH-2014 -0.94173 0.43968 -2.142 0.032206 *
## EnvMAH-2015 -0.01805 0.44683 -0.040 0.967778
## EnvSCH-2016 1.24843 0.37153 3.360 0.000779 ***
## EnvSIL-2014 -0.15597 0.45916 -0.340 0.734089
## EnvSIL-2015 0.61895 0.42558 1.454 0.145847
## GenoSS -0.26820 0.20970 -1.279 0.200893
## EnvALD-2014:GenoSS -0.07596 0.35969 -0.211 0.832754
## EnvALD-2015:GenoSS 0.34721 0.23224 1.495 0.134895
## EnvGTH-2016:GenoSS 1.05173 0.23880 4.404 1.06e-05 ***
## EnvMAH-2014:GenoSS 0.10735 0.25520 0.421 0.674017
## EnvMAH-2015:GenoSS 0.50360 0.28238 1.783 0.074516 .
## EnvSCH-2016:GenoSS 0.75182 0.24825 3.028 0.002458 **
## EnvSIL-2014:GenoSS -0.05341 0.26036 -0.205 0.837475
## EnvSIL-2015:GenoSS -0.68525 0.25734 -2.663 0.007749 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(s3)
## [1] 6860
## random effects tests
s3b<-glmer(Surv2 ~ Env + Geno + Geno*Env + (1|Fam), data=use.t, family=binomial(link="logit"), control=glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
s3c<-glmer(Surv2 ~ Env + Geno + Geno*Env + (1|Blk), data=use.t, family=binomial(link="logit"), control=glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
anova(s3, s3b)
## Data: use.t
## Models:
## s3b: Surv2 ~ Env + Geno + Geno * Env + (1 | Fam)
## s3: Surv2 ~ Env + Geno + Geno * Env + (1 | Blk) + (1 | Fam)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## s3b 19 7847.7 7977.5 -3904.8 7809.7
## s3 20 7134.1 7270.7 -3547.0 7094.1 715.59 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(s3, s3c)
## Data: use.t
## Models:
## s3c: Surv2 ~ Env + Geno + Geno * Env + (1 | Blk)
## s3: Surv2 ~ Env + Geno + Geno * Env + (1 | Blk) + (1 | Fam)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## s3c 19 7159.1 7288.9 -3560.5 7121.1
## s3 20 7134.1 7270.7 -3547.0 7094.1 27.033 1 2e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Wald test
Anova(s3, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: Surv2
## Chisq Df Pr(>Chisq)
## (Intercept) 4.4527 1 0.03485 *
## Env 198.1527 8 < 2.2e-16 ***
## Geno 1.6359 1 0.20089
## Env:Geno 61.4704 8 2.397e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## parametric bootstrap
# parametric bootstrap#
# method is very slow,#
# i.e. hours, so they #
# are commented out. ##
# delete #s to run it #
### define reduced models eliminating effects
#### drop geno*site interaction
s3.gs<-glmer(Surv2 ~ Env + Geno + (1|Blk) + (1|Fam), data=use.t, family=binomial(link="logit"), control=glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
#### drop genotype (main effect + interaction)
s3.g<-glmer(Surv2 ~ Env + (1|Blk) + (1|Fam), data=use.t, family=binomial(link="logit"), control=glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
#### drop site (main effect + interaction)
s3.s<-glmer(Surv2 ~ Geno + (1|Blk) + (1|Fam), data=use.t, family=binomial(link="logit"), control=glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
### run bootstrapping in parallel
nc<-detectCores()
cl<-makeCluster(rep("localhost", nc))
#### test geno*site
# s3.pb.gs<-PBmodcomp(s3, s3.gs, cl=cl)
# s3.pb.gs
#### test effect of geno
# s3.pb.g<-PBmodcomp(s3, s3.g, cl=cl)
# s3.pb.g
#### test effect of site
# s3.pb.s<-PBmodcomp(s3, s3.s, cl=cl)
# s3.pb.s
### turn off cluster
stopCluster(cl)
# Figure 2C (right pane)
## extract least squares means
emmeans(s3, pairwise ~ Geno* Env, type="response")$emmeans
## Geno Env prob SE df asymp.LCL asymp.UCL
## LL ALD-2013 0.654 0.0683 Inf 0.5113 0.774
## SS ALD-2013 0.591 0.0709 Inf 0.4488 0.720
## LL ALD-2014 0.904 0.0404 Inf 0.7907 0.959
## SS ALD-2014 0.869 0.0495 Inf 0.7391 0.940
## LL ALD-2015 0.372 0.0645 Inf 0.2559 0.504
## SS ALD-2015 0.390 0.0650 Inf 0.2726 0.522
## LL GTH-2016 0.112 0.0246 Inf 0.0721 0.170
## SS GTH-2016 0.216 0.0388 Inf 0.1499 0.302
## LL MAH-2014 0.424 0.0848 Inf 0.2719 0.593
## SS MAH-2014 0.386 0.0804 Inf 0.2441 0.550
## LL MAH-2015 0.650 0.0790 Inf 0.4846 0.786
## SS MAH-2015 0.702 0.0722 Inf 0.5446 0.822
## LL SCH-2016 0.868 0.0279 Inf 0.8033 0.914
## SS SCH-2016 0.914 0.0193 Inf 0.8682 0.945
## LL SIL-2014 0.618 0.0877 Inf 0.4385 0.770
## SS SIL-2014 0.540 0.0899 Inf 0.3660 0.704
## LL SIL-2015 0.778 0.0551 Inf 0.6525 0.868
## SS SIL-2015 0.575 0.0754 Inf 0.4250 0.712
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
## pairwise contrasts within environments
emmeans(s3, pairwise ~ Geno | Env, type="response")$contrasts
## Env = ALD-2013:
## contrast odds.ratio SE df z.ratio p.value
## LL / SS 1.308 0.274 Inf 1.279 0.2009
##
## Env = ALD-2014:
## contrast odds.ratio SE df z.ratio p.value
## LL / SS 1.411 0.485 Inf 1.000 0.3172
##
## Env = ALD-2015:
## contrast odds.ratio SE df z.ratio p.value
## LL / SS 0.924 0.168 Inf -0.434 0.6643
##
## Env = GTH-2016:
## contrast odds.ratio SE df z.ratio p.value
## LL / SS 0.457 0.092 Inf -3.889 0.0001
##
## Env = MAH-2014:
## contrast odds.ratio SE df z.ratio p.value
## LL / SS 1.175 0.273 Inf 0.692 0.4889
##
## Env = MAH-2015:
## contrast odds.ratio SE df z.ratio p.value
## LL / SS 0.790 0.191 Inf -0.972 0.3311
##
## Env = SCH-2016:
## contrast odds.ratio SE df z.ratio p.value
## LL / SS 0.617 0.130 Inf -2.293 0.0219
##
## Env = SIL-2014:
## contrast odds.ratio SE df z.ratio p.value
## LL / SS 1.379 0.328 Inf 1.351 0.1767
##
## Env = SIL-2015:
## contrast odds.ratio SE df z.ratio p.value
## LL / SS 2.595 0.551 Inf 4.491 <.0001
##
## Tests are performed on the log odds ratio scaleThese analyses (Supplementary Table 5) were conducted using standard least squares linear models in JMP.
Genetic variation in morphology in response to drought:
library(car)
dry$Rack<-as.character(dry$Rack)
dry$Fam<-as.character(dry$Fam)
# Supplementary Table 6
## separate statistical models for each morphological response
### leaf length 2 (final height)
len<-lmer(LL_Length2_mm ~ Geno + Trmt + Geno*Trmt + LL1_mm + (1|Rack) + (1|Fam), data=dry, control=lmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
summary(len)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LL_Length2_mm ~ Geno + Trmt + Geno * Trmt + LL1_mm + (1 | Rack) +
## (1 | Fam)
## Data: dry
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## REML criterion at convergence: 7363.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5996 -0.4916 0.0870 0.6153 3.5471
##
## Random effects:
## Groups Name Variance Std.Dev.
## Rack (Intercept) 23.040 4.80
## Fam (Intercept) 6.551 2.56
## Residual 122.681 11.08
## Number of obs: 954, groups: Rack, 32; Fam, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 28.96138 1.89243 36.53828 15.304 < 2e-16 ***
## GenoMM -1.70151 1.91204 10.66410 -0.890 0.39316
## TrmtWet -3.15346 1.04385 922.77176 -3.021 0.00259 **
## LL1_mm 0.53371 0.03178 926.63889 16.794 < 2e-16 ***
## GenoMM:TrmtWet 7.73979 1.48079 917.95306 5.227 2.14e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) GenoMM TrmtWt LL1_mm
## GenoMM -0.499
## TrmtWet -0.181 0.261
## LL1_mm -0.537 -0.012 -0.151
## GnMM:TrmtWt 0.083 -0.367 -0.711 0.189
Anova(len, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: LL_Length2_mm
## Chisq Df Pr(>Chisq)
## (Intercept) 234.2056 1 < 2.2e-16 ***
## Geno 0.7919 1 0.373525
## Trmt 9.1264 1 0.002519 **
## LL1_mm 282.0298 1 < 2.2e-16 ***
## Geno:Trmt 27.3194 1 1.725e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(len)
## [1] 954
ranova(len)
## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## LL_Length2_mm ~ Geno + Trmt + LL1_mm + (1 | Rack) + (1 | Fam) +
## Geno:Trmt
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 8 -3681.6 7379.3
## (1 | Rack) 7 -3727.6 7469.2 91.868 1 < 2.2e-16 ***
## (1 | Fam) 7 -3693.5 7400.9 23.615 1 1.177e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### total leaves
tlf<-lmer(TotalLvs ~ Geno + Trmt + Geno*Trmt + (1|Rack) + (1|Fam) + LL1_mm, data=dry, control=lmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
summary(tlf)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: TotalLvs ~ Geno + Trmt + Geno * Trmt + (1 | Rack) + (1 | Fam) +
## LL1_mm
## Data: dry
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## REML criterion at convergence: 4667.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5549 -0.5734 -0.0095 0.4870 12.1697
##
## Random effects:
## Groups Name Variance Std.Dev.
## Rack (Intercept) 0.38291 0.6188
## Fam (Intercept) 0.04476 0.2116
## Residual 7.32034 2.7056
## Number of obs: 958, groups: Rack, 32; Fam, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 10.384840 0.326204 98.442935 31.835 <2e-16 ***
## GenoMM -0.272180 0.281167 19.964449 -0.968 0.3446
## TrmtWet 0.505336 0.252187 941.935067 2.004 0.0454 *
## LL1_mm 0.150553 0.007403 719.815659 20.338 <2e-16 ***
## GenoMM:TrmtWet 1.048174 0.359182 937.492029 2.918 0.0036 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) GenoMM TrmtWt LL1_mm
## GenoMM -0.409
## TrmtWet -0.269 0.433
## LL1_mm -0.719 -0.029 -0.142
## GnMM:TrmtWt 0.130 -0.610 -0.710 0.182
Anova(tlf, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: TotalLvs
## Chisq Df Pr(>Chisq)
## (Intercept) 1013.4948 1 < 2e-16 ***
## Geno 0.9371 1 0.33303
## Trmt 4.0153 1 0.04509 *
## LL1_mm 413.6266 1 < 2e-16 ***
## Geno:Trmt 8.5160 1 0.00352 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(tlf)
## [1] 958
ranova(tlf)
## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## TotalLvs ~ Geno + Trmt + LL1_mm + (1 | Rack) + (1 | Fam) + Geno:Trmt
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 8 -2333.9 4683.9
## (1 | Rack) 7 -2342.7 4699.3 17.4105 1 3.012e-05 ***
## (1 | Fam) 7 -2334.4 4682.8 0.9143 1 0.339
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### LMA
lma<-lmer(LMAest ~ Geno + Trmt + Geno*Trmt + (1|Rack)+ (1|Fam) + LL1_mm, data=dry, control=lmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
summary(lma)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LMAest ~ Geno + Trmt + Geno * Trmt + (1 | Rack) + (1 | Fam) +
## LL1_mm
## Data: dry
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## REML criterion at convergence: -3393.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1426 -0.3309 -0.1003 0.1852 13.6523
##
## Random effects:
## Groups Name Variance Std.Dev.
## Rack (Intercept) 4.061e-06 0.002015
## Fam (Intercept) 1.103e-06 0.001050
## Residual 1.300e-04 0.011401
## Number of obs: 568, groups: Rack, 32; Fam, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.979e-02 1.559e-03 6.274e+01 12.694 < 2e-16 ***
## GenoMM -8.462e-04 1.244e-03 9.808e+00 -0.680 0.51226
## TrmtWet 3.215e-03 1.852e-03 4.252e+02 1.736 0.08327 .
## LL1_mm 1.053e-04 3.876e-05 4.256e+02 2.716 0.00688 **
## GenoMM:TrmtWet 3.319e-03 2.704e-03 5.472e+02 1.227 0.22025
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) GenoMM TrmtWt LL1_mm
## GenoMM -0.366
## TrmtWet -0.104 0.247
## LL1_mm -0.794 -0.041 -0.109
## GnMM:TrmtWt 0.036 -0.334 -0.639 0.122
Anova(lma, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: LMAest
## Chisq Df Pr(>Chisq)
## (Intercept) 161.1296 1 < 2.2e-16 ***
## Geno 0.4624 1 0.496526
## Trmt 3.0141 1 0.082545 .
## LL1_mm 7.3763 1 0.006609 **
## Geno:Trmt 1.5062 1 0.219726
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(lma)
## [1] 568
ranova(lma)
## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## LMAest ~ Geno + Trmt + LL1_mm + (1 | Rack) + (1 | Fam) + Geno:Trmt
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 8 1697.0 -3377.9
## (1 | Rack) 7 1695.3 -3376.7 3.1922 1 0.07399 .
## (1 | Fam) 7 1696.6 -3379.2 0.6450 1 0.42191
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### leaf width
wid<-lmer(LL_Width2_mm ~ Geno + Trmt + Geno*Trmt + LL1_mm + (1|Rack) + (1|Fam), data=dry, control=lmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
summary(wid)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LL_Width2_mm ~ Geno + Trmt + Geno * Trmt + LL1_mm + (1 | Rack) +
## (1 | Fam)
## Data: dry
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## REML criterion at convergence: 3762.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8023 -0.5186 0.0553 0.6007 3.0446
##
## Random effects:
## Groups Name Variance Std.Dev.
## Rack (Intercept) 0.5363 0.7323
## Fam (Intercept) 0.1952 0.4418
## Residual 2.7882 1.6698
## Number of obs: 951, groups: Rack, 32; Fam, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 6.339402 0.302091 30.180346 20.985 < 2e-16 ***
## GenoMM -0.250119 0.318963 10.134307 -0.784 0.4509
## TrmtWet -0.390530 0.157776 918.670530 -2.475 0.0135 *
## LL1_mm 0.061854 0.004812 928.831192 12.855 < 2e-16 ***
## GenoMM:TrmtWet 1.128147 0.223651 913.548309 5.044 5.49e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) GenoMM TrmtWt LL1_mm
## GenoMM -0.522
## TrmtWet -0.171 0.237
## LL1_mm -0.508 -0.012 -0.154
## GnMM:TrmtWt 0.079 -0.333 -0.712 0.189
Anova(wid, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: LL_Width2_mm
## Chisq Df Pr(>Chisq)
## (Intercept) 440.3730 1 < 2.2e-16 ***
## Geno 0.6149 1 0.43294
## Trmt 6.1267 1 0.01332 *
## LL1_mm 165.2502 1 < 2.2e-16 ***
## Geno:Trmt 25.4443 1 4.553e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(wid)
## [1] 951
ranova(wid)
## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## LL_Width2_mm ~ Geno + Trmt + LL1_mm + (1 | Rack) + (1 | Fam) +
## Geno:Trmt
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 8 -1881.4 3778.9
## (1 | Rack) 7 -1923.5 3860.9 84.066 1 < 2.2e-16 ***
## (1 | Fam) 7 -1898.2 3810.4 33.562 1 6.904e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### growth rate
gr2<-lmer(GrowthRt ~ Geno + Trmt + Geno*Trmt + (1|Rack) + (1|Fam) + LL1_mm, data=dry, control=lmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
summary(gr2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: GrowthRt ~ Geno + Trmt + Geno * Trmt + (1 | Rack) + (1 | Fam) +
## LL1_mm
## Data: dry
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## REML criterion at convergence: 1877.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5996 -0.4916 0.0870 0.6153 3.5471
##
## Random effects:
## Groups Name Variance Std.Dev.
## Rack (Intercept) 0.07111 0.2667
## Fam (Intercept) 0.02022 0.1422
## Residual 0.37865 0.6153
## Number of obs: 954, groups: Rack, 32; Fam, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.608965 0.105135 36.538313 15.304 < 2e-16 ***
## GenoMM -0.094528 0.106225 10.664104 -0.890 0.39316
## TrmtWet -0.175192 0.057991 922.771759 -3.021 0.00259 **
## LL1_mm -0.025905 0.001766 926.638887 -14.672 < 2e-16 ***
## GenoMM:TrmtWet 0.429988 0.082266 917.953056 5.227 2.14e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) GenoMM TrmtWt LL1_mm
## GenoMM -0.499
## TrmtWet -0.181 0.261
## LL1_mm -0.537 -0.012 -0.151
## GnMM:TrmtWt 0.083 -0.367 -0.711 0.189
Anova(gr2, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: GrowthRt
## Chisq Df Pr(>Chisq)
## (Intercept) 234.2057 1 < 2.2e-16 ***
## Geno 0.7919 1 0.373525
## Trmt 9.1264 1 0.002519 **
## LL1_mm 215.2767 1 < 2.2e-16 ***
## Geno:Trmt 27.3194 1 1.725e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(gr2)
## [1] 954
ranova(gr2)
## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## GrowthRt ~ Geno + Trmt + LL1_mm + (1 | Rack) + (1 | Fam) + Geno:Trmt
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 8 -938.68 1893.4
## (1 | Rack) 7 -984.62 1983.2 91.868 1 < 2.2e-16 ***
## (1 | Fam) 7 -950.49 1915.0 23.615 1 1.177e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### green leaves
glf<-lmer(NumLvs_Grn ~ Geno + Trmt + Geno*Trmt + (1|Rack) + (1|Fam) + LL1_mm, data=dry, control=lmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
summary(glf)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: NumLvs_Grn ~ Geno + Trmt + Geno * Trmt + (1 | Rack) + (1 | Fam) +
## LL1_mm
## Data: dry
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## REML criterion at convergence: 4437.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.6275 -0.5629 -0.0521 0.5536 4.5457
##
## Random effects:
## Groups Name Variance Std.Dev.
## Rack (Intercept) 0.20639 0.4543
## Fam (Intercept) 0.08516 0.2918
## Residual 5.74036 2.3959
## Number of obs: 959, groups: Rack, 32; Fam, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 9.771777 0.300931 61.771991 32.472 < 2e-16 ***
## GenoMM -0.920397 0.285998 15.092490 -3.218 0.00571 **
## TrmtWet -0.485014 0.222697 943.990125 -2.178 0.02966 *
## LL1_mm 0.047072 0.006592 841.281137 7.141 2.01e-12 ***
## GenoMM:TrmtWet 1.516285 0.317653 939.977577 4.773 2.10e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) GenoMM TrmtWt LL1_mm
## GenoMM -0.455
## TrmtWet -0.253 0.376
## LL1_mm -0.692 -0.027 -0.146
## GnMM:TrmtWt 0.120 -0.529 -0.710 0.186
Anova(glf, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: NumLvs_Grn
## Chisq Df Pr(>Chisq)
## (Intercept) 1054.4218 1 < 2.2e-16 ***
## Geno 10.3568 1 0.00129 **
## Trmt 4.7433 1 0.02941 *
## LL1_mm 50.9936 1 9.267e-13 ***
## Geno:Trmt 22.7854 1 1.811e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(glf)
## [1] 959
ranova(glf)
## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## NumLvs_Grn ~ Geno + Trmt + LL1_mm + (1 | Rack) + (1 | Fam) +
## Geno:Trmt
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 8 -2218.8 4453.7
## (1 | Rack) 7 -2224.0 4461.9 10.2554 1 0.001363 **
## (1 | Fam) 7 -2220.8 4455.5 3.8712 1 0.049120 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### LWC
lwc<-lmer(LWC_Ros ~ Geno + Trmt + Geno*Trmt + (1|Rack) + (1|Fam) + LL1_mm, data=dry, control=lmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
summary(lwc)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LWC_Ros ~ Geno + Trmt + Geno * Trmt + (1 | Rack) + (1 | Fam) +
## LL1_mm
## Data: dry
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## REML criterion at convergence: -1275
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.4728 -0.4046 0.0925 0.5552 3.0604
##
## Random effects:
## Groups Name Variance Std.Dev.
## Rack (Intercept) 3.073e-03 5.543e-02
## Fam (Intercept) 1.002e-18 1.001e-09
## Residual 5.077e-03 7.125e-02
## Number of obs: 569, groups: Rack, 32; Fam, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 6.049e-01 1.347e-02 8.356e+01 44.911 <2e-16 ***
## GenoMM -6.411e-03 6.663e-03 5.372e+02 -0.962 0.336
## TrmtWet 1.220e-01 1.246e-02 5.527e+02 9.796 <2e-16 ***
## LL1_mm -2.655e-03 2.478e-04 5.467e+02 -10.714 <2e-16 ***
## GenoMM:TrmtWet 3.015e-02 1.702e-02 5.344e+02 1.772 0.077 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) GenoMM TrmtWt LL1_mm
## GenoMM -0.220
## TrmtWet -0.084 0.278
## LL1_mm -0.591 -0.043 -0.075
## GnMM:TrmtWt 0.027 -0.397 -0.608 0.120
## convergence code: 0
## boundary (singular) fit: see ?isSingular
Anova(lwc, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: LWC_Ros
## Chisq Df Pr(>Chisq)
## (Intercept) 2017.0092 1 < 2e-16 ***
## Geno 0.9257 1 0.33598
## Trmt 95.9684 1 < 2e-16 ***
## LL1_mm 114.7917 1 < 2e-16 ***
## Geno:Trmt 3.1395 1 0.07642 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(lwc)
## [1] 569
ranova(lwc)
## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## LWC_Ros ~ Geno + Trmt + LL1_mm + (1 | Rack) + (1 | Fam) + Geno:Trmt
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 8 637.50 -1259.0
## (1 | Rack) 7 565.18 -1116.3 144.66 1 <2e-16 ***
## (1 | Fam) 7 637.50 -1261.0 0.00 1 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Morphological predictors of water use under drought:
library(car)
# Supplementary Table 7
## Supp Table 7A: size predictors
wu<-lmer(LSM_daily_norack ~ Geno + (1|Rack) + LL_Length2_mm + Geno*LL_Length2_mm + LL1_mm + (1|Fam), data=dry, control=lmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
summary(wu)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LSM_daily_norack ~ Geno + (1 | Rack) + LL_Length2_mm + Geno *
## LL_Length2_mm + LL1_mm + (1 | Fam)
## Data: dry
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## REML criterion at convergence: 467.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7835 -0.6379 -0.0872 0.5889 4.0839
##
## Random effects:
## Groups Name Variance Std.Dev.
## Rack (Intercept) 0.051055 0.2260
## Fam (Intercept) 0.009389 0.0969
## Residual 0.121463 0.3485
## Number of obs: 486, groups: Rack, 32; Fam, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.681253 0.105799 116.568927 53.699 < 2e-16 ***
## GenoMM -0.039813 0.129640 78.073446 -0.307 0.759585
## LL_Length2_mm 0.007380 0.002061 471.062504 3.582 0.000377 ***
## LL1_mm 0.012136 0.001668 458.176040 7.274 1.52e-12 ***
## GenoMM:LL_Length2_mm -0.001726 0.002429 455.324507 -0.711 0.477632
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) GenoMM LL_L2_ LL1_mm
## GenoMM -0.646
## LL_Lngth2_m -0.665 0.537
## LL1_mm -0.091 -0.018 -0.460
## GnMM:LL_L2_ 0.549 -0.845 -0.615 0.001
Anova(wu, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: LSM_daily_norack
## Chisq Df Pr(>Chisq)
## (Intercept) 2883.5450 1 < 2.2e-16 ***
## Geno 0.0943 1 0.7587674
## LL_Length2_mm 12.8290 1 0.0003413 ***
## LL1_mm 52.9137 1 3.485e-13 ***
## Geno:LL_Length2_mm 0.5051 1 0.4772678
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(wu)
## [1] 486
ranova(wu)
## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## LSM_daily_norack ~ Geno + LL_Length2_mm + LL1_mm + (1 | Rack) +
## (1 | Fam) + Geno:LL_Length2_mm
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 8 -233.57 483.15
## (1 | Rack) 7 -274.35 562.71 81.558 1 < 2.2e-16 ***
## (1 | Fam) 7 -240.04 494.08 12.934 1 0.0003227 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Supp Table 7B: leaf number
wu2<-lmer(LSM_daily_norack ~ Geno + (1|Rack) + TotalLvs + Geno*TotalLvs + (1|Fam), data=dry, control=lmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
summary(wu2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LSM_daily_norack ~ Geno + (1 | Rack) + TotalLvs + Geno * TotalLvs +
## (1 | Fam)
## Data: dry
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## REML criterion at convergence: 521.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0183 -0.6387 -0.0592 0.5886 3.2906
##
## Random effects:
## Groups Name Variance Std.Dev.
## Rack (Intercept) 0.05144 0.2268
## Fam (Intercept) 0.01954 0.1398
## Residual 0.14007 0.3743
## Number of obs: 487, groups: Rack, 32; Fam, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.677153 0.135205 106.297637 41.989 < 2e-16 ***
## GenoMM -0.133408 0.198389 115.703617 -0.672 0.503
## TotalLvs 0.048018 0.007232 462.263441 6.640 8.81e-11 ***
## GenoMM:TotalLvs 0.001593 0.011543 457.947492 0.138 0.890
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) GenoMM TtlLvs
## GenoMM -0.594
## TotalLvs -0.816 0.522
## GnMM:TtlLvs 0.480 -0.878 -0.588
Anova(wu2, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: LSM_daily_norack
## Chisq Df Pr(>Chisq)
## (Intercept) 1763.0923 1 < 2.2e-16 ***
## Geno 0.4522 1 0.5013
## TotalLvs 44.0913 1 3.134e-11 ***
## Geno:TotalLvs 0.0190 1 0.8902
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(wu2)
## [1] 487
ranova(wu2)
## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## LSM_daily_norack ~ Geno + TotalLvs + (1 | Rack) + (1 | Fam) +
## Geno:TotalLvs
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 7 -260.68 535.37
## (1 | Rack) 6 -296.46 604.93 71.564 1 < 2.2e-16 ***
## (1 | Fam) 6 -275.65 563.30 29.934 1 4.47e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Supp Table 7C: LMA
### exclude LMA extreme outliers
dry.use<-filter(dry, LMAest<.1)
### fit model
wu3<-lmer(LSM_daily_norack ~ Geno + (1|Rack) + LMAest + Geno*LMAest + (1|Fam), data=dry.use, control=lmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
summary(wu3)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LSM_daily_norack ~ Geno + (1 | Rack) + LMAest + Geno * LMAest +
## (1 | Fam)
## Data: dry.use
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## REML criterion at convergence: 523.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6013 -0.6387 -0.0511 0.5754 3.3613
##
## Random effects:
## Groups Name Variance Std.Dev.
## Rack (Intercept) 0.06547 0.2559
## Fam (Intercept) 0.02452 0.1566
## Residual 0.14977 0.3870
## Number of obs: 477, groups: Rack, 32; Fam, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.94063 0.13275 70.79656 44.752 < 2e-16 ***
## GenoMM 0.05163 0.18322 65.84865 0.282 0.779
## LMAest 21.59883 4.55808 445.24385 4.739 2.9e-06 ***
## GenoMM:LMAest -8.08351 6.74018 438.38550 -1.199 0.231
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) GenoMM LMAest
## GenoMM -0.615
## LMAest -0.754 0.511
## GenMM:LMAst 0.478 -0.817 -0.633
Anova(wu3, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: LSM_daily_norack
## Chisq Df Pr(>Chisq)
## (Intercept) 2002.7414 1 < 2.2e-16 ***
## Geno 0.0794 1 0.7781
## LMAest 22.4542 1 2.152e-06 ***
## Geno:LMAest 1.4383 1 0.2304
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(wu3)
## [1] 477
ranova(wu3)
## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## LSM_daily_norack ~ Geno + LMAest + (1 | Rack) + (1 | Fam) + Geno:LMAest
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 7 -261.95 537.90
## (1 | Rack) 6 -300.86 613.71 77.804 1 < 2.2e-16 ***
## (1 | Fam) 6 -279.55 571.10 35.192 1 2.987e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Supp Table 7D: all predictors together
wu4<-lmer(LSM_daily_norack ~ Geno + (1|Rack) + LL_Length2_mm + TotalLvs + LMAest + LL1_mm + (1|Fam), data=dry, control=lmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun=2e5)))
summary(wu4)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LSM_daily_norack ~ Geno + (1 | Rack) + LL_Length2_mm + TotalLvs +
## LMAest + LL1_mm + (1 | Fam)
## Data: dry
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## REML criterion at convergence: 450.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8906 -0.6486 -0.0963 0.5863 4.3035
##
## Random effects:
## Groups Name Variance Std.Dev.
## Rack (Intercept) 0.050236 0.22413
## Fam (Intercept) 0.009146 0.09563
## Residual 0.120494 0.34712
## Number of obs: 478, groups: Rack, 32; Fam, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.527581 0.121694 174.944551 45.422 < 2e-16 ***
## GenoMM -0.116958 0.068845 7.616420 -1.699 0.130
## LL_Length2_mm 0.007819 0.001787 458.644204 4.375 1.50e-05 ***
## TotalLvs 0.007588 0.007518 458.892479 1.009 0.313
## LMAest 3.494983 1.770551 449.056322 1.974 0.049 *
## LL1_mm 0.010221 0.002008 451.196513 5.091 5.24e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) GenoMM LL_L2_ TtlLvs LMAest
## GenoMM -0.322
## LL_Lngth2_m -0.393 0.039
## TotalLvs -0.551 0.033 -0.165
## LMAest -0.300 0.036 0.350 -0.150
## LL1_mm 0.296 -0.056 -0.467 -0.445 -0.255
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
Anova(wu4, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: LSM_daily_norack
## Chisq Df Pr(>Chisq)
## (Intercept) 2063.1495 1 < 2.2e-16 ***
## Geno 2.8861 1 0.08935 .
## LL_Length2_mm 19.1444 1 1.212e-05 ***
## TotalLvs 1.0186 1 0.31285
## LMAest 3.8965 1 0.04839 *
## LL1_mm 25.9175 1 3.563e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(wu4)
## [1] 478
ranova(wu4)
## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## LSM_daily_norack ~ Geno + LL_Length2_mm + TotalLvs + LMAest +
## LL1_mm + (1 | Rack) + (1 | Fam)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 9 -225.15 468.30
## (1 | Rack) 8 -264.98 545.96 79.661 1 < 2.2e-16 ***
## (1 | Fam) 8 -231.50 478.99 12.696 1 0.0003664 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Genetic correlation between LWC and BC-ratio:
# Supplementary Table 8A
gwas.drt.reg<-lm(LWC_Dry_LSM ~ BCratio + Group, data=rp)
summary(gwas.drt.reg)
##
## Call:
## lm(formula = LWC_Dry_LSM ~ BCratio + Group, data = rp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0100583 -0.0030995 -0.0005596 0.0025897 0.0163209
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.237e-01 1.330e-03 469.139 < 2e-16 ***
## BCratio -2.625e-05 9.181e-06 -2.858 0.00465 **
## GroupCOL 9.843e-04 1.386e-03 0.710 0.47848
## GroupNOR 1.650e-03 1.332e-03 1.238 0.21685
## GroupUTA 1.490e-03 1.314e-03 1.134 0.25815
## GroupWES 2.611e-03 1.366e-03 1.912 0.05716 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.004347 on 231 degrees of freedom
## (325 observations deleted due to missingness)
## Multiple R-squared: 0.04238, Adjusted R-squared: 0.02165
## F-statistic: 2.044 on 5 and 231 DF, p-value: 0.07337
Anova(gwas.drt.reg, type="3")
## Anova Table (Type III tests)
##
## Response: LWC_Dry_LSM
## Sum Sq Df F value Pr(>F)
## (Intercept) 4.1580 1 2.2009e+05 < 2.2e-16 ***
## BCratio 0.0002 1 8.1709e+00 0.004646 **
## Group 0.0001 4 1.2381e+00 0.295475
## Residuals 0.0044 231
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(gwas.drt.reg)
## [1] 237
# Supplementary Table 8B
gwas.drt.reg.2<-lm(LWC_Dry_LSM ~ BCratio + gPCA1+gPCA2+gPCA3+gPCA4+gPCA5, data=rp)
summary(gwas.drt.reg.2)
##
## Call:
## lm(formula = LWC_Dry_LSM ~ BCratio + gPCA1 + gPCA2 + gPCA3 +
## gPCA4 + gPCA5, data = rp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0099812 -0.0029282 -0.0005765 0.0027871 0.0171346
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.250e-01 6.275e-04 996.054 <2e-16 ***
## BCratio -2.180e-05 8.412e-06 -2.591 0.0101 *
## gPCA1 -1.794e-03 8.976e-03 -0.200 0.8418
## gPCA2 -4.061e-03 5.621e-03 -0.722 0.4707
## gPCA3 -8.980e-03 5.768e-03 -1.557 0.1208
## gPCA4 -2.559e-03 6.387e-03 -0.401 0.6890
## gPCA5 -5.023e-03 5.542e-03 -0.906 0.3656
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.00432 on 248 degrees of freedom
## (307 observations deleted due to missingness)
## Multiple R-squared: 0.03447, Adjusted R-squared: 0.01111
## F-statistic: 1.475 on 6 and 248 DF, p-value: 0.187
Anova(gwas.drt.reg.2, type="3")
## Anova Table (Type III tests)
##
## Response: LWC_Dry_LSM
## Sum Sq Df F value Pr(>F)
## (Intercept) 18.5115 1 9.9212e+05 < 2e-16 ***
## BCratio 0.0001 1 6.7133e+00 0.01014 *
## gPCA1 0.0000 1 3.9900e-02 0.84176
## gPCA2 0.0000 1 5.2180e-01 0.47074
## gPCA3 0.0000 1 2.4239e+00 0.12077
## gPCA4 0.0000 1 1.6060e-01 0.68898
## gPCA5 0.0000 1 8.2160e-01 0.36559
## Residuals 0.0046 248
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(gwas.drt.reg.2)
## [1] 255
# Supplementary Table 8C
gwas.drt.reg.3<-lm(LWC_Dry_LSM ~ BCratio + colutaPCA1+colutaPCA2+colutaPCA3+colutaPCA4+colutaPCA5, data=rp)
summary(gwas.drt.reg.3)
##
## Call:
## lm(formula = LWC_Dry_LSM ~ BCratio + colutaPCA1 + colutaPCA2 +
## colutaPCA3 + colutaPCA4 + colutaPCA5, data = rp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.009775 -0.003266 -0.000705 0.002878 0.014638
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.249e-01 8.644e-04 722.962 <2e-16 ***
## BCratio -2.456e-05 1.277e-05 -1.924 0.0572 .
## colutaPCA1 4.516e-03 5.258e-03 0.859 0.3924
## colutaPCA2 1.020e-03 5.113e-03 0.200 0.8422
## colutaPCA3 9.411e-03 5.049e-03 1.864 0.0651 .
## colutaPCA4 3.561e-03 5.147e-03 0.692 0.4906
## colutaPCA5 -5.994e-03 5.016e-03 -1.195 0.2349
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.004719 on 103 degrees of freedom
## (452 observations deleted due to missingness)
## Multiple R-squared: 0.1116, Adjusted R-squared: 0.05985
## F-statistic: 2.156 on 6 and 103 DF, p-value: 0.0532
Anova(gwas.drt.reg.3, type="3")
## Anova Table (Type III tests)
##
## Response: LWC_Dry_LSM
## Sum Sq Df F value Pr(>F)
## (Intercept) 11.6405 1 5.2267e+05 < 2e-16 ***
## BCratio 0.0001 1 3.7000e+00 0.05717 .
## colutaPCA1 0.0000 1 7.3760e-01 0.39242
## colutaPCA2 0.0000 1 3.9800e-02 0.84222
## colutaPCA3 0.0001 1 3.4751e+00 0.06515 .
## colutaPCA4 0.0000 1 4.7860e-01 0.49061
## colutaPCA5 0.0000 1 1.4276e+00 0.23490
## Residuals 0.0023 103
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(gwas.drt.reg.3)
## [1] 110Effects of climate PCs on BC-ratio:
library(car)
# Supplementary Table 9A
clim.use<-filter(rp, !is.na(Group))
clim.pca<-lm(BCratio ~ Group + Latitude + Longitude + Elev_m + Prin1 + Prin2 + Prin3 +Prin4 + Prin5, data=clim.use)
summary(clim.pca)
##
## Call:
## lm(formula = BCratio ~ Group + Latitude + Longitude + Elev_m +
## Prin1 + Prin2 + Prin3 + Prin4 + Prin5, data = clim.use)
##
## Residuals:
## Min 1Q Median 3Q Max
## -69.774 -23.979 9.857 23.473 54.352
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 767.37533 317.83140 2.414 0.016426 *
## GroupCOL -22.88750 13.80455 -1.658 0.098485 .
## GroupNOR 22.43580 11.47151 1.956 0.051522 .
## GroupUTA 2.69566 9.97427 0.270 0.787165
## GroupWES 34.34359 12.21289 2.812 0.005283 **
## Latitude -3.20089 2.81812 -1.136 0.257037
## Longitude 3.88218 2.24420 1.730 0.084796 .
## Elev_m -0.05521 0.02316 -2.384 0.017829 *
## Prin1 6.04188 1.55452 3.887 0.000128 ***
## Prin2 0.43460 1.67307 0.260 0.795247
## Prin3 -2.20075 2.02478 -1.087 0.278047
## Prin4 -4.37423 2.46018 -1.778 0.076528 .
## Prin5 -1.38109 2.58700 -0.534 0.593879
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30.94 on 270 degrees of freedom
## (201 observations deleted due to missingness)
## Multiple R-squared: 0.249, Adjusted R-squared: 0.2156
## F-statistic: 7.459 on 12 and 270 DF, p-value: 7.317e-12
Anova(clim.pca, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 5579 1 5.8294 0.016426 *
## Group 11399 4 2.9778 0.019735 *
## Latitude 1235 1 1.2901 0.257037
## Longitude 2864 1 2.9925 0.084796 .
## Elev_m 5438 1 5.6823 0.017829 *
## Prin1 14456 1 15.1060 0.000128 ***
## Prin2 65 1 0.0675 0.795247
## Prin3 1131 1 1.1814 0.278047
## Prin4 3025 1 3.1613 0.076528 .
## Prin5 273 1 0.2850 0.593879
## Residuals 258389 270
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(clim.pca)
## [1] 283
# Supplementary Table 9B
clim.pca.2<-lm(BCratio ~ gPCA1+gPCA2+gPCA3+gPCA4+gPCA5 + Latitude + Longitude + Elev_m + Prin1 + Prin2 + Prin3 +Prin4 + Prin5, data=rp)
summary(clim.pca.2)
##
## Call:
## lm(formula = BCratio ~ gPCA1 + gPCA2 + gPCA3 + gPCA4 + gPCA5 +
## Latitude + Longitude + Elev_m + Prin1 + Prin2 + Prin3 + Prin4 +
## Prin5, data = rp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -79.67 -22.80 9.02 22.39 58.01
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 280.22771 216.93295 1.292 0.197476
## gPCA1 -68.89008 72.28922 -0.953 0.341403
## gPCA2 -9.50791 56.25149 -0.169 0.865896
## gPCA3 -36.25943 55.71415 -0.651 0.515689
## gPCA4 -75.94926 46.63565 -1.629 0.104502
## gPCA5 -32.33279 37.33405 -0.866 0.387192
## Latitude -0.02286 2.35447 -0.010 0.992261
## Longitude 0.75566 1.66746 0.453 0.650762
## Elev_m -0.05055 0.02170 -2.330 0.020523 *
## Prin1 5.41360 1.47403 3.673 0.000286 ***
## Prin2 1.96863 1.70342 1.156 0.248768
## Prin3 -3.89549 1.77502 -2.195 0.028992 *
## Prin4 -3.25136 2.39034 -1.360 0.174831
## Prin5 -1.87265 2.72882 -0.686 0.493110
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 31.32 on 287 degrees of freedom
## (261 observations deleted due to missingness)
## Multiple R-squared: 0.2199, Adjusted R-squared: 0.1846
## F-statistic: 6.223 on 13 and 287 DF, p-value: 2.65e-10
Anova(clim.pca.2, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 1636 1 1.6687 0.1974757
## gPCA1 891 1 0.9082 0.3414028
## gPCA2 28 1 0.0286 0.8658960
## gPCA3 415 1 0.4236 0.5156886
## gPCA4 2601 1 2.6522 0.1045020
## gPCA5 736 1 0.7500 0.3871915
## Latitude 0 1 0.0001 0.9922615
## Longitude 201 1 0.2054 0.6507621
## Elev_m 5322 1 5.4269 0.0205230 *
## Prin1 13227 1 13.4883 0.0002863 ***
## Prin2 1310 1 1.3356 0.2487680
## Prin3 4723 1 4.8164 0.0289919 *
## Prin4 1814 1 1.8502 0.1748309
## Prin5 462 1 0.4709 0.4931099
## Residuals 281444 287
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(clim.pca.2)
## [1] 301
# Supplementary Table 9C
clim.pca.3<-lm(BCratio ~ colutaPCA1+colutaPCA2+colutaPCA3+colutaPCA4+colutaPCA5 + Latitude + Longitude + Elev_m + Prin1 + Prin2 + Prin3 +Prin4 + Prin5, data=rp)
summary(clim.pca.3)
##
## Call:
## lm(formula = BCratio ~ colutaPCA1 + colutaPCA2 + colutaPCA3 +
## colutaPCA4 + colutaPCA5 + Latitude + Longitude + Elev_m +
## Prin1 + Prin2 + Prin3 + Prin4 + Prin5, data = rp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -67.414 -27.121 4.959 24.451 63.798
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1056.9635 467.1186 2.263 0.02548 *
## colutaPCA1 111.6806 80.2442 1.392 0.16661
## colutaPCA2 84.9290 96.2540 0.882 0.37939
## colutaPCA3 -75.0512 45.2579 -1.658 0.09991 .
## colutaPCA4 -42.8002 46.0758 -0.929 0.35483
## colutaPCA5 66.0869 47.0720 1.404 0.16296
## Latitude -6.8074 4.9414 -1.378 0.17093
## Longitude 4.0521 3.9138 1.035 0.30262
## Elev_m -0.1031 0.0431 -2.392 0.01835 *
## Prin1 10.0592 3.7850 2.658 0.00896 **
## Prin2 -0.2182 3.2787 -0.067 0.94704
## Prin3 -6.9895 4.1252 -1.694 0.09284 .
## Prin4 -10.1140 5.1457 -1.966 0.05170 .
## Prin5 -4.4414 5.6652 -0.784 0.43462
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 33.74 on 118 degrees of freedom
## (430 observations deleted due to missingness)
## Multiple R-squared: 0.3266, Adjusted R-squared: 0.2525
## F-statistic: 4.403 on 13 and 118 DF, p-value: 4.879e-06
Anova(clim.pca.3, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 5830 1 5.1200 0.025481 *
## colutaPCA1 2205 1 1.9370 0.166613
## colutaPCA2 886 1 0.7785 0.379386
## colutaPCA3 3131 1 2.7500 0.099913 .
## colutaPCA4 982 1 0.8629 0.354833
## colutaPCA5 2244 1 1.9711 0.162960
## Latitude 2161 1 1.8978 0.170930
## Longitude 1221 1 1.0720 0.302620
## Elev_m 6513 1 5.7204 0.018349 *
## Prin1 8042 1 7.0631 0.008959 **
## Prin2 5 1 0.0044 0.947040
## Prin3 3269 1 2.8708 0.092840 .
## Prin4 4399 1 3.8634 0.051700 .
## Prin5 700 1 0.6146 0.434617
## Residuals 134357 118
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(clim.pca.3)
## [1] 132Axis loading in Supplementary Table 10 analyzed in JMP.
Effects of raw climate predicors on BC-ratio:
library(car)
# Supplementary Table 11
## max temp. in warmest month
### model A
bc5<-lm(BCratio ~ bio5_12 + Group, data=clim.use)
summary(bc5)
##
## Call:
## lm(formula = BCratio ~ bio5_12 + Group, data = clim.use)
##
## Residuals:
## Min 1Q Median 3Q Max
## -64.632 -25.666 6.684 24.537 49.753
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 86.48899 22.77335 3.798 0.000179 ***
## bio5_12 -0.11655 0.09834 -1.185 0.236964
## GroupCOL -14.60503 9.21810 -1.584 0.114247
## GroupNOR 5.50130 8.86755 0.620 0.535514
## GroupUTA -0.02333 8.81616 -0.003 0.997891
## GroupWES 30.75290 9.20246 3.342 0.000947 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 32.1 on 277 degrees of freedom
## (201 observations deleted due to missingness)
## Multiple R-squared: 0.1704, Adjusted R-squared: 0.1554
## F-statistic: 11.38 on 5 and 277 DF, p-value: 5.361e-10
Anova(bc5, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 14862 1 14.4234 0.0001794 ***
## bio5_12 1447 1 1.4046 0.2369645
## Group 58590 4 14.2150 1.447e-10 ***
## Residuals 285429 277
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc5)
## [1] 283
### model B
bc5b<-lm(BCratio ~ bio5_12 + gPCA1+gPCA2+gPCA3+gPCA4+gPCA5 , data=rp)
summary(bc5b)
##
## Call:
## lm(formula = BCratio ~ bio5_12 + gPCA1 + gPCA2 + gPCA3 + gPCA4 +
## gPCA5, data = rp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -82.12 -13.50 13.75 23.36 52.56
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 71.03179 22.80899 3.114 0.002025 **
## bio5_12 -0.03121 0.10072 -0.310 0.756861
## gPCA1 -167.34469 61.77015 -2.709 0.007139 **
## gPCA2 -78.02524 39.83432 -1.959 0.051081 .
## gPCA3 -141.42467 41.65034 -3.396 0.000779 ***
## gPCA4 -9.20119 46.11283 -0.200 0.841980
## gPCA5 -50.53263 39.64434 -1.275 0.203433
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 33.71 on 296 degrees of freedom
## (259 observations deleted due to missingness)
## Multiple R-squared: 0.06937, Adjusted R-squared: 0.0505
## F-statistic: 3.677 on 6 and 296 DF, p-value: 0.00154
Anova(bc5b, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 11022 1 9.6983 0.0020250 **
## bio5_12 109 1 0.0960 0.7568614
## gPCA1 8341 1 7.3395 0.0071385 **
## gPCA2 4360 1 3.8367 0.0510815 .
## gPCA3 13103 1 11.5296 0.0007786 ***
## gPCA4 45 1 0.0398 0.8419802
## gPCA5 1847 1 1.6247 0.2034331
## Residuals 336407 296
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc5b)
## [1] 303
### model C
bc5c<-lm(BCratio ~ bio5_12 + colutaPCA1+colutaPCA2+colutaPCA3+colutaPCA4+colutaPCA5 , data=rp)
summary(bc5c)
##
## Call:
## lm(formula = BCratio ~ bio5_12 + colutaPCA1 + colutaPCA2 + colutaPCA3 +
## colutaPCA4 + colutaPCA5, data = rp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -69.47 -36.99 11.97 30.74 59.31
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 101.3733 45.2398 2.241 0.0268 *
## bio5_12 -0.2083 0.2015 -1.034 0.3033
## colutaPCA1 91.4868 37.1469 2.463 0.0151 *
## colutaPCA2 50.0404 36.7856 1.360 0.1762
## colutaPCA3 -103.6384 43.0082 -2.410 0.0174 *
## colutaPCA4 -61.0502 37.4362 -1.631 0.1055
## colutaPCA5 -6.9976 36.7991 -0.190 0.8495
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 36.79 on 125 degrees of freedom
## (430 observations deleted due to missingness)
## Multiple R-squared: 0.1523, Adjusted R-squared: 0.1116
## F-statistic: 3.743 on 6 and 125 DF, p-value: 0.001853
Anova(bc5c, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 6795 1 5.0212 0.02680 *
## bio5_12 1446 1 1.0683 0.30332
## colutaPCA1 8208 1 6.0656 0.01515 *
## colutaPCA2 2504 1 1.8505 0.17618
## colutaPCA3 7858 1 5.8068 0.01742 *
## colutaPCA4 3599 1 2.6594 0.10545
## colutaPCA5 49 1 0.0362 0.84950
## Residuals 169147 125
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc5c)
## [1] 132
## annual precipitation
### model A
bc12<-lm(BCratio ~ bio12_12 + Group, data=clim.use)
summary(bc12)
##
## Call:
## lm(formula = BCratio ~ bio12_12 + Group, data = clim.use)
##
## Residuals:
## Min 1Q Median 3Q Max
## -68.14 -24.23 9.05 24.26 55.12
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.98641 13.54149 1.624 0.105592
## bio12_12 0.07421 0.02084 3.561 0.000435 ***
## GroupCOL -14.71369 9.01905 -1.631 0.103942
## GroupNOR 10.83660 8.83265 1.227 0.220910
## GroupUTA -0.64720 8.60981 -0.075 0.940134
## GroupWES 32.58201 8.94739 3.642 0.000323 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 31.47 on 277 degrees of freedom
## (201 observations deleted due to missingness)
## Multiple R-squared: 0.2027, Adjusted R-squared: 0.1883
## F-statistic: 14.08 on 5 and 277 DF, p-value: 2.821e-12
Anova(bc12, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 2611 1 2.6362 0.1055916
## bio12_12 12555 1 12.6778 0.0004354 ***
## Group 65541 4 16.5453 3.598e-12 ***
## Residuals 274322 277
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc12)
## [1] 283
### model B
bc12b<-lm(BCratio ~ bio12_12 + gPCA1+gPCA2+gPCA3+gPCA4+gPCA5, data=rp)
summary(bc12b)
##
## Call:
## lm(formula = BCratio ~ bio12_12 + gPCA1 + gPCA2 + gPCA3 + gPCA4 +
## gPCA5, data = rp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -76.86 -17.50 12.58 23.50 51.82
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30.91277 11.02639 2.804 0.00539 **
## bio12_12 0.06638 0.02176 3.051 0.00249 **
## gPCA1 -181.64707 60.89438 -2.983 0.00309 **
## gPCA2 -124.36306 42.03333 -2.959 0.00334 **
## gPCA3 -164.22040 41.42907 -3.964 9.25e-05 ***
## gPCA4 -30.01382 45.89501 -0.654 0.51364
## gPCA5 -48.59144 38.96171 -1.247 0.21333
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 33.2 on 296 degrees of freedom
## (259 observations deleted due to missingness)
## Multiple R-squared: 0.09744, Adjusted R-squared: 0.07915
## F-statistic: 5.326 on 6 and 296 DF, p-value: 3.09e-05
Anova(bc12b, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 8663 1 7.8598 0.005389 **
## bio12_12 10257 1 9.3057 0.002491 **
## gPCA1 9808 1 8.8982 0.003092 **
## gPCA2 9649 1 8.7538 0.003339 **
## gPCA3 17319 1 15.7124 9.25e-05 ***
## gPCA4 471 1 0.4277 0.513641
## gPCA5 1714 1 1.5554 0.213325
## Residuals 326259 296
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc12b)
## [1] 303
### model C
bc12c<-lm(BCratio ~ bio12_12 + colutaPCA1+colutaPCA2+colutaPCA3+colutaPCA4+colutaPCA5, data=rp)
summary(bc12c)
##
## Call:
## lm(formula = BCratio ~ bio12_12 + colutaPCA1 + colutaPCA2 + colutaPCA3 +
## colutaPCA4 + colutaPCA5, data = rp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -67.899 -33.183 8.784 31.353 58.354
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -18.16546 25.23954 -0.720 0.47304
## bio12_12 0.13934 0.04788 2.910 0.00427 **
## colutaPCA1 90.27632 35.77871 2.523 0.01288 *
## colutaPCA2 24.14510 36.86889 0.655 0.51374
## colutaPCA3 -59.25154 42.59940 -1.391 0.16673
## colutaPCA4 -19.98069 37.59934 -0.531 0.59608
## colutaPCA5 -26.24634 36.42508 -0.721 0.47253
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 35.75 on 125 degrees of freedom
## (430 observations deleted due to missingness)
## Multiple R-squared: 0.1993, Adjusted R-squared: 0.1609
## F-statistic: 5.186 on 6 and 125 DF, p-value: 8.547e-05
Anova(bc12c, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 662 1 0.5180 0.473040
## bio12_12 10826 1 8.4704 0.004275 **
## colutaPCA1 8137 1 6.3665 0.012884 *
## colutaPCA2 548 1 0.4289 0.513741
## colutaPCA3 2473 1 1.9346 0.166727
## colutaPCA4 361 1 0.2824 0.596077
## colutaPCA5 664 1 0.5192 0.472528
## Residuals 159766 125
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc12c)
## [1] 132
## precip. in wettest month
### model A
bc13<-lm(BCratio ~ bio13_12 + Group, data=clim.use)
summary(bc13)
##
## Call:
## lm(formula = BCratio ~ bio13_12 + Group, data = clim.use)
##
## Residuals:
## Min 1Q Median 3Q Max
## -64.21 -29.50 9.78 23.94 49.34
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.4261 12.1107 2.595 0.00997 **
## bio13_12 0.5100 0.1571 3.246 0.00131 **
## GroupCOL -14.3482 9.0559 -1.584 0.11424
## GroupNOR 7.5600 8.7325 0.866 0.38739
## GroupUTA -0.2106 8.6444 -0.024 0.98058
## GroupWES 26.1274 8.9596 2.916 0.00383 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 31.59 on 277 degrees of freedom
## (201 observations deleted due to missingness)
## Multiple R-squared: 0.1967, Adjusted R-squared: 0.1822
## F-statistic: 13.57 on 5 and 277 DF, p-value: 7.555e-12
Anova(bc13, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 6718 1 6.7336 0.009966 **
## bio13_12 10511 1 10.5355 0.001315 **
## Group 46721 4 11.7070 8.467e-09 ***
## Residuals 276365 277
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc13)
## [1] 283
### model B
bc13b<-lm(BCratio ~ bio13_12 + gPCA1+gPCA2+gPCA3+gPCA4+gPCA5, data=rp)
summary(bc13b)
##
## Call:
## lm(formula = BCratio ~ bio13_12 + gPCA1 + gPCA2 + gPCA3 + gPCA4 +
## gPCA5, data = rp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -77.96 -21.38 10.20 22.89 55.10
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.5987 9.4592 2.072 0.03914 *
## bio13_12 0.7694 0.1604 4.797 2.56e-06 ***
## gPCA1 -176.3879 59.3996 -2.970 0.00323 **
## gPCA2 -121.1120 39.2882 -3.083 0.00224 **
## gPCA3 -122.3305 39.9222 -3.064 0.00238 **
## gPCA4 -59.2207 45.6079 -1.298 0.19513
## gPCA5 -51.7694 38.1163 -1.358 0.17544
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 32.48 on 296 degrees of freedom
## (259 observations deleted due to missingness)
## Multiple R-squared: 0.1362, Adjusted R-squared: 0.1187
## F-statistic: 7.779 on 6 and 296 DF, p-value: 8.75e-08
Anova(bc13b, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 4528 1 4.2929 0.039138 *
## bio13_12 24271 1 23.0087 2.56e-06 ***
## gPCA1 9302 1 8.8180 0.003227 **
## gPCA2 10024 1 9.5028 0.002245 **
## gPCA3 9905 1 9.3895 0.002383 **
## gPCA4 1779 1 1.6860 0.195135
## gPCA5 1946 1 1.8447 0.175436
## Residuals 312245 296
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc13b)
## [1] 303
### model C
bc13c<-lm(BCratio ~ bio13_12 + colutaPCA1+colutaPCA2+colutaPCA3+colutaPCA4+colutaPCA5, data=rp)
summary(bc13c)
##
## Call:
## lm(formula = BCratio ~ bio13_12 + colutaPCA1 + colutaPCA2 + colutaPCA3 +
## colutaPCA4 + colutaPCA5, data = rp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -68.963 -30.347 5.277 30.042 59.500
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -18.2133 19.7722 -0.921 0.358742
## bio13_12 1.2838 0.3438 3.734 0.000285 ***
## colutaPCA1 71.5247 35.2574 2.029 0.044619 *
## colutaPCA2 60.5407 35.1439 1.723 0.087425 .
## colutaPCA3 -63.1587 38.9497 -1.622 0.107421
## colutaPCA4 -14.7675 36.5713 -0.404 0.687049
## colutaPCA5 -26.4076 35.4656 -0.745 0.457911
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 35.04 on 125 degrees of freedom
## (430 observations deleted due to missingness)
## Multiple R-squared: 0.2308, Adjusted R-squared: 0.1939
## F-statistic: 6.252 on 6 and 125 DF, p-value: 9.11e-06
Anova(bc13c, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 1042 1 0.8485 0.3587423
## bio13_12 17118 1 13.9423 0.0002851 ***
## colutaPCA1 5053 1 4.1154 0.0446191 *
## colutaPCA2 3644 1 2.9675 0.0874249 .
## colutaPCA3 3228 1 2.6294 0.1074207
## colutaPCA4 200 1 0.1631 0.6870488
## colutaPCA5 681 1 0.5544 0.4579112
## Residuals 153474 125
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc13c)
## [1] 132
## precip. in driest month
### model A
bc14<-lm(BCratio ~ bio14_12 + Group, data=clim.use)
summary(bc14)
##
## Call:
## lm(formula = BCratio ~ bio14_12 + Group, data = clim.use)
##
## Residuals:
## Min 1Q Median 3Q Max
## -67.044 -25.250 7.081 24.648 56.509
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 38.4368 13.6267 2.821 0.005138 **
## bio14_12 0.6967 0.3375 2.064 0.039938 *
## GroupCOL -12.9648 9.2224 -1.406 0.160906
## GroupNOR 9.5391 9.1198 1.046 0.296485
## GroupUTA 1.3543 8.8079 0.154 0.877915
## GroupWES 35.7158 9.6086 3.717 0.000244 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 31.94 on 277 degrees of freedom
## (201 observations deleted due to missingness)
## Multiple R-squared: 0.1788, Adjusted R-squared: 0.164
## F-statistic: 12.06 on 5 and 277 DF, p-value: 1.398e-10
Anova(bc14, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 8115 1 7.9564 0.005138 **
## bio14_12 4346 1 4.2606 0.039938 *
## Group 60517 4 14.8330 5.388e-11 ***
## Residuals 282531 277
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc14)
## [1] 283
### model B
bc14b<-lm(BCratio ~ bio14_12 + gPCA1+gPCA2+gPCA3+gPCA4+gPCA5, data=rp)
summary(bc14b)
##
## Call:
## lm(formula = BCratio ~ bio14_12 + gPCA1 + gPCA2 + gPCA3 + gPCA4 +
## gPCA5, data = rp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -82.06 -13.58 13.93 23.23 52.08
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 62.53316 9.46275 6.608 1.81e-10 ***
## bio14_12 0.05315 0.33684 0.158 0.8747
## gPCA1 -166.87764 61.88091 -2.697 0.0074 **
## gPCA2 -78.61588 41.28409 -1.904 0.0578 .
## gPCA3 -142.25518 44.38452 -3.205 0.0015 **
## gPCA4 -7.87691 46.18717 -0.171 0.8647
## gPCA5 -49.29633 39.66441 -1.243 0.2149
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 33.72 on 296 degrees of freedom
## (259 observations deleted due to missingness)
## Multiple R-squared: 0.06914, Adjusted R-squared: 0.05028
## F-statistic: 3.664 on 6 and 296 DF, p-value: 0.001587
Anova(bc14b, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 49644 1 43.6703 1.806e-10 ***
## bio14_12 28 1 0.0249 0.874729
## gPCA1 8267 1 7.2725 0.007403 **
## gPCA2 4122 1 3.6262 0.057845 .
## gPCA3 11678 1 10.2724 0.001498 **
## gPCA4 33 1 0.0291 0.864700
## gPCA5 1756 1 1.5446 0.214912
## Residuals 336488 296
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc14b)
## [1] 303
### model C
bc14c<-lm(BCratio ~ bio14_12 + colutaPCA1+colutaPCA2+colutaPCA3+colutaPCA4+colutaPCA5, data=rp)
summary(bc14c)
##
## Call:
## lm(formula = BCratio ~ bio14_12 + colutaPCA1 + colutaPCA2 + colutaPCA3 +
## colutaPCA4 + colutaPCA5, data = rp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -68.62 -37.15 11.88 29.97 63.92
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.2221 25.1162 1.004 0.31722
## bio14_12 1.0042 0.8477 1.185 0.23845
## colutaPCA1 101.4081 38.9316 2.605 0.01031 *
## colutaPCA2 10.5721 49.7907 0.212 0.83219
## colutaPCA3 -107.3120 40.2074 -2.669 0.00862 **
## colutaPCA4 -57.9553 36.8986 -1.571 0.11879
## colutaPCA5 -5.1293 36.7422 -0.140 0.88920
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 36.74 on 125 degrees of freedom
## (430 observations deleted due to missingness)
## Multiple R-squared: 0.1545, Adjusted R-squared: 0.114
## F-statistic: 3.808 on 6 and 125 DF, p-value: 0.001611
Anova(bc14c, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 1361 1 1.0084 0.317215
## bio14_12 1894 1 1.4031 0.238445
## colutaPCA1 9157 1 6.7849 0.010308 *
## colutaPCA2 61 1 0.0451 0.832195
## colutaPCA3 9614 1 7.1233 0.008618 **
## colutaPCA4 3329 1 2.4670 0.118789
## colutaPCA5 26 1 0.0195 0.889200
## Residuals 168699 125
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc14c)
## [1] 132
## precip. in wettest quarter
### model A
bc16<-lm(BCratio ~ bio16_12 + Group, data=clim.use)
summary(bc16)
##
## Call:
## lm(formula = BCratio ~ bio16_12 + Group, data = clim.use)
##
## Residuals:
## Min 1Q Median 3Q Max
## -64.664 -29.912 9.722 23.467 48.310
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.00080 11.68921 2.738 0.00659 **
## bio16_12 0.18621 0.05498 3.387 0.00081 ***
## GroupCOL -15.24902 9.03630 -1.688 0.09263 .
## GroupNOR 7.58900 8.71502 0.871 0.38462
## GroupUTA -0.47207 8.62853 -0.055 0.95641
## GroupWES 26.45400 8.93369 2.961 0.00333 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 31.54 on 277 degrees of freedom
## (201 observations deleted due to missingness)
## Multiple R-squared: 0.1993, Adjusted R-squared: 0.1849
## F-statistic: 13.79 on 5 and 277 DF, p-value: 4.912e-12
Anova(bc16, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 7453 1 7.4946 0.0065889 **
## bio16_12 11407 1 11.4703 0.0008096 ***
## Group 50361 4 12.6602 1.783e-09 ***
## Residuals 275470 277
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc16)
## [1] 283
### model B
bc16b<-lm(BCratio ~ bio16_12 + gPCA1+gPCA2+gPCA3+gPCA4+gPCA5, data=rp)
summary(bc16b)
##
## Call:
## lm(formula = BCratio ~ bio16_12 + gPCA1 + gPCA2 + gPCA3 + gPCA4 +
## gPCA5, data = rp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -77.55 -21.46 10.57 23.10 54.71
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.7671 9.0143 2.415 0.016353 *
## bio16_12 0.2702 0.0563 4.798 2.54e-06 ***
## gPCA1 -181.1240 59.4423 -3.047 0.002519 **
## gPCA2 -120.6882 39.2658 -3.074 0.002312 **
## gPCA3 -143.0147 39.7630 -3.597 0.000378 ***
## gPCA4 -62.7583 45.7826 -1.371 0.171479
## gPCA5 -52.5443 38.1174 -1.378 0.169095
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 32.48 on 296 degrees of freedom
## (259 observations deleted due to missingness)
## Multiple R-squared: 0.1363, Adjusted R-squared: 0.1187
## F-statistic: 7.782 on 6 and 296 DF, p-value: 8.689e-08
Anova(bc16b, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 6151 1 5.8309 0.0163529 *
## bio16_12 24288 1 23.0253 2.539e-06 ***
## gPCA1 9794 1 9.2846 0.0025193 **
## gPCA2 9965 1 9.4471 0.0023118 **
## gPCA3 13645 1 12.9361 0.0003776 ***
## gPCA4 1982 1 1.8791 0.1714793
## gPCA5 2004 1 1.9002 0.1690949
## Residuals 312229 296
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc16b)
## [1] 303
### model C
bc16c<-lm(BCratio ~ bio16_12 + colutaPCA1+colutaPCA2+colutaPCA3+colutaPCA4+colutaPCA5, data=rp)
summary(bc16c)
##
## Call:
## lm(formula = BCratio ~ bio16_12 + colutaPCA1 + colutaPCA2 + colutaPCA3 +
## colutaPCA4 + colutaPCA5, data = rp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -69.695 -30.620 5.397 29.917 60.002
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -20.1810 19.3294 -1.044 0.298475
## bio16_12 0.4825 0.1229 3.924 0.000143 ***
## colutaPCA1 89.4113 34.8671 2.564 0.011520 *
## colutaPCA2 33.6537 35.1160 0.958 0.339733
## colutaPCA3 -68.4244 37.8858 -1.806 0.073313 .
## colutaPCA4 -3.0381 37.1866 -0.082 0.935016
## colutaPCA5 -22.5663 35.1136 -0.643 0.521618
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 34.86 on 125 degrees of freedom
## (430 observations deleted due to missingness)
## Multiple R-squared: 0.2388, Adjusted R-squared: 0.2023
## F-statistic: 6.536 on 6 and 125 DF, p-value: 5.057e-06
Anova(bc16c, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 1324 1 1.0901 0.2984745
## bio16_12 18711 1 15.3993 0.0001428 ***
## colutaPCA1 7990 1 6.5759 0.0115202 *
## colutaPCA2 1116 1 0.9184 0.3397328
## colutaPCA3 3963 1 3.2619 0.0733133 .
## colutaPCA4 8 1 0.0067 0.9350160
## colutaPCA5 502 1 0.4130 0.5216177
## Residuals 151881 125
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc16c)
## [1] 132
## precip. in driest quarter
### model A
bc17<-lm(BCratio ~ bio17_12 + Group, data=clim.use)
summary(bc17)
##
## Call:
## lm(formula = BCratio ~ bio17_12 + Group, data = clim.use)
##
## Residuals:
## Min 1Q Median 3Q Max
## -67.214 -23.017 7.181 24.505 60.344
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.0806 14.2954 2.384 0.0178 *
## bio17_12 0.2483 0.1085 2.288 0.0229 *
## GroupCOL -13.4026 9.1740 -1.461 0.1452
## GroupNOR 9.6519 9.0543 1.066 0.2874
## GroupUTA 0.5041 8.7451 0.058 0.9541
## GroupWES 35.8184 9.5011 3.770 0.0002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 31.88 on 277 degrees of freedom
## (201 observations deleted due to missingness)
## Multiple R-squared: 0.1816, Adjusted R-squared: 0.1669
## F-statistic: 12.3 on 5 and 277 DF, p-value: 8.863e-11
Anova(bc17, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 5777 1 5.6836 0.0178 *
## bio17_12 5320 1 5.2339 0.0229 *
## Group 61719 4 15.1799 3.103e-11 ***
## Residuals 281557 277
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc17)
## [1] 283
### model B
bc17b<-lm(BCratio ~ bio17_12 + gPCA1+gPCA2+gPCA3+gPCA4+gPCA5, data=rp)
summary(bc17b)
##
## Call:
## lm(formula = BCratio ~ bio17_12 + gPCA1 + gPCA2 + gPCA3 + gPCA4 +
## gPCA5, data = rp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -82.26 -13.52 13.91 23.16 52.18
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.314e+01 1.031e+01 6.127 2.85e-09 ***
## bio17_12 8.892e-03 1.059e-01 0.084 0.93315
## gPCA1 -1.664e+02 6.180e+01 -2.692 0.00750 **
## gPCA2 -7.790e+01 4.177e+01 -1.865 0.06316 .
## gPCA3 -1.407e+02 4.308e+01 -3.267 0.00122 **
## gPCA4 -8.145e+00 4.618e+01 -0.176 0.86013
## gPCA5 -4.960e+01 3.960e+01 -1.252 0.21144
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 33.72 on 296 degrees of freedom
## (259 observations deleted due to missingness)
## Multiple R-squared: 0.06909, Adjusted R-squared: 0.05022
## F-statistic: 3.661 on 6 and 296 DF, p-value: 0.001598
Anova(bc17b, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 42678 1 37.5408 2.848e-09 ***
## bio17_12 8 1 0.0070 0.933146
## gPCA1 8240 1 7.2477 0.007503 **
## gPCA2 3955 1 3.4786 0.063157 .
## gPCA3 12131 1 10.6703 0.001217 **
## gPCA4 35 1 0.0311 0.860128
## gPCA5 1783 1 1.5683 0.211441
## Residuals 336508 296
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc17b)
## [1] 303
### model C
bc17c<-lm(BCratio ~ bio17_12 + colutaPCA1+colutaPCA2+colutaPCA3+colutaPCA4+colutaPCA5, data=rp)
summary(bc17c)
##
## Call:
## lm(formula = BCratio ~ bio17_12 + colutaPCA1 + colutaPCA2 + colutaPCA3 +
## colutaPCA4 + colutaPCA5, data = rp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -69.30 -37.29 12.32 31.80 58.77
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.7665 25.9872 1.299 0.1962
## bio17_12 0.2040 0.2509 0.813 0.4178
## colutaPCA1 89.8805 37.1306 2.421 0.0169 *
## colutaPCA2 36.4796 40.6205 0.898 0.3709
## colutaPCA3 -106.9157 44.1366 -2.422 0.0169 *
## colutaPCA4 -51.0097 37.0126 -1.378 0.1706
## colutaPCA5 -8.2230 36.9518 -0.223 0.8243
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 36.85 on 125 degrees of freedom
## (430 observations deleted due to missingness)
## Multiple R-squared: 0.1495, Adjusted R-squared: 0.1087
## F-statistic: 3.663 on 6 and 125 DF, p-value: 0.002195
Anova(bc17c, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 2292 1 1.6883 0.19622
## bio17_12 897 1 0.6609 0.41779
## colutaPCA1 7955 1 5.8596 0.01693 *
## colutaPCA2 1095 1 0.8065 0.37088
## colutaPCA3 7966 1 5.8679 0.01685 *
## colutaPCA4 2578 1 1.8993 0.17061
## colutaPCA5 67 1 0.0495 0.82426
## Residuals 169695 125
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc17c)
## [1] 132
## precip. in coldest quarter
### model A
bc19<-lm(BCratio ~ bio19_12 + Group, data=clim.use)
summary(bc19)
##
## Call:
## lm(formula = BCratio ~ bio19_12 + Group, data = clim.use)
##
## Residuals:
## Min 1Q Median 3Q Max
## -65.060 -26.411 8.437 24.632 54.141
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 49.74063 9.95284 4.998 1.03e-06 ***
## bio19_12 0.08309 0.04290 1.937 0.05380 .
## GroupCOL -14.30358 9.17478 -1.559 0.12014
## GroupNOR 6.31626 8.83765 0.715 0.47540
## GroupUTA -1.21043 8.74631 -0.138 0.89003
## GroupWES 29.01108 9.02699 3.214 0.00146 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 31.97 on 277 degrees of freedom
## (201 observations deleted due to missingness)
## Multiple R-squared: 0.1773, Adjusted R-squared: 0.1625
## F-statistic: 11.94 on 5 and 277 DF, p-value: 1.776e-10
Anova(bc19, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 25521 1 24.9763 1.030e-06 ***
## bio19_12 3832 1 3.7506 0.0538 .
## Group 56035 4 13.7095 3.261e-10 ***
## Residuals 283044 277
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc19)
## [1] 283
### model B
bc19b<-lm(BCratio ~ bio19_12 + gPCA1+gPCA2+gPCA3+gPCA4+gPCA5, data=rp)
summary(bc19b)
##
## Call:
## lm(formula = BCratio ~ bio19_12 + gPCA1 + gPCA2 + gPCA3 + gPCA4 +
## gPCA5, data = rp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -77.15 -19.29 11.87 23.98 54.12
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 47.88087 6.21147 7.708 1.94e-13 ***
## bio19_12 0.12179 0.04444 2.740 0.006510 **
## gPCA1 -168.83018 60.86565 -2.774 0.005892 **
## gPCA2 -101.70364 40.19061 -2.531 0.011907 *
## gPCA3 -160.96162 41.49379 -3.879 0.000129 ***
## gPCA4 -35.38125 46.53140 -0.760 0.447636
## gPCA5 -49.58843 39.07503 -1.269 0.205418
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 33.3 on 296 degrees of freedom
## (259 observations deleted due to missingness)
## Multiple R-squared: 0.0921, Adjusted R-squared: 0.0737
## F-statistic: 5.004 on 6 and 296 DF, p-value: 6.664e-05
Anova(bc19b, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 65882 1 59.4203 1.944e-13 ***
## bio19_12 8326 1 7.5097 0.0065096 **
## gPCA1 8531 1 7.6941 0.0058922 **
## gPCA2 7100 1 6.4036 0.0119072 *
## gPCA3 16684 1 15.0480 0.0001292 ***
## gPCA4 641 1 0.5782 0.4476365
## gPCA5 1786 1 1.6105 0.2054178
## Residuals 328190 296
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc19b)
## [1] 303
### model C
bc19c<-lm(BCratio ~ bio19_12 + colutaPCA1+colutaPCA2+colutaPCA3+colutaPCA4+colutaPCA5, data=rp)
summary(bc19c)
##
## Call:
## lm(formula = BCratio ~ bio19_12 + colutaPCA1 + colutaPCA2 + colutaPCA3 +
## colutaPCA4 + colutaPCA5, data = rp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -68.11 -37.24 11.64 30.79 57.92
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 38.8541 18.6604 2.082 0.0394 *
## bio19_12 0.1172 0.1356 0.864 0.3894
## colutaPCA1 82.2193 37.1110 2.215 0.0285 *
## colutaPCA2 35.7888 40.5232 0.883 0.3788
## colutaPCA3 -108.1497 42.6199 -2.538 0.0124 *
## colutaPCA4 -39.1684 40.5749 -0.965 0.3362
## colutaPCA5 -21.9217 41.2226 -0.532 0.5958
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 36.83 on 125 degrees of freedom
## (430 observations deleted due to missingness)
## Multiple R-squared: 0.1501, Adjusted R-squared: 0.1093
## F-statistic: 3.68 on 6 and 125 DF, p-value: 0.002119
Anova(bc19c, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 5882 1 4.3354 0.03937 *
## bio19_12 1012 1 0.7460 0.38941
## colutaPCA1 6659 1 4.9084 0.02854 *
## colutaPCA2 1058 1 0.7800 0.37884
## colutaPCA3 8736 1 6.4391 0.01239 *
## colutaPCA4 1264 1 0.9319 0.33624
## colutaPCA5 384 1 0.2828 0.59582
## Residuals 169580 125
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc19c)
## [1] 132Validate p-values from linear models using permutation tests:
# Supplementary Table 11, sections A
set.seed(1)
nperm<-10000
var.list<-c('bio5_12', 'bio12_12', 'bio13_12', 'bio14_12', 'bio16_12', 'bio17_12', 'bio19_12')
mod.list<-list(bc5, bc12, bc13, bc14, bc16, bc17, bc19)
name.list<-c('Max. T in warmest month', 'Annual precipitation', 'Precip. in wettest month', 'Precip. in driest month', 'Precip. in wettest quarter', 'Precip. in driest quarter', 'Precip. in coldest quarter')
i<-1
j<-1
start<-Sys.time()
for (i in 1:length(var.list)){
perm.out.climate<-data.frame(matrix(NA, nrow=nperm, ncol=4))
colnames(perm.out.climate)<-c("Sum.Sq", "Df","F.value", "P")
perm.out.group<-perm.out.climate
for(j in 1:nperm){
perm.BC<-sample(clim.use$BCratio)
tmp<-cbind(clim.use, perm.BC)
model<-lm(paste('perm.BC ~', var.list[i], ' + Group'), data=tmp)
test<-data.frame(Anova(model, type="3"))
perm.out.climate[j,]<-test[2,]
perm.out.group[j,]<-test[3,]
}
# calculate empirical cdfs from permuted test statistics
pct.clim <- ecdf(perm.out.climate$F.value) # get CDF within null distribution
pct.group<-ecdf(perm.out.group$F.value)
# calculate observed test statistics from linear models
obs.F.clim<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=clim.use)), type=3)", sep="")))[2,3]
obs.F.group<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=clim.use)), type=3)", sep="")))[3,3]
# calculate p-values for climate and group from permutation tests
p.clim <-(1-pct.clim(obs.F.clim)) # 1 tailed test
p.group<-(1-pct.group(obs.F.group)) # 1 tailed test
# calculate location of 95% tail in climate F distribution
p.95.clim<-arrange(perm.out.climate, F.value)$F.value[0.95*nrow(perm.out.climate)]
# plot observed F onto null distribution
p<-ggplot(data=perm.out.climate) +
geom_histogram(aes(x = F.value),
color="lightgray", fill="lightgray", bins = 50) +
labs(title=paste(name.list[i], " (P = ",
sprintf("%0.4f", p.clim),")", sep=""),
x = "Permuted F-value", y = "Frequency") +
geom_segment(mapping=aes(x=obs.F.clim, y=900, xend=obs.F.clim, yend=50),
lineend = 'butt', linejoin = 'mitre',
arrow=arrow(length=unit(.1, "inches")),
size=.5, color="red") +
annotate("text", x=obs.F.clim, y=1000,
label="Observed F", hjust=0.5, color="red", size=4)+
geom_vline(xintercept=p.95.clim, colour="black", linetype="dashed")+
theme(plot.title = element_text(size=6))+
theme_bw()
# save plot
svg(filename = paste("SupFig4-PermTest-", var.list[i], ".svg", sep=""),
width = 4, height = 2, pointsize = 6)
print(p)
dev.off()
# print output for table
print(paste("Climate predictor: ", var.list[i], " (", name.list[i], ")", sep=""))
print(paste("Climate: F = ", sprintf("%0.4f", obs.F.clim), ", P = ", sprintf("%0.6f", p.clim), sep=""))
print(paste("Group: F = ", sprintf("%0.4f", obs.F.group), ", P = ", sprintf("%0.6f", p.group), sep=""))
print(paste("Elapsed time:", round((Sys.time()-start), digits=2)))
}
## [1] "Climate predictor: bio5_12 (Max. T in warmest month)"
## [1] "Climate: F = 1.4046, P = 0.230600"
## [1] "Group: F = 14.2150, P = 0.000000"
## [1] "Elapsed time: 1.59"
## [1] "Climate predictor: bio12_12 (Annual precipitation)"
## [1] "Climate: F = 12.6778, P = 0.000900"
## [1] "Group: F = 16.5453, P = 0.000000"
## [1] "Elapsed time: 3.2"
## [1] "Climate predictor: bio13_12 (Precip. in wettest month)"
## [1] "Climate: F = 10.5355, P = 0.001400"
## [1] "Group: F = 11.7070, P = 0.000000"
## [1] "Elapsed time: 4.81"
## [1] "Climate predictor: bio14_12 (Precip. in driest month)"
## [1] "Climate: F = 4.2606, P = 0.038800"
## [1] "Group: F = 14.8330, P = 0.000000"
## [1] "Elapsed time: 6.29"
## [1] "Climate predictor: bio16_12 (Precip. in wettest quarter)"
## [1] "Climate: F = 11.4703, P = 0.000500"
## [1] "Group: F = 12.6602, P = 0.000000"
## [1] "Elapsed time: 7.77"
## [1] "Climate predictor: bio17_12 (Precip. in driest quarter)"
## [1] "Climate: F = 5.2339, P = 0.022200"
## [1] "Group: F = 15.1799, P = 0.000000"
## [1] "Elapsed time: 9.25"
## [1] "Climate predictor: bio19_12 (Precip. in coldest quarter)"
## [1] "Climate: F = 3.7506, P = 0.054000"
## [1] "Group: F = 13.7095, P = 0.000000"
## [1] "Elapsed time: 10.72"# Supplementary Table 11, sections B
set.seed(1)
nperm<-10000
var.list<-c('bio5_12', 'bio12_12', 'bio13_12', 'bio14_12', 'bio16_12', 'bio17_12', 'bio19_12')
mod.list<-list(bc5b, bc12b, bc13b, bc14b, bc16b, bc17b, bc19b)
name.list<-c('Max. T in warmest month', 'Annual precipitation', 'Precip. in wettest month', 'Precip. in driest month', 'Precip. in wettest quarter', 'Precip. in driest quarter', 'Precip. in coldest quarter')
i<-1
j<-1
start<-Sys.time()
for (i in 1:length(var.list)){
perm.out.climate<-data.frame(matrix(NA, nrow=nperm, ncol=4))
colnames(perm.out.climate)<-c("Sum.Sq", "Df","F.value", "P")
perm.out.pca1<-perm.out.climate
perm.out.pca2<-perm.out.climate
perm.out.pca3<-perm.out.climate
perm.out.pca4<-perm.out.climate
perm.out.pca5<-perm.out.climate
for(j in 1:nperm){
perm.BC<-sample(rp$BCratio)
tmp<-cbind(rp, perm.BC)
model<-lm(paste('perm.BC ~', var.list[i], ' + gPCA1 + gPCA2 + gPCA3 + gPCA4 + gPCA5'), data=tmp)
test<-data.frame(Anova(model, type="3"))
perm.out.climate[j,]<-test[2,]
perm.out.pca1[j,]<-test[3,]
perm.out.pca2[j,]<-test[4,]
perm.out.pca3[j,]<-test[5,]
perm.out.pca4[j,]<-test[6,]
perm.out.pca5[j,]<-test[7,]
}
# calculate empirical cdfs from permuted test statistics
pct.clim <- ecdf(perm.out.climate$F.value) # get CDF within null distribution
pct.pca1<-ecdf(perm.out.pca1$F.value)
pct.pca2<-ecdf(perm.out.pca2$F.value)
pct.pca3<-ecdf(perm.out.pca3$F.value)
pct.pca4<-ecdf(perm.out.pca4$F.value)
pct.pca5<-ecdf(perm.out.pca5$F.value)
# calculate observed test statistics from linear models
obs.F.clim<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[2,3]
obs.F.pca1<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[3,3]
obs.F.pca2<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[4,3]
obs.F.pca3<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[5,3]
obs.F.pca4<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[6,3]
obs.F.pca5<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[7,3]
# calculate p-values for climate and group from permutation tests
p.clim <-(1-pct.clim(obs.F.clim)) # 1 tailed test
p.pca1<-(1-pct.pca1(obs.F.pca1)) # 1 tailed test
p.pca2<-(1-pct.pca2(obs.F.pca2)) # 1 tailed test
p.pca3<-(1-pct.pca3(obs.F.pca3)) # 1 tailed test
p.pca4<-(1-pct.pca4(obs.F.pca4)) # 1 tailed test
p.pca5<-(1-pct.pca5(obs.F.pca5)) # 1 tailed test
# calculate location of 95% tail in climate F distribution
p.95.clim<-arrange(perm.out.climate, F.value)$F.value[0.95*nrow(perm.out.climate)]
# plot observed F onto null distribution
p<-ggplot(data=perm.out.climate) +
geom_histogram(aes(x = F.value),
color="lightgray", fill="lightgray", bins = 50) +
labs(title=paste(name.list[i], " (P = ",
sprintf("%0.4f", p.clim),")", sep=""),
x = "Permuted F-value", y = "Frequency") +
geom_segment(mapping=aes(x=obs.F.clim, y=900, xend=obs.F.clim, yend=50),
lineend = 'butt', linejoin = 'mitre',
arrow=arrow(length=unit(.1, "inches")),
size=.5, color="red") +
annotate("text", x=obs.F.clim, y=1000,
label="Observed F", hjust=0.5, color="red", size=4)+
geom_vline(xintercept=p.95.clim, colour="black", linetype="dashed")+
theme(plot.title = element_text(size=6))+
theme_bw()
# save plot
svg(filename = paste("SupFig4-PermTest-", var.list[i], "-GeneticPCA-all.svg", sep=""),
width = 4, height = 2, pointsize = 6)
print(p)
dev.off()
# print output for table
print(paste("Climate predictor: ", var.list[i], " (", name.list[i], ")", sep=""))
print(paste("Climate: F = ", sprintf("%0.4f", obs.F.clim), ", P = ", sprintf("%0.6f", p.clim), sep=""))
print(paste("Genetic PCA1: F = ", sprintf("%0.4f", obs.F.pca1), ", P = ", sprintf("%0.6f", p.pca1), sep=""))
print(paste("Genetic PCA2: F = ", sprintf("%0.4f", obs.F.pca2), ", P = ", sprintf("%0.6f", p.pca2), sep=""))
print(paste("Genetic PCA3: F = ", sprintf("%0.4f", obs.F.pca3), ", P = ", sprintf("%0.6f", p.pca3), sep=""))
print(paste("Genetic PCA4: F = ", sprintf("%0.4f", obs.F.pca4), ", P = ", sprintf("%0.6f", p.pca4), sep=""))
print(paste("Genetic PCA5: F = ", sprintf("%0.4f", obs.F.pca5), ", P = ", sprintf("%0.6f", p.pca5), sep=""))
print(paste("Elapsed time:", round((Sys.time()-start), digits=2)))
}
## [1] "Climate predictor: bio5_12 (Max. T in warmest month)"
## [1] "Climate: F = 0.0960, P = 0.759100"
## [1] "Genetic PCA1: F = 7.3395, P = 0.006800"
## [1] "Genetic PCA2: F = 3.8367, P = 0.049900"
## [1] "Genetic PCA3: F = 11.5296, P = 0.000700"
## [1] "Genetic PCA4: F = 0.0398, P = 0.840900"
## [1] "Genetic PCA5: F = 1.6247, P = 0.197400"
## [1] "Elapsed time: 2.58"
## [1] "Climate predictor: bio12_12 (Annual precipitation)"
## [1] "Climate: F = 9.3057, P = 0.002600"
## [1] "Genetic PCA1: F = 8.8982, P = 0.002700"
## [1] "Genetic PCA2: F = 8.7538, P = 0.003100"
## [1] "Genetic PCA3: F = 15.7124, P = 0.000200"
## [1] "Genetic PCA4: F = 0.4277, P = 0.510800"
## [1] "Genetic PCA5: F = 1.5554, P = 0.210900"
## [1] "Elapsed time: 5.21"
## [1] "Climate predictor: bio13_12 (Precip. in wettest month)"
## [1] "Climate: F = 23.0087, P = 0.000000"
## [1] "Genetic PCA1: F = 8.8180, P = 0.003200"
## [1] "Genetic PCA2: F = 9.5028, P = 0.001900"
## [1] "Genetic PCA3: F = 9.3895, P = 0.002500"
## [1] "Genetic PCA4: F = 1.6860, P = 0.197700"
## [1] "Genetic PCA5: F = 1.8447, P = 0.178900"
## [1] "Elapsed time: 7.88"
## [1] "Climate predictor: bio14_12 (Precip. in driest month)"
## [1] "Climate: F = 0.0249, P = 0.872400"
## [1] "Genetic PCA1: F = 7.2725, P = 0.007800"
## [1] "Genetic PCA2: F = 3.6262, P = 0.063600"
## [1] "Genetic PCA3: F = 10.2724, P = 0.002000"
## [1] "Genetic PCA4: F = 0.0291, P = 0.866900"
## [1] "Genetic PCA5: F = 1.5446, P = 0.214500"
## [1] "Elapsed time: 10.48"
## [1] "Climate predictor: bio16_12 (Precip. in wettest quarter)"
## [1] "Climate: F = 23.0253, P = 0.000000"
## [1] "Genetic PCA1: F = 9.2846, P = 0.002700"
## [1] "Genetic PCA2: F = 9.4471, P = 0.002100"
## [1] "Genetic PCA3: F = 12.9361, P = 0.000500"
## [1] "Genetic PCA4: F = 1.8791, P = 0.168600"
## [1] "Genetic PCA5: F = 1.9002, P = 0.162300"
## [1] "Elapsed time: 13.16"
## [1] "Climate predictor: bio17_12 (Precip. in driest quarter)"
## [1] "Climate: F = 0.0070, P = 0.932700"
## [1] "Genetic PCA1: F = 7.2477, P = 0.008100"
## [1] "Genetic PCA2: F = 3.4786, P = 0.065300"
## [1] "Genetic PCA3: F = 10.6703, P = 0.001400"
## [1] "Genetic PCA4: F = 0.0311, P = 0.866100"
## [1] "Genetic PCA5: F = 1.5683, P = 0.210900"
## [1] "Elapsed time: 15.79"
## [1] "Climate predictor: bio19_12 (Precip. in coldest quarter)"
## [1] "Climate: F = 7.5097, P = 0.007100"
## [1] "Genetic PCA1: F = 7.6941, P = 0.005300"
## [1] "Genetic PCA2: F = 6.4036, P = 0.011400"
## [1] "Genetic PCA3: F = 15.0480, P = 0.000200"
## [1] "Genetic PCA4: F = 0.5782, P = 0.455600"
## [1] "Genetic PCA5: F = 1.6105, P = 0.214400"
## [1] "Elapsed time: 18.38"# Supplementary Table 11, sections C
set.seed(1)
nperm<-10000
var.list<-c('bio5_12', 'bio12_12', 'bio13_12', 'bio14_12', 'bio16_12', 'bio17_12', 'bio19_12')
mod.list<-list(bc5c, bc12c, bc13c, bc14c, bc16c, bc17c, bc19c)
name.list<-c('Max. T in warmest month', 'Annual precipitation', 'Precip. in wettest month', 'Precip. in driest month', 'Precip. in wettest quarter', 'Precip. in driest quarter', 'Precip. in coldest quarter')
i<-1
j<-1
start<-Sys.time()
for (i in 1:length(var.list)){
perm.out.climate<-data.frame(matrix(NA, nrow=nperm, ncol=4))
colnames(perm.out.climate)<-c("Sum.Sq", "Df","F.value", "P")
perm.out.pca1<-perm.out.climate
perm.out.pca2<-perm.out.climate
perm.out.pca3<-perm.out.climate
perm.out.pca4<-perm.out.climate
perm.out.pca5<-perm.out.climate
for(j in 1:nperm){
perm.BC<-sample(rp$BCratio)
tmp<-cbind(rp, perm.BC)
model<-lm(paste('perm.BC ~', var.list[i], ' + colutaPCA1 + colutaPCA2 + colutaPCA3 + colutaPCA4 + colutaPCA5'), data=tmp)
test<-data.frame(Anova(model, type="3"))
perm.out.climate[j,]<-test[2,]
perm.out.pca1[j,]<-test[3,]
perm.out.pca2[j,]<-test[4,]
perm.out.pca3[j,]<-test[5,]
perm.out.pca4[j,]<-test[6,]
perm.out.pca5[j,]<-test[7,]
}
# calculate empirical cdfs from permuted test statistics
pct.clim <- ecdf(perm.out.climate$F.value) # get CDF within null distribution
pct.pca1<-ecdf(perm.out.pca1$F.value)
pct.pca2<-ecdf(perm.out.pca2$F.value)
pct.pca3<-ecdf(perm.out.pca3$F.value)
pct.pca4<-ecdf(perm.out.pca4$F.value)
pct.pca5<-ecdf(perm.out.pca5$F.value)
# calculate observed test statistics from linear models
obs.F.clim<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[2,3]
obs.F.pca1<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[3,3]
obs.F.pca2<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[4,3]
obs.F.pca3<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[5,3]
obs.F.pca4<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[6,3]
obs.F.pca5<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[7,3]
# calculate p-values for climate and group from permutation tests
p.clim <-(1-pct.clim(obs.F.clim)) # 1 tailed test
p.pca1<-(1-pct.pca1(obs.F.pca1)) # 1 tailed test
p.pca2<-(1-pct.pca2(obs.F.pca2)) # 1 tailed test
p.pca3<-(1-pct.pca3(obs.F.pca3)) # 1 tailed test
p.pca4<-(1-pct.pca4(obs.F.pca4)) # 1 tailed test
p.pca5<-(1-pct.pca5(obs.F.pca5)) # 1 tailed test
# calculate location of 95% tail in climate F distribution
p.95.clim<-arrange(perm.out.climate, F.value)$F.value[0.95*nrow(perm.out.climate)]
# plot observed F onto null distribution
p<-ggplot(data=perm.out.climate) +
geom_histogram(aes(x = F.value),
color="lightgray", fill="lightgray", bins = 50) +
labs(title=paste(name.list[i], " (P = ",
sprintf("%0.4f", p.clim),")", sep=""),
x = "Permuted F-value", y = "Frequency") +
geom_segment(mapping=aes(x=obs.F.clim, y=900, xend=obs.F.clim, yend=50),
lineend = 'butt', linejoin = 'mitre',
arrow=arrow(length=unit(.1, "inches")),
size=.5, color="red") +
annotate("text", x=obs.F.clim, y=1000,
label="Observed F", hjust=0.5, color="red", size=4)+
geom_vline(xintercept=p.95.clim, colour="black", linetype="dashed")+
theme(plot.title = element_text(size=6))+
theme_bw()
# save plot
svg(filename = paste("SupFig4-PermTest-", var.list[i], "-GeneticPCA-coluta.svg", sep=""),
width = 4, height = 2, pointsize = 6)
print(p)
dev.off()
# print output for table
print(paste("Climate predictor: ", var.list[i], " (", name.list[i], ")", sep=""))
print(paste("Climate: F = ", sprintf("%0.4f", obs.F.clim), ", P = ", sprintf("%0.6f", p.clim), sep=""))
print(paste("Genetic PCA1: F = ", sprintf("%0.4f", obs.F.pca1), ", P = ", sprintf("%0.6f", p.pca1), sep=""))
print(paste("Genetic PCA2: F = ", sprintf("%0.4f", obs.F.pca2), ", P = ", sprintf("%0.6f", p.pca2), sep=""))
print(paste("Genetic PCA3: F = ", sprintf("%0.4f", obs.F.pca3), ", P = ", sprintf("%0.6f", p.pca3), sep=""))
print(paste("Genetic PCA4: F = ", sprintf("%0.4f", obs.F.pca4), ", P = ", sprintf("%0.6f", p.pca4), sep=""))
print(paste("Genetic PCA5: F = ", sprintf("%0.4f", obs.F.pca5), ", P = ", sprintf("%0.6f", p.pca5), sep=""))
print(paste("Elapsed time:", round((Sys.time()-start), digits=2)))
}
## [1] "Climate predictor: bio5_12 (Max. T in warmest month)"
## [1] "Climate: F = 1.0683, P = 0.303800"
## [1] "Genetic PCA1: F = 6.0656, P = 0.017600"
## [1] "Genetic PCA2: F = 1.8505, P = 0.180500"
## [1] "Genetic PCA3: F = 5.8068, P = 0.015100"
## [1] "Genetic PCA4: F = 2.6594, P = 0.102000"
## [1] "Genetic PCA5: F = 0.0362, P = 0.846100"
## [1] "Elapsed time: 2.6"
## [1] "Climate predictor: bio12_12 (Annual precipitation)"
## [1] "Climate: F = 8.4704, P = 0.004800"
## [1] "Genetic PCA1: F = 6.3665, P = 0.012700"
## [1] "Genetic PCA2: F = 0.4289, P = 0.521200"
## [1] "Genetic PCA3: F = 1.9346, P = 0.166400"
## [1] "Genetic PCA4: F = 0.2824, P = 0.599600"
## [1] "Genetic PCA5: F = 0.5192, P = 0.469600"
## [1] "Elapsed time: 5.24"
## [1] "Climate predictor: bio13_12 (Precip. in wettest month)"
## [1] "Climate: F = 13.9423, P = 0.000200"
## [1] "Genetic PCA1: F = 4.1154, P = 0.047200"
## [1] "Genetic PCA2: F = 2.9675, P = 0.091100"
## [1] "Genetic PCA3: F = 2.6294, P = 0.102600"
## [1] "Genetic PCA4: F = 0.1631, P = 0.681700"
## [1] "Genetic PCA5: F = 0.5544, P = 0.454700"
## [1] "Elapsed time: 7.84"
## [1] "Climate predictor: bio14_12 (Precip. in driest month)"
## [1] "Climate: F = 1.4031, P = 0.242600"
## [1] "Genetic PCA1: F = 6.7849, P = 0.011000"
## [1] "Genetic PCA2: F = 0.0451, P = 0.829800"
## [1] "Genetic PCA3: F = 7.1233, P = 0.007000"
## [1] "Genetic PCA4: F = 2.4670, P = 0.120800"
## [1] "Genetic PCA5: F = 0.0195, P = 0.890000"
## [1] "Elapsed time: 10.43"
## [1] "Climate predictor: bio16_12 (Precip. in wettest quarter)"
## [1] "Climate: F = 15.3993, P = 0.000100"
## [1] "Genetic PCA1: F = 6.5759, P = 0.013100"
## [1] "Genetic PCA2: F = 0.9184, P = 0.347400"
## [1] "Genetic PCA3: F = 3.2619, P = 0.067700"
## [1] "Genetic PCA4: F = 0.0067, P = 0.935700"
## [1] "Genetic PCA5: F = 0.4130, P = 0.528100"
## [1] "Elapsed time: 13.11"
## [1] "Climate predictor: bio17_12 (Precip. in driest quarter)"
## [1] "Climate: F = 0.6609, P = 0.423000"
## [1] "Genetic PCA1: F = 5.8596, P = 0.016600"
## [1] "Genetic PCA2: F = 0.8065, P = 0.367500"
## [1] "Genetic PCA3: F = 5.8679, P = 0.017100"
## [1] "Genetic PCA4: F = 1.8993, P = 0.165300"
## [1] "Genetic PCA5: F = 0.0495, P = 0.823700"
## [1] "Elapsed time: 15.7"
## [1] "Climate predictor: bio19_12 (Precip. in coldest quarter)"
## [1] "Climate: F = 0.7460, P = 0.394400"
## [1] "Genetic PCA1: F = 4.9084, P = 0.029900"
## [1] "Genetic PCA2: F = 0.7800, P = 0.384100"
## [1] "Genetic PCA3: F = 6.4391, P = 0.009700"
## [1] "Genetic PCA4: F = 0.9319, P = 0.342800"
## [1] "Genetic PCA5: F = 0.2828, P = 0.585900"
## [1] "Elapsed time: 18.3"Structural variants in BCMA3 and transcribed amino acid sequences were determined from .fasta sequence data using scripts in ‘/BCMA archive.R1/Sanger data/Sequence analysis/sanger_seq_analysis-v2.Rmd’.
Distribution of and covariance among alleles:
# Contingency tables of alleles
table(rp$Genotype, rp$Group)
##
## Admixed COL NOR UTA WES
## Complete exons + No deletion 5 18 18 14 22
## Premature stop + 8bp deletion 1 16 0 3 0
## Premature stop + No deletion 0 2 0 1 0
table(rp$aa.concat, rp$Group)
##
## Admixed COL NOR UTA WES
## ** 1 16 0 3 0
## LV 2 10 3 5 0
## VM 0 0 2 0 14
## VV 3 10 13 10 8
table(rp$Genotype, rp$aa.concat, rp$Group)
## , , = Admixed
##
##
## ** LV VM VV
## Complete exons + No deletion 0 2 0 3
## Premature stop + 8bp deletion 1 0 0 0
## Premature stop + No deletion 0 0 0 0
##
## , , = COL
##
##
## ** LV VM VV
## Complete exons + No deletion 0 8 0 10
## Premature stop + 8bp deletion 16 0 0 0
## Premature stop + No deletion 0 2 0 0
##
## , , = NOR
##
##
## ** LV VM VV
## Complete exons + No deletion 0 3 2 13
## Premature stop + 8bp deletion 0 0 0 0
## Premature stop + No deletion 0 0 0 0
##
## , , = UTA
##
##
## ** LV VM VV
## Complete exons + No deletion 0 4 0 10
## Premature stop + 8bp deletion 3 0 0 0
## Premature stop + No deletion 0 1 0 0
##
## , , = WES
##
##
## ** LV VM VV
## Complete exons + No deletion 0 0 14 8
## Premature stop + 8bp deletion 0 0 0 0
## Premature stop + No deletion 0 0 0 0
table(rp$Genotype, rp$aa.concat)
##
## ** LV VM VV
## Complete exons + No deletion 0 18 18 51
## Premature stop + 8bp deletion 20 0 0 0
## Premature stop + No deletion 0 3 0 0
table(rp$Deletion, rp$Full_Length)
##
## Complete exons Premature stop
## 8bp deletion 0 20
## No deletion 87 3
# Measure geographic distances of alleles
full.length<-rp %>% filter(Genotype=="Complete exons + No deletion") %>%dplyr::select(Longitude, Latitude)
stop.only<-rp %>% filter(Genotype =="Premature stop + No deletion") %>% dplyr::select(Latitude, Longitude)
stop.del<-rp%>% filter(Genotype=="Premature stop + 8bp deletion") %>% dplyr::select(Latitude, Longitude)
georange(full.length)/1000
## minimum maximum
## 0.3648679 1216.6627968
georange(stop.only)/1000
## minimum maximum
## 104.7254 328.2988
georange(stop.del)/1000
## minimum maximum
## 0.4168024 729.0933943Effects of BCMA3 variants on BC-ratio:
library(car)
# Data subset with only full-length alleles
final.fl<-filter(rp, Full_Length=="Complete exons")
# Supplementary Table 12A:
met.sep.1<-lm(BCratio~Group+Full_Length+Deletion, data=rp)
summary(met.sep.1)
##
## Call:
## lm(formula = BCratio ~ Group + Full_Length + Deletion, data = rp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -70.522 -2.136 3.575 15.882 35.852
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 70.27 30.34 2.316 0.02327 *
## GroupCOL 10.92 12.57 0.869 0.38777
## GroupNOR 11.44 12.64 0.905 0.36858
## GroupUTA 16.58 12.98 1.277 0.20558
## GroupWES 35.37 12.36 2.862 0.00546 **
## Full_LengthPremature stop -71.98 27.57 -2.610 0.01091 *
## DeletionNo deletion -13.29 28.74 -0.463 0.64504
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26.61 on 75 degrees of freedom
## (480 observations deleted due to missingness)
## Multiple R-squared: 0.5283, Adjusted R-squared: 0.4905
## F-statistic: 14 on 6 and 75 DF, p-value: 1.288e-10
Anova(met.sep.1, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 3798 1 5.3655 0.02327 *
## Group 9763 4 3.4483 0.01215 *
## Full_Length 4823 1 6.8141 0.01091 *
## Deletion 151 1 0.2139 0.64504
## Residuals 53089 75
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(met.sep.1)
## [1] 82
# Supplementary Table 12B:
met.sep.2<-lm(BCratio~Group+aa.concat, data=final.fl)
summary(met.sep.2)
##
## Call:
## lm(formula = BCratio ~ Group + aa.concat, data = final.fl)
##
## Residuals:
## Min 1Q Median 3Q Max
## -70.740 -1.448 4.764 14.395 43.028
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 48.123 13.942 3.452 0.00101 **
## GroupCOL 12.243 15.225 0.804 0.42440
## GroupNOR 8.441 14.638 0.577 0.56626
## GroupUTA 17.877 15.209 1.175 0.24433
## GroupWES 26.352 16.124 1.634 0.10725
## aa.concatVM 20.502 14.056 1.459 0.14973
## aa.concatVV 13.251 9.336 1.419 0.16081
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 28.55 on 62 degrees of freedom
## (18 observations deleted due to missingness)
## Multiple R-squared: 0.1913, Adjusted R-squared: 0.1131
## F-statistic: 2.445 on 6 and 62 DF, p-value: 0.0347
Anova(met.sep.2, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 9710 1 11.9137 0.001009 **
## Group 3368 4 1.0331 0.397388
## aa.concat 2070 2 1.2698 0.288089
## Residuals 50532 62
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(met.sep.2)
## [1] 69
# Supplementary Table 12C:
met.sep.3<-lm(BCratio~Group+AA148+AA268, data=final.fl)
summary(met.sep.3)
##
## Call:
## lm(formula = BCratio ~ Group + AA148 + AA268, data = final.fl)
##
## Residuals:
## Min 1Q Median 3Q Max
## -70.740 -1.448 4.764 14.395 43.028
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 55.373 17.842 3.104 0.00288 **
## GroupCOL 12.243 15.225 0.804 0.42440
## GroupNOR 8.441 14.638 0.577 0.56626
## GroupUTA 17.877 15.209 1.175 0.24433
## GroupWES 26.352 16.124 1.634 0.10725
## AA148V 13.251 9.336 1.419 0.16081
## AA268V -7.250 10.902 -0.665 0.50851
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 28.55 on 62 degrees of freedom
## (18 observations deleted due to missingness)
## Multiple R-squared: 0.1913, Adjusted R-squared: 0.1131
## F-statistic: 2.445 on 6 and 62 DF, p-value: 0.0347
Anova(met.sep.3, type="3")
## Anova Table (Type III tests)
##
## Response: BCratio
## Sum Sq Df F value Pr(>F)
## (Intercept) 7851 1 9.6324 0.002879 **
## Group 3368 4 1.0331 0.397388
## AA148 1642 1 2.0145 0.160808
## AA268 360 1 0.4422 0.508505
## Residuals 50532 62
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(met.sep.3)
## [1] 69FST was calculated using methods and scripts described in Wang et al. 2019.
Pairwise LD among SNPs in the BCMA CFR-HIF region was calculated using the software PLINK (Purcell et al. 2007). PLINK output was modified for plotting heatmaps using a custom script (ld2matrix.py).
Pi for BCMA3 and other genes in the B. stricta genome were calculated using DnaSP, after correcting for possible sequence errors in regions supported by a single read (Python scripts: calculate_per_site_div.lable_singleton.py and calculate_per_window_diversity.masked_singleton.py).
Filter pi data:
# designate observed pi BCMA3 using DnaSP
pi.BCMA3 <- 0.00591
# Only compare pi to genes with a similar length to BCMA3, since length is correlated with nucleotide diversity (Supplementary Figure 9). Also, filter out genes with pi = 0.
min.num.sites <- 360
nuc.div.OK <- filter(nuc.div, Pi >= 0 & site >= min.num.sites)
# Three main regions of the B. stricta genome have previously been identified as having elevated nucleotide diversity due to historical balancing selection (Wang et al. 2019). We also exclude these regions from comparisons testing for balancing selection on BCMA3, since we already know they have experienced balancing selection.
## by default, assign all genes to the region "background"
nuc.div.OK$reg <- "background"
## flag the regions of elevated polymorphism on LG1, LG5, and LG7 as identified in Wang et al. 2019
nuc.div.OK$reg[nuc.div.OK$start >= 14423807 &
nuc.div.OK$end <= 24375392 & nuc.div.OK$LG ==1] <- "LG1"
nuc.div.OK$reg[nuc.div.OK$start >= 471982 &
nuc.div.OK$end <= 3400076 & nuc.div.OK$LG == 5] <- "LG5"
nuc.div.OK$reg[nuc.div.OK$start >= 451506 &
nuc.div.OK$end <= 5083168 & nuc.div.OK$LG == 7] <- "LG7"
## filter to focus only on background regions
nuc.div.OK.2 <- filter(nuc.div.OK, reg == "background")
nrow(nuc.div.OK.2) # this is the number of comparable genes retained
## [1] 1689Compare BCMA to the genome-wide distribution of pi:
# calculate the proportion of comparable genes for which pi > BCMA3
x <- nuc.div.OK.2$Pi > pi.BCMA3
p.val <- mean(x)
p.val
## [1] 0.02901125Average pi in comparable genes:
mean(nuc.div.OK.2$Pi)
## [1] 0.001690409Get mean and SD of herbivore damage in each transplant environment:
# subset to transplant experiments
field.use<-filter(field, Expt=="transplants")
# mean herbivory for each geno in each env
h.means<-field.use %>% dplyr::select(Geno, Env, LAR2) %>% drop_na() %>%
group_by(Geno, Env) %>%
summarise(Mean_dam=mean(LAR2))
# reshape means so there's one column for each geno
h.means<-spread(h.means, Geno, Mean_dam)
# SD of herbivory for each geno in each env
h.sd<-field.use%>%dplyr::select(Geno, Env, LAR2) %>% drop_na() %>%
group_by(Env) %>%
summarise(SD_dam=sd(LAR2))Calculate effect size:
# vector of differences in genotypic means/SD per environment
## initialize vector
diff.vec<-rep(NA, nrow(h.sd))
## calculate effect size in each environment
i<-1
for(i in 1:length(diff.vec)){
mean.LL<-h.means$LL[i]
mean.SS<-h.means$SS[i]
stdev<-h.sd$SD_dam[i]
std.diff<-abs(mean.LL-mean.SS)/stdev
diff.vec[i]<-std.diff
}
# calculate mean effect size across environments
mean(diff.vec)
## [1] 0.1715397
# calculate min and max effect sizes
min(diff.vec)
## [1] 0.01260256
max(diff.vec)
## [1] 0.4329623
# view differences by environment
data.frame(cbind(h.means$Env, as.numeric(diff.vec)))
## X1 X2
## 1 ALD-2013 0.3853977028103
## 2 ALD-2014 0.0594271726528837
## 3 ALD-2015 0.0854734188355042
## 4 GTH-2016 0.432962264183009
## 5 MAH-2014 0.0126025590252014
## 6 MAH-2015 0.15228075876984
## 7 SCH-2016 0.144273660559002
## 8 SIL-2014 0.0771454958018908
## 9 SIL-2015 0.194294400846806Package versions, etc.
sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] geodist_0.0.4 ggplot2_3.3.2 pbkrtest_0.4-8.6 emmeans_1.5.0
## [5] car_3.0-9 carData_3.0-4 lmerTest_3.1-2 lme4_1.1-23
## [9] Matrix_1.2-18 tidyr_1.1.1 dplyr_1.0.2
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.5 mvtnorm_1.1-1 lattice_0.20-41
## [4] zoo_1.8-8 digest_0.6.25 plyr_1.8.6
## [7] R6_2.4.1 cellranger_1.1.0 evaluate_0.14
## [10] coda_0.19-3 pillar_1.4.6 rlang_0.4.7
## [13] curl_4.3 multcomp_1.4-13 readxl_1.3.1
## [16] rstudioapi_0.11 minqa_1.2.4 data.table_1.13.0
## [19] nloptr_1.2.2.2 rmarkdown_2.3 labeling_0.3
## [22] splines_4.0.2 statmod_1.4.34 stringr_1.4.0
## [25] foreign_0.8-80 munsell_0.5.0 compiler_4.0.2
## [28] numDeriv_2016.8-1.1 xfun_0.16 pkgconfig_2.0.3
## [31] htmltools_0.5.0 tidyselect_1.1.0 tibble_3.0.3
## [34] rio_0.5.16 codetools_0.2-16 crayon_1.3.4
## [37] withr_2.2.0 MASS_7.3-52 grid_4.0.2
## [40] nlme_3.1-148 xtable_1.8-4 gtable_0.3.0
## [43] lifecycle_0.2.0 magrittr_1.5 scales_1.1.1
## [46] zip_2.0.4 estimability_1.3 stringi_1.4.6
## [49] farver_2.0.3 ellipsis_0.3.1 generics_0.0.2
## [52] vctrs_0.3.2 boot_1.3-25 sandwich_2.5-1
## [55] openxlsx_4.1.5 TH.data_1.0-10 tools_4.0.2
## [58] forcats_0.5.0 glue_1.4.1 purrr_0.3.4
## [61] hms_0.5.3 abind_1.4-5 survival_3.2-3
## [64] yaml_2.2.1 colorspace_1.4-1 knitr_1.29
## [67] haven_2.3.1