This notebook explores some of the data analyzed in Carlson et al., Annals of Botany 117:195-207; 2016 (doi: [https://dx.doi.org/10.1093/aob/mcv146]). Specifically, the focus here is on trait-environment relationships using data measured on seedlings in the Kirstenbosch Botanical Garden in 2013 and 2014. The data file here is the same one (or should be the same one) available at [https://github.com/kholsinger/Protea-repens-physiology/releases/tag/v1.0].
## clear memory of any existing analyses to avoid conflicts
##
rm(list=ls())
dat <- read.csv("traits-environment-pca.csv", header=TRUE, na.strings=c("NA","."))This analysis will use only a subset of the data, so the first step is to reduce the data frame to include only the relevant columns.
dat <- dat[,c("stomatal_density","Prin1_temp","Prin2_dry","Prin3_map","year")]
summary(dat) stomatal_density Prin1_temp Prin2_dry Prin3_map year
Min. :22.22 Min. :-2.0057 Min. :-1.65822 Min. :-1.59709 Min. :2013
1st Qu.:37.04 1st Qu.:-0.9930 1st Qu.:-0.57682 1st Qu.:-0.84424 1st Qu.:2013
Median :44.44 Median : 0.3598 Median : 0.00273 Median :-0.02303 Median :2014
Mean :45.19 Mean :-0.1211 Mean : 0.13461 Mean :-0.02465 Mean :2014
3rd Qu.:51.85 3rd Qu.: 0.7217 3rd Qu.: 1.34120 3rd Qu.: 0.56104 3rd Qu.:2014
Max. :79.63 Max. : 1.2202 Max. : 1.51035 Max. : 2.05629 Max. :2014
NA's :1130 NA's :2 NA's :2 NA's :2 Stomatal density will be our response variable in the analysis, so there’s no reason to retain any rows where stomatal density is NA. Year should also be treated as a factor, rather than as a numeric value.
dat <- subset(dat, !is.na(stomatal_density))
dat$year <- factor(dat$year)
summary(dat) stomatal_density Prin1_temp Prin2_dry Prin3_map year
Min. :22.22 Min. :-2.00573 Min. :-1.65822 Min. :-1.59709 2013:233
1st Qu.:37.04 1st Qu.:-0.75779 1st Qu.:-0.59996 1st Qu.:-0.36235 2014:277
Median :44.44 Median : 0.50219 Median :-0.08574 Median :-0.02303
Mean :45.19 Mean : 0.03688 Mean : 0.04012 Mean : 0.04309
3rd Qu.:51.85 3rd Qu.: 0.73215 3rd Qu.: 0.84284 3rd Qu.: 0.56104
Max. :79.63 Max. : 1.22024 Max. : 1.51035 Max. : 2.05629 Carlson et al. reported a positive relationship between principal component axis 1 and stomatal density and a negative relationship between axis 2 and stomatal density.
library(rstanarm)Loading required package: Rcpp
rstanarm (Version 2.17.2, packaged: 2017-12-20 23:59:28 UTC)
- Do not expect the default priors to remain the same in future rstanarm versions.
Thus, R scripts should specify priors explicitly, even if they are just the defaults.
- For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores())
- Plotting theme set to bayesplot::theme_default().options(mc.cores = parallel::detectCores())
dat_lm <- stan_lmer(stomatal_density ~ Prin1_temp + Prin2_dry + Prin3_map + (1|year),
data=dat,
adapt_delta=0.999)starting worker pid=28686 on localhost:11903 at 20:05:35.555
starting worker pid=28694 on localhost:11903 at 20:05:35.735
starting worker pid=28702 on localhost:11903 at 20:05:35.909
starting worker pid=28710 on localhost:11903 at 20:05:36.082
Manually removed progress report
opt_old <- options(width=180)
summary(dat_lm, digits=3)
Model Info:
function: stan_lmer
family: gaussian [identity]
formula: stomatal_density ~ Prin1_temp + Prin2_dry + Prin3_map + (1 |
year)
algorithm: sampling
priors: see help('prior_summary')
sample: 4000 (posterior sample size)
observations: 510
groups: year (2)
Estimates:
mean sd 2.5% 25% 50% 75% 97.5%
(Intercept) 45.068 8.994 26.685 42.254 45.229 48.192 62.266
Prin1_temp 2.422 0.420 1.597 2.137 2.423 2.712 3.216
Prin2_dry -2.125 0.428 -2.980 -2.410 -2.117 -1.851 -1.277
Prin3_map 1.317 0.399 0.538 1.041 1.314 1.585 2.099
b[(Intercept) year:2013] 3.449 8.998 -13.696 0.357 3.262 6.259 21.890
b[(Intercept) year:2014] -2.791 8.994 -19.903 -5.932 -2.967 0.046 15.220
sigma 8.865 0.291 8.309 8.665 8.854 9.063 9.457
Sigma[year:(Intercept),(Intercept)] 125.560 316.250 4.211 16.411 42.338 112.450 746.986
mean_PPD 45.189 0.549 44.137 44.829 45.185 45.564 46.270
log-posterior -1851.093 2.270 -1856.443 -1852.358 -1850.746 -1849.412 -1847.772
Diagnostics:
mcse Rhat n_eff
(Intercept) 0.450 1.009 400
Prin1_temp 0.007 0.999 4000
Prin2_dry 0.007 1.001 4000
Prin3_map 0.006 1.001 4000
b[(Intercept) year:2013] 0.450 1.008 400
b[(Intercept) year:2014] 0.450 1.009 400
sigma 0.005 1.000 4000
Sigma[year:(Intercept),(Intercept)] 10.479 1.002 911
mean_PPD 0.009 1.001 4000
log-posterior 0.059 1.002 1489
For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).options(opt_old)This analysis detects the same relationships Carlson et al. reported, but it also detects a relationship with axis 2 that they did not detect. The difference arises because the Carlson et al. model treated stomatal density as one component of a vector including other traits. Presumably the relationship between stomatal density and axis 2 is detected here because stomatal density is negatively associated with leaf area. The difference in results suggests that the association between axis 2 and stomatal density detected here is not causal.