PEcAn.allometry Vignette
install package from Github
Only needs to be done the first time
Define species groups
library(PEcAn.allometry)
## define the Plant Functional Types you want to fit based on USDA acronyms and USFS species codes (spcd)
## can involve one or many species
## multiple PFTs can be fit at once by putting them in a list
## Note that the name used in the PFT list defines the name that the code will match against when making predictions for new individuals
pfts = list(FAGR = data.frame(spcd=531,acronym="FAGR"))
## example of a PFT with multiple species (LH = late hardwood)
## note that if you're just using Jenkins the acronym column is optional
pfts = list(LH = data.frame(spcd = c(531,318),acronym=c("FAGR","ACSA3")))Run the Bayesian allometry model
Note the side effects of this function are that it will create two files in your working directory for every PFT and component pair, a pdf file of diagnostics and a RData file saving the full output. These files will be named based on the PFT and the component.
The return from this function will be a summary table of statistics. The function will actually run two variants of the model, a “global” modelthat fits a single equation to all equations and a ‘hierarchical’ model that accounts for the variability among equations. This function also print out DIC statistics on which fit was better (lowest score wins): DIC is the hierarchical model, DICg is the global model.
## [1] "writing output to"
## [2] "/tmp/RtmpOsxXzQ/Rbuild115576e7ec7b/PEcAn.allometry/vignettes"
## $spcd
## [1] 531 318
##
## $acronym
## [1] "FAGR" "ACSA3"
##
## [[3]]
## [1] 6## Warning in nu(allom$parm$Component.ID): NAs introduced by coercion## [1] "Dropping allom rows: "
## | | | 0%
## [1] "Dropping allom rows: "
## | | | 0%
## [1] "Dropping allom rows: "
## | | | 0%
## [1] "saving MCMC output to"
## [2] "/tmp/RtmpOsxXzQ/Rbuild115576e7ec7b/PEcAn.allometry/vignettes/Allom.LH.6.Rdata"
## [1] "saving diagnostic graphs to"
## [2] "/tmp/RtmpOsxXzQ/Rbuild115576e7ec7b/PEcAn.allometry/vignettes/Allom.LH.6.MCMC.pdf"## [1] "DIC" "400.66525424069" "pD" "239.296347468632"
## [1] "DICg" "474.35233564325" "pDg" "234.806964646513"If you want to run with a response variable other than the default (e.g. components = 6; stem biomass), look up the relevant component IDs in data(allom.components). The default component is 3 (total aboveground biomass). Note that if you specify multiple PFTs (as a list) and multiple components (as a vector) then AllomAve will generate allometries for all PFT x component combinations
## [1] "writing output to"
## [2] "/tmp/RtmpOsxXzQ/Rbuild115576e7ec7b/PEcAn.allometry/vignettes"
## $spcd
## [1] 531 318
##
## $acronym
## [1] "FAGR" "ACSA3"
##
## [[3]]
## [1] 3## Warning in nu(allom$parm$Component.ID): NAs introduced by coercion## [1] "Dropping allom rows: "
## | | | 0%
## [1] "Dropping allom rows: "
## | | | 0%
## [1] "Dropping allom rows: "
## | | | 0%
## [1] "saving MCMC output to"
## [2] "/tmp/RtmpOsxXzQ/Rbuild115576e7ec7b/PEcAn.allometry/vignettes/Allom.LH.3.Rdata"
## [1] "saving diagnostic graphs to"
## [2] "/tmp/RtmpOsxXzQ/Rbuild115576e7ec7b/PEcAn.allometry/vignettes/Allom.LH.3.MCMC.pdf"## [1] "DIC" "407.754500784945" "pD" "186.541641869388"
## [1] "DICg" "462.505292742636" "pDg" "215.992457641081"
## $spcd
## [1] 531 318
##
## $acronym
## [1] "FAGR" "ACSA3"
##
## [[3]]
## [1] 6## Warning in nu(allom$parm$Component.ID): NAs introduced by coercion## [1] "Dropping allom rows: "
## | | | 0%
## [1] "Dropping allom rows: "
## | | | 0%
## [1] "Dropping allom rows: "
## | | | 0%
## [1] "saving MCMC output to"
## [2] "/tmp/RtmpOsxXzQ/Rbuild115576e7ec7b/PEcAn.allometry/vignettes/Allom.LH.6.Rdata"
## [1] "saving diagnostic graphs to"
## [2] "/tmp/RtmpOsxXzQ/Rbuild115576e7ec7b/PEcAn.allometry/vignettes/Allom.LH.6.MCMC.pdf"## [1] "DIC" "416.958962387703" "pD" "247.702517225291"
## [1] "DICg" "478.628908367374" "pDg" "233.505057022347"Predict for individual trees
allom.fit = load.allom(getwd())
dbh = 1:50
pred = allom.predict(allom.fit,dbh = dbh,pft = "LH",component = 3,use = "Bg",interval = "prediction")
conf = allom.predict(allom.fit,dbh = dbh,pft = "LH",component = 3,use = "Bg",interval = "confidence")
PI = apply(pred,2,quantile,c(0.025,0.5,0.975),na.rm=TRUE)
CI = apply(conf,2,quantile,c(0.025,0.5,0.975),na.rm=TRUE)
plot(dbh,CI[2,],type='l',lwd=3,ylim=range(PI),ylab="Biomass (kg)")
lines(dbh,CI[1,],lty=2,col="blue")
lines(dbh,CI[3,],lty=2,col="blue")
lines(dbh,PI[1,],lty=3,col="red")
lines(dbh,PI[3,],lty=3,col="red")
