diff --git a/.gitignore b/.gitignore index 4213df5fcdd0f67a102c62fa4b1bc866bc6681ef..d6523264e7523bd51f13d8be266b6da5f805d1f4 100644 --- a/.gitignore +++ b/.gitignore @@ -14,3 +14,4 @@ figs *.rds docs/metadata/sites .Rproj.user +GROWTH.MODEL.* diff --git a/R/analysis/fun.analysis.R b/R/analysis/JAGS.fun.R similarity index 93% rename from R/analysis/fun.analysis.R rename to R/analysis/JAGS.fun.R index 892cc18dd02148d241535fc171296fa11816d829..82ae982998f6550b631fa3014f768bdc42bdae3a 100644 --- a/R/analysis/fun.analysis.R +++ b/R/analysis/JAGS.fun.R @@ -1,7 +1,8 @@ ########################### ########################### ##### FUNCTIONS USED IN FORMATTING DATA FOR JAGS AND RUN ANALYSIS - +source("R/packages.R") +check_packages(c("rjags", "R2jags")) ###### ###### @@ -165,7 +166,8 @@ return(jags.data) ######################################### ### WRITE JAGS MODELS - +write.jags.model.HD <- function(file = "GROWTH.MODEL.LOGLIN.HD"){ + ### LOG LIN HD GROWTH.MODEL.LOGLIN.HD <- "################################################################################ @@ -233,7 +235,10 @@ sigma.alpha ~ dunif(0.000001,5) } # End of the jags model " -cat(GROWTH.MODEL.LOGLIN.HD , file = "GROWTH.MODEL.LOGLIN.HD", sep=" ", fill = FALSE, labels = NULL, append = FALSE) +cat(GROWTH.MODEL.LOGLIN.HD , file = file, sep=" ", fill = FALSE, labels = NULL, append = FALSE) +} + +write.jags.model.AD <- function(file = "GROWTH.MODEL.LOGLIN.AD"){ ### LOG LIN HD GROWTH.MODEL.LOGLIN.AD <- @@ -302,8 +307,8 @@ sigma.alpha ~ dunif(0.000001,5) } # End of the jags model " -cat(GROWTH.MODEL.LOGLIN.AD , file = "GROWTH.MODEL.LOGLIN.AD", sep=" ", fill = FALSE, labels = NULL, append = FALSE) - +cat(GROWTH.MODEL.LOGLIN.AD , file = file, sep=" ", fill = FALSE, labels = NULL, append = FALSE) +} ################################################ ########################### @@ -455,3 +460,22 @@ n.chain <- dim(jags.output$BUGSoutput$sims.array)[2] length.chain <- dim(jags.output$BUGSoutput$sims.array)[1] plot(x,y,col=rep(1:n.chain,each=length.chain),xlab=param.name.1,ylab=param.name.2) } + + +get.ecoregions <- function(site){ + sub(paste("./output/processed/",site,"/",sep=""),"",list.dirs(paste("./output/processed/",site,sep="")))[-1] +} + +run.jags.for.trait <- function(site, trait, A1, A2, ecoregions = get.ecoregions(site)){ + for (e in ecoregions) + try(fun.run.4models.jags(site=site,ecoregion.name=e,min.obs=10, + name.trait=trait, + A1=A1, + A2=A2, + n.chains=3, + n.iter=3000, + n.burnin=1000, + n.thin=1, + sample.size=2000, + na.action.traits="fill with global mean")) +} diff --git a/R/analysis/run.local.R b/R/analysis/run.local.R index a808bb9f807e3dc4c40f224fefa71f5b999d2de5..2c3da7f147e2bda61a0ce26c41d723550d7e84f0 100644 --- a/R/analysis/run.local.R +++ b/R/analysis/run.local.R @@ -1,20 +1,24 @@ ##### -## run jags code - -source("./R/analysis/fun.analysis.R") -site <- "NVS" - -ecoregion.code <- sub(paste("./output/processed/",site,"/",sep=""),"",list.dirs(paste("./output/processed/",site,sep="")))[-1] - -for (i in ecoregion.code[1]){ -fun.run.4models.jags(site=site,ecoregion.name=i,min.obs=10, - name.trait="SLA", - A1="MAT", - A2="MAP", - n.chains=3, - n.iter=3000, - n.burnin=1000, - n.thin=1, - sample.size=NA, - na.action.traits="fill with global mean") -} + +source("R/analysis/JAGS.fun.R") + +write.jags.model.HD() +write.jags.model.AD() + +traits <- c("SLA", "Leaf.N","Seed.mass","SLA","Wood.density","Max.height") + +for(this.trait in traits) + run.jags.for.trait("France", trait = this.trait, A1 = "sgdd", A2 = "WB.s") + +for(this.trait in traits) + run.jags.for.trait("NVS", trait = this.trait, A1 = "MAT", A2 = "MAP") + +for(this.trait in traits) + run.jags.for.trait("US", trait = this.trait, A1 = "MAT", A2 = "MAP") + +for(this.trait in traits) + run.jags.for.trait("Spain", trait = this.trait, A1 = "MAT", A2 = "PP") + +for(this.trait in traits) + run.jags.for.trait("Canada", trait = this.trait, A1 = "MAT", A2 = "MAP") + diff --git a/readme.md b/readme.md index b9e7b4e7020734ad225b23161998f7d1ad38bf55..3b09209bfcb621b759333c3beba0d794678124d1 100644 --- a/readme.md +++ b/readme.md @@ -33,6 +33,7 @@ We are making use of the following tools: - [*markdown*](http://daringfireball.net/projects/markdown/) - The [*knitr*](http://cran.r-project.org/web/packages/knitr/index.html) package in R, for generating documents with R and markdown. For a n intro see [here](http://nicercode.github.io/guides/reports/) - version control via git (see next section) +- JAGS - Just another Gibbs Sampler. Install [here](http://sourceforge.net/projects/mcmc-jags/files/) ### Git