library("PEcAn.MA")
library("tibble", include.only = "tribble")
set.seed(20251222)An example of a prior distributions table, demonstrating the expected format.
prior_distns <- tribble(
~row_name , ~distn , ~parama , ~paramb , ~n , #nolint
"growth_resp_factor" , "beta" , 2.63 , 6.52 , 0 , #nolint
"leaf_turnover_rate" , "weibull" , 1.37 , 1.43 , 363 , #nolint
"root_respiration_rate" , "unif" , 0.00 , 100.00 , NA , #nolint
"root_turnover_rate" , "unif" , 0.00 , 10.00 , NA , #nolint
"Amax" , "unif" , 0.00 , 40.00 , NA , #nolint
"leaf_respiration_rate_m2" , "weibull" , 2.00 , 6.00 , NA , #nolint
"SLA" , "lnorm" , 1.89 , 0.61 , 455 , #nolint
"leafC" , "norm" , 50.60 , 1.32 , 291 , #nolint
"Vm_low_temp" , "norm" , 0.00 , 3.00 , NA , #nolint
"AmaxFrac" , "unif" , 0.60 , 0.90 , NA , #nolint
"psnTOpt" , "unif" , 5.00 , 40.00 , NA , #nolint
"stem_respiration_rate" , "unif" , 0.00 , 100.00 , NA , #nolint
"extinction_coefficient" , "unif" , 0.38 , 0.62 , NA , #nolint
"half_saturation_PAR" , "unif" , 4.00 , 27.00 , NA , #nolint
"dVPDSlope" , "unif" , 0.01 , 0.25 , NA , #nolint
"dVpdExp" , "unif" , 1.00 , 3.00 , NA , #nolint
"veg_respiration_Q10" , "unif" , 1.40 , 2.60 , NA , #nolint
"fine_root_respiration_Q10" , "unif" , 1.40 , 5.00 , NA , #nolint
"coarse_root_respiration_Q10" , "unif" , 1.40 , 5.00 , NA , #nolint
)PEcAn.MA expects the code to be a base R
data.frame with trait names as row names.
priors <- as.data.frame(prior_distns)
rownames(priors) <- priors[["row_name"]]
priors[["row_name"]] <- NULLThe format of the trait_data is a named list of
data.frames, with each data.frame containing
data for the corresponding trait. PEcAn database queries return all the
columns the meta-analysis needs, in the expected format. If you are
providing your own trait data, you will need the following columns:
-
name(character) – A string description of the trait -
mean(numeric) – The mean value of the trait measurement (or the only value, if only one value is given) -
statname(character) – error statistic type (e.g., “SE” for standard error) -
stat(numeric) – value of error statistic -
n(integer or NA) – sample size -
greenhouse(boolean) –TRUEif from a greenhouse;FALSEif not (e.g., natural setting) -
control(boolean) – IfTRUE, this is the “control” part of an experiment (or there is no experiment). IfFALSE, this is the treatment. -
site_id(integer) – site ID (used for grouping) -
specie_id(integer) – species ID (as above) -
citation_id(integer) – citation ID (as above) -
treatment_id(integer) – treatment ID (as above) -
cultivar_id(integer / NA) – cultivar ID (as above) -
date(datetime / NA) – date of trait observation (as above) -
time(datetime / NA) – time of trait observation (as above)
To avoid having to load custom data, here is some simulated data that fits these criteria.
simulate_data <- function(n_rows, mean_lo, mean_hi, se_lo, se_hi) {
return(data.frame(
name = sprintf("t%02d", seq_len(n_rows)),
mean = runif(n_rows, mean_lo, mean_hi),
statname = "SE",
stat = suppressWarnings(runif(n_rows, se_lo, se_hi)),
n = sample(1:100, size = n_rows, replace = TRUE),
greenhouse = sample(
c(TRUE, FALSE),
size = n_rows,
replace = TRUE,
prob = c(0.1, 0.9)
),
control = TRUE,
site_id = sample(1:10, size = n_rows, replace = TRUE),
specie_id = sample(1:4, size = n_rows, replace = TRUE),
citation_id = sample(1:8, size = n_rows, replace = TRUE),
treatment_id = sample(1:3, size = n_rows, replace = TRUE),
cultivar_id = NA,
date = NA,
time = NA
))
}
n_rows <- 40
trait_data <- list()
trait_data[["Amax"]] <- simulate_data(n_rows, 4.0, 20.0, 0.13, 4.0)
trait_data[["SLA"]] <- simulate_data(n_rows, 3.0, 15.0, 0.06, 1.7)
trait_data[["leaf_turnover_rate"]] <- simulate_data(n_rows, 0.1, 0.5, 0.01, 0.05)Run the meta analysis, performing prior and posterior checks and summarizing the results.
ma_result <- meta_analysis_standalone(trait_data, priors, iterations = 1000)## Each meta-analysis will be run with:
## 1000 total iterations,
## 4 chains,
## a burnin of 500 samples,
## ,
## thus the total number of samples will be 2000
## ################################################
## ------------------------------------------------
## starting meta-analysis for:
##
## Amax
##
## ------------------------------------------------
## prior for Amax
## (using R parameterization):
## unif(0, 40)
## data max: 19.7234384641051
## data min: 4.49936119094491
## mean: 11.1
## n: 40
## stem plot of data points
##
## The decimal point is at the |
##
## 4 | 57047
## 6 | 267123356
## 8 | 17288
## 10 | 6847
## 12 | 225712
## 14 | 36
## 16 | 0483
## 18 | 04457
##
## stem plot of obs.prec:
##
## The decimal point is 1 digit(s) to the left of the |
##
## 0 | 00000000000000000000000000000000011267
## 1 | 4
## 2 |
## 3 |
## 4 |
## 5 |
## 6 | 3
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 80
## Unobserved stochastic nodes: 14
## Total graph size: 334
##
## Initializing model
##
##
## Iterations = 2002:12000
## Thinning interval = 2
## Number of chains = 4
## Sample size per chain = 5000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## beta.ghs[2] -1.4420 0.47136 0.0033330 0.003826
## beta.o 11.8257 1.14627 0.0081053 0.067789
## beta.site[10] -2.8659 1.35183 0.0095589 0.063092
## beta.site[1] 3.1129 1.22795 0.0086829 0.066724
## beta.site[2] 6.5100 1.27668 0.0090275 0.067806
## beta.site[3] 2.0148 1.19730 0.0084662 0.067373
## beta.site[4] -1.2421 1.24042 0.0087711 0.065162
## beta.site[5] -1.6030 1.24873 0.0088299 0.065828
## beta.site[6] -0.4874 1.18008 0.0083444 0.067667
## beta.site[7] -1.3037 1.22139 0.0086365 0.065208
## beta.site[8] -1.9042 1.23062 0.0087018 0.068134
## beta.site[9] -2.0757 1.20788 0.0085410 0.065506
## sd.site 3.2815 0.93140 0.0065860 0.022511
## sd.y 6.7515 0.07951 0.0005622 0.000585
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## beta.ghs[2] -2.3644 -1.756 -1.4453 -1.1219 -0.5168
## beta.o 9.6659 11.115 11.7845 12.4724 14.2770
## beta.site[10] -5.6968 -3.699 -2.8192 -1.9781 -0.3225
## beta.site[1] 0.5441 2.373 3.1510 3.8892 5.4946
## beta.site[2] 3.8897 5.716 6.5221 7.3173 9.0094
## beta.site[3] -0.5210 1.303 2.0601 2.7797 4.2818
## beta.site[4] -3.8537 -1.984 -1.2026 -0.4387 1.1217
## beta.site[5] -4.2917 -2.337 -1.5561 -0.8056 0.7598
## beta.site[6] -3.0183 -1.161 -0.4479 0.2548 1.7447
## beta.site[7] -3.9233 -2.019 -1.2599 -0.5334 1.0110
## beta.site[8] -4.5498 -2.619 -1.8534 -1.1217 0.3796
## beta.site[9] -4.6903 -2.775 -2.0306 -1.3053 0.2113
## sd.site 2.0109 2.633 3.1048 3.7290 5.5556
## sd.y 6.5960 6.697 6.7520 6.8050 6.9055## Warning in FUN(X[[i]], ...): start value not changed
## Warning in FUN(X[[i]], ...): start value not changed
## Warning in FUN(X[[i]], ...): start value not changed
## Warning in FUN(X[[i]], ...): start value not changed## ################################################
## ------------------------------------------------
## starting meta-analysis for:
##
## SLA
##
## ------------------------------------------------
## prior for SLA
## (using R parameterization):
## lnorm(1.89, 0.61)
## data max: 14.8496243376285
## data min: 3.24236821196973
## mean: 9.44
## n: 40
## stem plot of data points
##
## The decimal point is at the |
##
## 2 | 246
## 4 | 33833
## 6 | 257
## 8 | 13345589023
## 10 | 45780239
## 12 | 12456
## 14 | 13788
##
## stem plot of obs.prec:
##
## The decimal point is at the |
##
## 0 | 000000000000000000000000000000000011122
## 1 |
## 2 |
## 3 |
## 4 |
## 5 |
## 6 |
## 7 |
## 8 |
## 9 | 8
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 80
## Unobserved stochastic nodes: 14
## Total graph size: 327
##
## Initializing model
##
##
## Iterations = 1002:2000
## Thinning interval = 2
## Number of chains = 4
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## beta.ghs[2] -2.1145 0.31930 0.0071397 0.007217
## beta.o 9.2904 0.72892 0.0162991 0.158876
## beta.site[10] -3.0759 0.76672 0.0171445 0.164020
## beta.site[1] 3.0925 0.74808 0.0167276 0.165726
## beta.site[2] 0.2960 0.74860 0.0167392 0.151693
## beta.site[3] -0.5173 0.75107 0.0167944 0.172714
## beta.site[4] -2.8397 0.76036 0.0170021 0.157784
## beta.site[5] -2.4925 0.75345 0.0168476 0.149637
## beta.site[6] 1.0989 0.73808 0.0165039 0.162774
## beta.site[7] 4.9682 1.02993 0.0230300 0.120173
## beta.site[8] 0.5767 0.80166 0.0179257 0.141919
## beta.site[9] 1.1812 0.76368 0.0170764 0.142086
## sd.site 2.8209 0.75496 0.0168815 0.019535
## sd.y 3.3058 0.03076 0.0006878 0.000788
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## beta.ghs[2] -2.7496 -2.326004 -2.1126 -1.90474 -1.4850
## beta.o 7.9663 8.720092 9.3383 9.83792 10.5471
## beta.site[10] -4.4915 -3.646018 -3.1092 -2.50264 -1.6687
## beta.site[1] 1.7966 2.533707 3.0524 3.66110 4.4651
## beta.site[2] -1.0476 -0.274346 0.2607 0.89234 1.6683
## beta.site[3] -1.8639 -1.090763 -0.5518 0.07715 0.8494
## beta.site[4] -4.2512 -3.403460 -2.8673 -2.27177 -1.4259
## beta.site[5] -3.8461 -3.069205 -2.5413 -1.92020 -1.0728
## beta.site[6] -0.1934 0.540058 1.0692 1.68641 2.4591
## beta.site[7] 2.9988 4.253963 4.9368 5.68109 7.0028
## beta.site[8] -0.9230 0.001037 0.5517 1.18345 2.0641
## beta.site[9] -0.2052 0.624986 1.1410 1.75005 2.6209
## sd.site 1.7698 2.305645 2.6838 3.17408 4.6230
## sd.y 3.2466 3.285053 3.3051 3.32605 3.3699## Warning in FUN(X[[i]], ...): start value not changed
## Warning in FUN(X[[i]], ...): start value not changed
## Warning in FUN(X[[i]], ...): start value not changed
## Warning in FUN(X[[i]], ...): start value not changed## ################################################
## ------------------------------------------------
## starting meta-analysis for:
##
## leaf_turnover_rate
##
## ------------------------------------------------
## prior for leaf_turnover_rate
## (using R parameterization):
## weibull(1.37, 1.43)
## data max: 0.494713499397039
## data min: 0.109222034458071
## mean: 0.264
## n: 40
## stem plot of data points
##
## The decimal point is 1 digit(s) to the left of the |
##
## 1 | 112223344
## 1 | 5578899
## 2 | 123
## 2 | 55689
## 3 | 000234
## 3 | 59
## 4 | 113
## 4 | 56799
##
## stem plot of obs.prec:
##
## The decimal point is 5 digit(s) to the right of the |
##
## 0 | 00000000000000000000000001111112234
## 0 | 56
## 1 | 3
## 1 | 5
## 2 |
## 2 |
## 3 |
## 3 |
## 4 |
## 4 | 5
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 80
## Unobserved stochastic nodes: 14
## Total graph size: 345
##
## Initializing model
##
##
## Iterations = 2002:12000
## Thinning interval = 2
## Number of chains = 4
## Sample size per chain = 5000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## beta.ghs[2] -0.02533 0.014808 1.047e-04 1.142e-04
## beta.o 0.61636 0.473198 3.346e-03 1.251e-01
## beta.site[10] -0.42505 0.473885 3.351e-03 1.238e-01
## beta.site[1] -0.37332 0.473397 3.347e-03 1.241e-01
## beta.site[2] -0.35250 0.473202 3.346e-03 1.213e-01
## beta.site[3] -0.41153 0.473525 3.348e-03 1.251e-01
## beta.site[4] -0.32480 0.473002 3.345e-03 1.258e-01
## beta.site[5] -0.41561 0.474450 3.355e-03 1.271e-01
## beta.site[6] -0.34330 0.473173 3.346e-03 1.250e-01
## beta.site[7] -0.28043 0.471975 3.337e-03 1.223e-01
## beta.site[8] -0.47762 0.474845 3.358e-03 1.268e-01
## beta.site[9] -0.20057 0.472311 3.340e-03 1.268e-01
## sd.site 0.44141 0.501753 3.548e-03 1.156e-01
## sd.y 0.16455 0.001919 1.357e-05 1.426e-05
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## beta.ghs[2] -0.05454 -0.03526 -0.025271 -0.015358 0.0035796
## beta.o 0.20732 0.25134 0.288883 1.002865 1.6855714
## beta.site[10] -1.49714 -0.81068 -0.099174 -0.058943 -0.0131513
## beta.site[1] -1.44359 -0.76129 -0.046640 -0.007973 0.0377670
## beta.site[2] -1.42030 -0.74014 -0.029628 0.012692 0.0626442
## beta.site[3] -1.48281 -0.79831 -0.085011 -0.046212 0.0002769
## beta.site[4] -1.39347 -0.71092 0.001051 0.040348 0.0861524
## beta.site[5] -1.48782 -0.80719 -0.090689 -0.049071 -0.0005336
## beta.site[6] -1.41447 -0.73086 -0.016927 0.021760 0.0675481
## beta.site[7] -1.35034 -0.66678 0.041786 0.082950 0.1336213
## beta.site[8] -1.54896 -0.86820 -0.151495 -0.110487 -0.0629440
## beta.site[9] -1.26741 -0.58683 0.124529 0.163466 0.2097895
## sd.site 0.05688 0.08023 0.114930 0.771905 1.6593078
## sd.y 0.16076 0.16327 0.164552 0.165842 0.1683330## Warning in FUN(X[[i]], ...): start value not changed
## Warning in FUN(X[[i]], ...): start value not changed
## Warning in FUN(X[[i]], ...): start value not changed
## Warning in FUN(X[[i]], ...): start value not changedThis returns three things:
names(ma_result)## [1] "trait.mcmc" "post.distns" "jagged.data"- The posterior distributions…
ma_result[["post.distns"]]## distn parama paramb n
## growth_resp_factor beta 2.630000 6.5200000 0
## leaf_turnover_rate weibull 1.370000 1.4300000 363
## root_respiration_rate unif 0.000000 100.0000000 NA
## root_turnover_rate unif 0.000000 10.0000000 NA
## Amax gamma 106.435422 9.0003222 NA
## leaf_respiration_rate_m2 weibull 2.000000 6.0000000 NA
## SLA norm 9.290374 0.7287377 455
## leafC norm 50.600000 1.3200000 291
## Vm_low_temp norm 0.000000 3.0000000 NA
## AmaxFrac unif 0.600000 0.9000000 NA
## psnTOpt unif 5.000000 40.0000000 NA
## stem_respiration_rate unif 0.000000 100.0000000 NA
## extinction_coefficient unif 0.380000 0.6200000 NA
## half_saturation_PAR unif 4.000000 27.0000000 NA
## dVPDSlope unif 0.010000 0.2500000 NA
## dVpdExp unif 1.000000 3.0000000 NA
## veg_respiration_Q10 unif 1.400000 2.6000000 NA
## fine_root_respiration_Q10 unif 1.400000 5.0000000 NA
## coarse_root_respiration_Q10 unif 1.400000 5.0000000 NA…(2) the full MCMC samples (as a list)…
tail(ma_result[["trait.mcmc"]][["Amax"]][[1]])## Markov Chain Monte Carlo (MCMC) output:
## Start = 11988
## End = 12000
## Thinning interval = 2
## beta.ghs[2] beta.o beta.site[10] beta.site[1] beta.site[2] beta.site[3]
## [1,] -1.1984728 12.72950 -2.704347 3.054231 6.224768 1.4579983
## [2,] -0.6427274 12.44421 -2.989694 1.848453 4.849856 1.4565529
## [3,] -1.5877333 12.67762 -3.618969 2.156694 6.339386 1.5205862
## [4,] -0.7963447 12.74144 -3.507818 2.562183 6.235358 0.7004895
## [5,] -1.2248751 12.71560 -3.335489 1.549860 6.733999 1.0221425
## [6,] -2.0792234 12.76844 -3.325667 2.753379 6.276371 0.9583105
## [7,] -2.9090424 12.78719 -5.119479 1.872165 6.252332 1.5218167
## beta.site[4] beta.site[5] beta.site[6] beta.site[7] beta.site[8]
## [1,] -1.805194 -2.273133 -2.0846222 -2.917668 -2.983224
## [2,] -3.063312 -2.885844 -1.1212161 -1.796047 -2.772800
## [3,] -2.586944 -2.743185 -1.2567432 -1.986648 -2.857963
## [4,] -2.631005 -2.834758 -0.9809788 -2.531084 -2.801400
## [5,] -2.922042 -3.162593 -1.5664136 -2.903866 -3.058427
## [6,] -1.608455 -1.730278 -1.0915844 -2.590668 -2.528091
## [7,] -2.172474 -1.756141 -1.0706952 -1.985520 -3.558391
## beta.site[9] sd.site sd.y
## [1,] -2.725720 3.705370 6.748434
## [2,] -2.679795 4.069632 6.742993
## [3,] -3.187307 3.929488 6.785414
## [4,] -2.718769 4.552949 6.760889
## [5,] -3.453120 2.803170 6.766290
## [6,] -3.336999 2.356853 6.718523
## [7,] -3.676124 3.534487 6.726614…(3) and the “JAGS-ified” trait data.
head(ma_result[["jagged.data"]][["Amax"]])## Y n site trt ghs obs.prec se cite greenhouse site_id
## 1 18.046454 74 1 control 1 0.001200041 3.3557231 6 1 1
## 2 18.480473 53 1 control 1 0.014015658 1.1602600 1 1 1
## 3 4.657720 39 1 control 1 0.001901748 3.6719033 4 1 1
## 4 13.236781 15 2 control 1 0.246123176 0.5204489 1 1 2
## 5 19.723438 92 2 control 1 0.003531625 1.7543603 4 1 2
## 6 9.803368 2 2 control 1 0.374649089 1.1552412 5 1 2
## treatment_id trt_name trt_num
## 1 2 t01 1
## 2 1 t19 1
## 3 2 t25 1
## 4 2 t02 1
## 5 3 t15 1
## 6 2 t36 1