diff --git a/R/analysis/JAGS.run.R b/R/analysis/JAGS.run.R index 2902c94e45497c0962028632e36c87a9cd546c1b..5e9292054779e00ea369c27c7a97360bb5cf7ab5 100644 --- a/R/analysis/JAGS.run.R +++ b/R/analysis/JAGS.run.R @@ -110,7 +110,8 @@ load.and.prepare.data.for.lmer <- function(trait, set, ecoregion, } -fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,BATOT,min.obs,min.BA.G=-50,max.BA.G=150){ +fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,perc.genus,BATOT,min.obs, + min.perc.genus=0.90,min.BA.G=-50,max.BA.G=150){ ## remove tree with NA data.tree <- subset(data.tree,subset=(!is.na(data.tree[["BA.G"]])) & (!is.na(data.tree[["D"]])) ) @@ -123,6 +124,8 @@ data.tree <- data.tree[(!is.na(data.tree[[abs.CWM.tntf]]) & # select species with minimum obs data.tree <- subset(data.tree,subset=data.tree[["sp"]] %in% names(table(factor(data.tree[["sp"]])))[table(factor(data.tree[["sp"]]))>min.obs]) +# select obs abov min perc.genusspecies with minimum obs +data.tree <- subset(data.tree,subset=data.tree[[perc.genus]] > min.perc.genus & !is.na(data.tree[[perc.genus]])) return(data.tree) } @@ -190,8 +193,9 @@ CWM.tn <- paste(trait,"CWM",type.filling,"log",sep=".") abs.CWM.tntf <- paste(trait,"abs.CWM",type.filling,"log",sep=".") tf <- paste(trait,"focal",sep=".") BATOT <- "BATOT.log" +perc.genus <- paste(trait,"perc.genus",sep=".") -data.tree <- fun.select.data.for.analysis(data.tree,abs.CWM.tntf,BATOT,min.obs) +data.tree <- fun.select.data.for.analysis(data.tree,abs.CWM.tntf,perc.genus,BATOT,min.obs) list.data <- fun.get.the.variables.for.lmer(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf) #= DATA LIST FOR JAGS diff --git a/R/analysis/jags.output.R b/R/analysis/jags.output.R index 9bd0ff9f8b8ab51c12ff84f967ccf698aca2dd79..d6e694636dadd5657846696bbe36ecfe5fb2606a 100644 --- a/R/analysis/jags.output.R +++ b/R/analysis/jags.output.R @@ -1,9 +1,9 @@ source("R/analysis/jags.output-fun.R") # summarises output for all jags models -## files <- list.files("output/jags", recursive=TRUE, full.names =TRUE, pattern = "rds") -## out <- rbind.fill(lapply(files, summarise.jags.output.table)) -## write.csv(out, file = "output/jags/results.csv", quote = FALSE, row.names = FALSE) +files <- list.files("output/lmer", recursive=TRUE, full.names =TRUE, pattern = "rds") +out <- rbind.fill(lapply(files, summarise.jags.output.table)) +write.csv(out, file = "output/jags/results.csv", quote = FALSE, row.names = FALSE) out <- read.csv(file = "output/results.csv", stringsAsFactors = FALSE) diff --git a/R/analysis/lmer.output-fun.R b/R/analysis/lmer.output-fun.R new file mode 100644 index 0000000000000000000000000000000000000000..cf9acf56f5d010fc2141465a21bfc8d50064beed --- /dev/null +++ b/R/analysis/lmer.output-fun.R @@ -0,0 +1,98 @@ +#### function to analyse lmer output + + +library(lme4) + + +read.lmer.output <- function(file.name){ + tryCatch(readRDS(file.name), error=function(cond)return(NULL)) # Choose a return value in case of error +} + + + +summarise.lmer.output <- function(x){ + + list( nobs = slot(x,'dims')['n'],# nobs(x) + R2m =Rsquared.glmm.lmer(x)$Marginal, + R2c =Rsquared.glmm.lmer(x)$Conditional, + AIC = AIC(x), + deviance = deviance(x), + fixed.coeff.E=fixef(x), + fixed.coeff.Std.Error=sqrt(diag(vcov(x)))) + +} + +## summarise.lmer.output <- function(x){ +## list( nobs = nobs(x),#slot(x,'dims')['n'] +## R2m =Rsquared.glmm.lmer(x)$Marginal, +## R2c =Rsquared.glmm.lmer(x)$Conditional, +## AIC = AIC(x), +## deviance = deviance(x), +## fixed.coeff.E=fixef(x), +## fixed.coeff.Std.Error=sqrt(diag(vcov(x)))) +## } + +summarise.lmer.output.list <- function(f ){ + output.lmer <- read.lmer.output(f) + if(!is.null(output.lmer)){ + res <- list(files.details=files.details(f), + lmer.summary=summarise.lmer.output( output.lmer)) + }else{ + res <- NULL + } + return(res) +} + + + +files.details <- function(x){ + s <- data.frame(t(strsplit(x, "/", fixed=TRUE)[[1]]), stringsAsFactors= FALSE) + names(s) <- c("d1", "d2", "model", "trait", "set", "ecocode", "filling","file" ) + s[-(1:2)] +} + + + +#### R squred functions + +# Function rsquared.glmm requires models to be input as a list (can include fixed- +# effects only models,but not a good idea to mix models of class "mer" with models +# of class "lme") FROM http://jslefche.wordpress.com/2013/03/13/r2-for-linear-mixed-effects-models/ + + +Rsquared.glmm.lmer <- function(i){ + # Get variance of fixed effects by multiplying coefficients by design matrix + VarF <- var(as.vector(fixef(i) %*% t(slot(i,'X')))) + # Get variance of random effects by extracting variance components + VarRand <- colSums(do.call(rbind,lapply(VarCorr(i),function(j) j[1]))) + # Get residual variance + VarResid <- attr(VarCorr(i),"sc")^2 + # Calculate marginal R-squared (fixed effects/total variance) + Rm <- VarF/(VarF+VarRand+VarResid) + # Calculate conditional R-squared (fixed effects+random effects/total variance) + Rc <- (VarF+VarRand)/(VarF+VarRand+VarResid) + # Bind R^2s into a matrix and return with AIC values + Rsquared.mat <- data.frame(Class=class(i),Family="Gaussian",Marginal=Rm, + Conditional=Rc,AIC=AIC(i)) +return(Rsquared.mat) +} + +## Rsquared.glmm.lmer <- function(i){ +## # Get variance of fixed effects by multiplying coefficients by design matrix +## VarF <- var(as.vector(fixef(i) %*% t(i@pp$X))) +## # Get variance of random effects by extracting variance components +## VarRand <- colSums(do.call(rbind,lapply(VarCorr(i),function(j) j[1]))) +## # Get residual variance +## VarResid <- attr(VarCorr(i),"sc")^2 +## # Calculate marginal R-squared (fixed effects/total variance) +## Rm <- VarF/(VarF+VarRand+VarResid) +## # Calculate conditional R-squared (fixed effects+random effects/total variance) +## Rc <- (VarF+VarRand)/(VarF+VarRand+VarResid) +## # Bind R^2s into a matrix and return with AIC values +## Rsquared.mat <- data.frame(Class=class(i),Family="Gaussian",Marginal=Rm, +## Conditional=Rc,AIC=AIC(update(i,REML=F))) +## return(Rsquared.mat) +## } + + + diff --git a/R/analysis/lmer.output.R b/R/analysis/lmer.output.R new file mode 100755 index 0000000000000000000000000000000000000000..15da452f722afd6118bcb25fc1444e5cf6792be3 --- /dev/null +++ b/R/analysis/lmer.output.R @@ -0,0 +1,8 @@ +#!/usr/bin/env Rscript + +source("R/analysis/lmer.output-fun.R") + +files <- list.files("output/lmer", recursive=TRUE, full.names =TRUE, pattern = "rds") +out <- lapply(files, summarise.lmer.output.list) +names(out) <- lapply(lapply(files,files.details),function(x) paste(as.vector(x[-6]),collapse="_")) +saveRDS(out,file='output/lmer/list.out.rds') diff --git a/R/analysis/run.local.R b/R/analysis/run.local.R index 4738b74e77e0cba545c2542dcb479e4720491e97..7f2cef85ae18862ba93d579e5a714ec65517a960 100644 --- a/R/analysis/run.local.R +++ b/R/analysis/run.local.R @@ -1,6 +1,6 @@ ##### source("R/analysis/lmer.run.R") -sets.inv <- c("US","Sweden","France","NSW","NVS","Canada","Swiss") #"Spain", +sets.inv <- c("US","Sweden","France","NSW","NVS","Canada","Spain","Swiss") # sets.trop <- c("BCI","Fushan","Paracou","Mbaiki","Luquillo") #========== # Test lmer diff --git a/R/analysis/theo.compet.model/theo.model.R b/R/analysis/theo.compet.model/theo.model.R new file mode 100644 index 0000000000000000000000000000000000000000..6018ccd2f0cb1349a78845cdb70ca385582d8dd4 --- /dev/null +++ b/R/analysis/theo.compet.model/theo.model.R @@ -0,0 +1,270 @@ +########################### +########################### +### explore the link between competitive effect/response and abs trait distance + +### G. Kunstler 21/09/2013 +setwd("../..") + +### explore models of link between traits and compet +## traits range between 0 and 1 + +min.t <- 0 +max.t <- 1 +a <- seq(min.t,max.t,length=100) +b <- seq(min.t,max.t,length=100) + +## fit to response/effect model +a.r <- runif(1000,min.t,max.t) +b.r <- runif(1000,min.t,max.t) +abs.r <- 1-abs(a.r-b.r) +hier.r <- 0.5*(a.r-b.r) +lm.abs.dist <- lm(abs.r~a.r*b.r) +lm.hier.dist <- lm(hier.r~a.r*b.r) +lm.abs.dist.nointer <- lm(abs.r~a.r+b.r) + +logistic <- function(x) 1/(1+exp(-(x))) + +mat.c1.dist.pred <- matrix(NA,100,100) +mat.c2.dist.pred <- matrix(NA,100,100) +mat.c3.dist.pred <- matrix(NA,100,100) +mat.c4.dist.pred <- matrix(NA,100,100) +mat.c5.dist.pred <- matrix(NA,100,100) +mat.c6.dist.pred<-matrix(NA,100,100) +mat.c7.dist.pred<-matrix(NA,100,100) +for (i in 1:100){ + for (j in 1:100){ + mat.c1.dist.pred[i,j] <- 0.5+0.5*a[i]-0.5*b[j] + mat.c2.dist.pred[i,j] <- 1-abs(a[i]-b[j]) + mat.c3.dist.pred[i,j] <- 0.65+0.04*a[i] -0.04*b[j] #1.272-1.206*a[i]-1.206*b[j]+2.4604*a[i]*b[j] + mat.c4.dist.pred[i,j] <-(0.769-1.568*b[j])*(0.769-1.568*a[i]) # (1.127-1.07*b[j])*(1.127-1.07*a[i]) # + mat.c5.dist.pred[i,j] <- (1-1*b[j])*(1*a[i]) + mat.c6.dist.pred[i,j] <- 0.5-0.5*b[j]+0.5*a[i]-0*a[i]*b[j] + mat.c7.dist.pred[i,j] <- (0.5+1.5*b[j])*(-0.5+1.5*a[i]) + + } + } + +png(file="./figs/compet_model.png",width=600,height=900) +m <- matrix(c(1:2,7,3,4,7,5,6,7), 3, 3,byrow=TRUE) +nclass <- 12 +layout(m, widths=c(3, 3,0.5),heights=c(1, 1,1)) +zrange <- range(mat.c1.dist.pred,mat.c2.dist.pred,mat.c3.dist.pred,mat.c4.dist.pred,mat.c5.dist.pred) +library(RColorBrewer) +breakss <- c(seq(zrange[1],0,length=4),seq(0,zrange[2],length=nclass+1+1-4)[-1]) +col.pal <- c(rev(brewer.pal(4-1,"Blues")),brewer.pal(nclass-4+1,"Reds")) +image(x=a,y=b,z=mat.c1.dist.pred,zlim=zrange,xlab="tn",ylab="tf",main="Hierarchical distance",col=col.pal,breaks=breakss) +contour(a,b,mat.c1.dist.pred,levels=0,add=TRUE,col="red") +image(a,b,mat.c2.dist.pred,zlim=zrange,xlab="tn",ylab="tf",main="Absolute distance",col=col.pal,breaks=breakss) +contour(a,b,mat.c2.dist.pred,levels=0,add=TRUE,col="red") +image(a,b,mat.c6.dist.pred,zlim=zrange,xlab="tn",ylab="tf",main="linear with linear increase of both and no interaction",col=col.pal,breaks=breakss) +contour(a,b,mat.c6.dist.pred,levels=0,add=TRUE,col="red") +image(a,b,mat.c3.dist.pred,zlim=zrange,xlab="tn",ylab="tf",main="linear with best fit to abs dist",col=col.pal,breaks=breakss) +contour(a,b,mat.c3.dist.pred,levels=0,add=TRUE,col="red") +image(a,b,mat.c5.dist.pred,zlim=zrange,xlab="tn",ylab="tf",main="multiplicative with linear increase of both",col=col.pal,breaks=breakss) +contour(a,b,mat.c5.dist.pred,levels=0,add=TRUE,col="red") +image(a,b,mat.c4.dist.pred,zlim=zrange,xlab="tn",ylab="tf",main="multiplicative with close to best fit to abs dist",col=col.pal,breaks=breakss) +contour(a,b,mat.c4.dist.pred,levels=0,add=TRUE,col="red") +par(mar=c(0,0,0,0)) +image(x=1,y=breakss,matrix(breakss[-1],1,nclass),col=col.pal,breaks=breakss) +dev.off() + + +png(file="./figs/compet_facill_abs.png",width=200,height=600) +par(mfrow=c(3,1)) +image(a,b,mat.c4.dist.pred,xlab="tn",ylab="tf",main="multiplicative with close best fit to abs dist",col=col.pal,breaks=breakss) +plot(a,-0.769+1.568*a,type="l",xlab="tn",ylab="competitive effect") +abline(h=0,lty=2) +plot(b,0.769-1.568*b,type="l",xlab="tf",ylab="competitive response") +abline(h=0,lty=2) +dev.off() + +#### 4 different multiplicative +png(file="./figs/multi4models.png",width=1200,height=900) +mat.c7.dist.pred<-matrix(NA,100,100) +for (i in 1:100){ + for (j in 1:100){ + mat.c7.dist.pred[i,j] <- (0+1*b[j])*(0+1*a[i]) + } + } + +par(mfcol=c(3,4)) +image(a,b,mat.c7.dist.pred,xlab="tn",ylab="tf",main="(0+b)(0+a)",col=col.pal,breaks=breakss) +plot(a,0+1*a,type="l",xlab="tn",ylab="competitive effect") +abline(h=0,lty=2) +plot(b,0+1*b,type="l",xlab="tf",ylab="competitive response") +abline(h=0,lty=2) + +mat.c7.dist.pred<-matrix(NA,100,100) +for (i in 1:100){ + for (j in 1:100){ + mat.c7.dist.pred[i,j] <- (1-1*b[j])*(0+1*a[i]) + } + } + +image(a,b,mat.c7.dist.pred,xlab="tn",ylab="tf",main="(1-b)(0+a)",col=col.pal,breaks=breakss) +plot(a,0+1*a,type="l",xlab="tn",ylab="competitive effect") +abline(h=0,lty=2) +plot(b,1-1*b,type="l",xlab="tf",ylab="competitive response") +abline(h=0,lty=2) + +mat.c7.dist.pred<-matrix(NA,100,100) +for (i in 1:100){ + for (j in 1:100){ + mat.c7.dist.pred[i,j] <- (1-1*b[j])*(1-1*a[i]) + } + } + +image(a,b,mat.c7.dist.pred,xlab="tn",ylab="tf",main="(1-b)(1-a)",col=col.pal,breaks=breakss) +plot(a,1-a,type="l",xlab="tn",ylab="competitive effect") +abline(h=0,lty=2) +plot(b,1-1*b,type="l",xlab="tf",ylab="competitive response") +abline(h=0,lty=2) + +mat.c7.dist.pred<-matrix(NA,100,100) +for (i in 1:100){ + for (j in 1:100){ + mat.c7.dist.pred[i,j] <- (0+1*b[j])*(1-1*a[i]) + } + } + +image(a,b,mat.c7.dist.pred,xlab="tn",ylab="tf",main="(0+b)(1-a)",col=col.pal,breaks=breakss) +plot(a,1-a,type="l",xlab="tn",ylab="competitive effect") +abline(h=0,lty=2) +plot(b,0+1*b,type="l",xlab="tf",ylab="competitive response") +abline(h=0,lty=2) +dev.off() + +## 4 different additive +png(file="./figs/additiv4models.png",width=1200,height=900) +mat.c7.dist.pred<-matrix(NA,100,100) +for (i in 1:100){ + for (j in 1:100){ + mat.c7.dist.pred[i,j] <- 0.5*(1*b[j]+1*a[i]) + } + } + +par(mfcol=c(3,4)) +image(a,b,mat.c7.dist.pred,xlab="tn",ylab="tf",main="(0.5*b+0.5*a)",col=col.pal,breaks=breakss) +plot(a,0+0.5*a,type="l",xlab="tn",ylab="competitive effect") +abline(h=0,lty=2) +plot(b,0+0.5*b,type="l",xlab="tf",ylab="competitive response") +abline(h=0,lty=2) + +mat.c7.dist.pred<-matrix(NA,100,100) +for (i in 1:100){ + for (j in 1:100){ + mat.c7.dist.pred[i,j] <- 0.5*(1-1*b[j]+1*a[i]) + } + } + +image(a,b,mat.c7.dist.pred,xlab="tn",ylab="tf",main="(0.5-0.5*b+0.5a)",col=col.pal,breaks=breakss) +plot(a,0+0.5*a,type="l",xlab="tn",ylab="competitive effect") +abline(h=0,lty=2) +plot(b,0.5-0.5*b,type="l",xlab="tf",ylab="competitive response") +abline(h=0,lty=2) + +mat.c7.dist.pred<-matrix(NA,100,100) +for (i in 1:100){ + for (j in 1:100){ + mat.c7.dist.pred[i,j] <- 0.5*(2-1*b[j]-1*a[i]) + } + } + + +image(a,b,mat.c7.dist.pred,xlab="tn",ylab="tf",main="(1-0.5*b-0.5a)",col=col.pal,breaks=breakss) +plot(a,0.5-0.5*a,type="l",xlab="tn",ylab="competitive effect") +abline(h=0,lty=2) +plot(b,0.5-0.5*b,type="l",xlab="tf",ylab="competitive response") +abline(h=0,lty=2) + +mat.c7.dist.pred<-matrix(NA,100,100) +for (i in 1:100){ + for (j in 1:100){ + mat.c7.dist.pred[i,j] <- 0.5*(1+1*b[j]-1*a[i]) + } + } + +image(a,b,mat.c7.dist.pred,xlab="tn",ylab="tf",main="(0.5+0.5*b-0.5a)",col=col.pal,breaks=breakss) +plot(a,0.5-0.5*a,type="l",xlab="tn",ylab="competitive effect") +abline(h=0,lty=2) +plot(b,0.5*b,type="l",xlab="tf",ylab="competitive response") +abline(h=0,lty=2) +dev.off() + + +plot(as.vector(mat.c2.dist.pred),as.vector(mat.c4.dist.pred),xlab="abs dist",ylab="multiplicative best fit") +plot(as.vector(mat.c2.dist.pred),as.vector(mat.c3.dist.pred),xlab="abs dist",ylab="multiplicative best fit") +### what is the meaning of negative competitive effect or response ? facilitation ? + + +## TRY TO UNDERSTAND HOW the multiplicative model reduce the space in comparison with full model +max.r <- 20 +a <- runif(100000,-max.r,max.r) +b <- runif(100000,-max.r,max.r) +c <- runif(100000,-max.r,max.r) +d <- runif(100000,-max.r,max.r) + +a.2 <- a*c +b.2 <- c*b +c.2 <- a*d +d.2 <- b*d + +mat.param <- cbind(a.2,b.2,c.2,d.2) +pairs(mat.param) +library(scatterplot3d) +png(file="./figs/param.space.png",width=20,height=20) +par(mfrow=c(2,2)) +scatterplot3d(c.2,b.2,d.2,highlight.3d=TRUE,cex.symbols=0.1) + +a.2.s <- b.2.s <- seq(-max.r,max.r,length=100) +mat.d <- outer(a.2.s,b.2.s,FUN=function(x,y,max.r) {x*y/(max.r*max.r)},max.r=max.r) +mat.d2 <- outer(a.2.s,b.2.s,FUN=function(x,y,max.r) {x*y/(-max.r*max.r)},max.r=max.r) +persp(a.2.s,b.2.s,mat.d, theta = 10, phi = 15, expand = 0.99, col = "lightblue") +persp(a.2.s,b.2.s,mat.d2, theta = 25, phi = 15, expand = 0.99, col = "lightblue") +scatterplot3d(c(c.2,c.2),c(b.2,b.2),c(d.2.a,d.2.b),color=c(rep("red",length(b.2)),rep("blue",length(b.2))),cex.symbols=0.1) +dev.off() + +x11() +par(mfrow=c(3,3)) + ss <- 1 +for (lim in c(-300,-200,-100,-50,0,50,100,200,300)){ + plot(c.2[b.2> (lim-ss) & b.2<(lim+ss)],d.2[b.2> (lim-ss) & b.2<(lim+ss)],xlab="c.2",ylab="d.2",main=paste("slice for b.2 + or - ",ss,"of",lim)) + lines((-400):400,(-400):400*lim/(max.r*max.r),col="red") + lines((-400):400,(-400):400*lim/(-max.r*max.r),col="red") +} + +#### DAVID CODE + +# response and effect coefficients +r0 = 1; r1 = .5; e0 = 1; e1 = 0.5 +#traits +t1 = 1:100 +t2 = 1:100 +#work out the response/effect model +l = matrix(data = 0, nrow = 100, ncol = 100) +for (i in 1:100) l[i,] = ((r0+ r1*t1[i])*(e0 + e1*t2)) +filled.contour(t1,t2, l, xlab = "trait1", ylab = "trait2", main = "hierarchical") +#work out absolute model +l1 = matrix(data = 0, nrow = 100, ncol = 100) +for (i in 1:100) l1[i,] = abs(t1[i]-t2) +filled.contour(t1,t2, l1, xlab = "trait1", ylab = "trait2", main = "absolute") + +### JOE + +# Joe's fiddling about with David's script +# response and effect coefficients +r0 = 1; r1 = .5; e0 = 1; e1 = 0.5 +#traits +l2 = matrix(data=0, nrow=10000, ncol=5) +l2[,1] = rep(1:100, times=100) +l2[,2] = rep(1:100, each=100) +l2[,3] = l2[,1]*l2[,2] +#work out absolute model +l2[,4] = abs(l2[,1]-l2[,2]) +#work out the response/effect model +l2[,5] = (r0+ r1*l2[,1])*(e0 + e1*l2[,2]) +# evaluate relationship between the absolute model & the components of the response/effect model +model <- lm(l2[,4]~1+l2[,1]+l2[,2]+l2[,3]) +summary(model) +plot(l2[,4],model$fitted.values) + diff --git a/R/format.data/BCI.R b/R/format.data/BCI.R index f6d4e96a6bf1d0c0a16038e5b8ce6b9ca8c0adea..3e6353bc1a6025b25f92041f612942594122cc9b 100644 --- a/R/format.data/BCI.R +++ b/R/format.data/BCI.R @@ -9,40 +9,52 @@ library(reshape,quietly=TRUE) ############# READ DATA read individuals tree data Requires careful formatting of 7 census ############ datasets The raw data is such that, once a tree dies in census X, then it no ########## longer exists in census X+1, X+2 etc... -data.bci1 <- read.table("data/raw/BCI/census1/PlotsDataReport.txt", header = TRUE, - stringsAsFactors = FALSE, sep = "\t") +data.bci1 <- read.csv("data/raw/BCI/BCI.all.plots.stem.census.1.csv", header = TRUE, + stringsAsFactors = FALSE) data.bci1$Date1 <- data.bci1$Date; data.bci1$Date <- NULL data.bci1$DBH1 <- data.bci1$DBH; data.bci1$DBH <- NULL +data.bci1$HOM1 <- data.bci1$HOM; data.bci1$HOM <- NULL +data.bci1$Codes1 <- data.bci1$Codes; data.bci1$Codes <- NULL data.bci1$Status1 <- data.bci1$Status big.bci <- NULL max_census = 7 # only use first census for now for (k in 2:max_census) { - new.directory <- paste("data/raw/BCI/census", k, "/PlotsDataReport.txt", - sep = "") - data.bci2 <- read.table(new.directory, header = TRUE, stringsAsFactors = FALSE, - sep = "\t"); + new.directory <- paste("data/raw/BCI/BCI.all.plots.stem.census",k,"csv",sep = ".") + data.bci2 <- read.csv(new.directory, header = TRUE, stringsAsFactors = FALSE); + data.bci2$ID <- paste(data.bci2$Plot,data.bci2$TreeID,data.bci2$StemID,sep="") + if (!is.null(big.bci)) { - sub.bci <- merge(data.bci1[, c("Latin","Quadrat","Census","gx","gy","TreeID","Tag","Date1","DBH1","Status1")], data.frame(TreeID = data.bci2[["TreeID"]], - DBH2 = data.bci2[["DBH"]], Date2 = data.bci2[["Date"]], dead = as.numeric(data.bci2[["Status"]] == - "dead"),Status2=data.bci2[["Status"]],stringsAsFactors=F), sort = T, by = "TreeID") ## Uses the Date1 as the census number + data.bci1$ID <- paste(data.bci1$Plot,data.bci1$TreeID,data.bci1$StemID,sep="") + sub.bci <- merge(data.bci1[, c("Latin","Plot","Quadrat","Census","gx","gy","TreeID", + "StemID","ID","Stem","Date1","DBH1","HOM1","Codes1","Status1")], + data.frame(ID = data.bci2[["ID"]], + DBH2 = data.bci2[["DBH"]],HOM2 = data.bci2[["HOM"]], Codes2 = data.bci2[["Codes"]], + Date2 = data.bci2[["Date"]], dead = as.numeric(data.bci2[["Status"]] == + "dead"),Status2=data.bci2[["Status"]],stringsAsFactors=F), sort = T, by = "ID") ## Uses the Date1 as the census number big.bci <- rbind(big.bci, sub.bci) } if (is.null(big.bci)) { - big.bci <- merge(data.bci1[, c("Latin","Quadrat","Census","gx","gy","TreeID","Tag","Date1","DBH1","Status1")], data.frame(TreeID = data.bci2[["TreeID"]], - DBH2 = data.bci2[["DBH"]], Date2 = data.bci2[["Date"]], dead = as.numeric(data.bci2[["Status"]] == - "dead"),Status2=data.bci2[["Status"]],stringsAsFactors=F), sort = T, by = "TreeID") + data.bci1$ID <- paste(data.bci1$TreeID,data.bci1$StemID,sep="") + big.bci <- merge(data.bci1[, c("Latin","Plot","Quadrat","Census","gx","gy","TreeID", + "StemID","ID","Stem","Date1","DBH1","HOM1","Codes1","Status1")], + data.frame(ID = data.bci2[["ID"]], + DBH2 = data.bci2[["DBH"]],HOM2 = data.bci2[["HOM"]], Codes2 = data.bci2[["Codes"]], + Date2 = data.bci2[["Date"]], dead = as.numeric(data.bci2[["Status"]] == + "dead"),Status2=data.bci2[["Status"]],stringsAsFactors=F), sort = T, by = "ID") } data.bci1 <- data.bci2 data.bci1$Date1 <- data.bci1$Date; data.bci1$Date <- NULL data.bci1$DBH1 <- data.bci1$DBH; data.bci1$DBH <- NULL + data.bci1$HOM1 <- data.bci1$HOM; data.bci1$HOM <- NULL + data.bci1$Codes1 <- data.bci1$Codes; data.bci1$Codes <- NULL data.bci1$Status1 <- data.bci1$Status - cat("Census", k, "now included\n") + cat("Census", k, "now included. Dim big.bci",dim(big.bci),"\n") # print(summary(big.bci$DBH1)); print(summary(big.bci$DBH2)) } rm(data.bci1, data.bci2, sub.bci) -big.bci <- big.bci[order(big.bci$TreeID), ] +big.bci <- big.bci[order(big.bci$ID), ] data.bci <- big.bci rm(big.bci) ### read species names @@ -65,37 +77,9 @@ data.bci <- merge(data.bci, data.frame(Latin = species.clean$Latin_name, sp = sp ## format sp as sp. data.bci$sp <- paste("sp",data.bci$sp,sep=".") -############# READ ADDITIONAL PLOT -data.bci1 <- read.csv("data/raw/BCI/censusextra/data.small.plot1.csv", header = TRUE, stringsAsFactors = FALSE) -data.bci1$Date1 <- data.bci1$Date; data.bci1$Date <- NULL -data.bci1$DBH1 <- data.bci1$DBH; data.bci1$DBH <- NULL -data.bci1$Status1 <- data.bci1$Status -big.bci <- NULL - -data.bci2 <- read.csv("data/raw/BCI/censusextra/data.small.plot2.csv", header = TRUE, stringsAsFactors = FALSE) -big.bci <- merge(data.bci1[, c("Latin","Quadrat","Census","gx","gy","TreeID","Tag","Date1","DBH1","Status1")], data.frame(TreeID = data.bci2[["TreeID"]], - DBH2 = data.bci2[["DBH"]], Date2 = data.bci2[["Date"]], dead = as.numeric(data.bci2[["Status"]] == "dead"), - Status2=data.bci2[["Status"]],stringsAsFactors=F), sort = T, by = "TreeID") ## Uses the Date1 as the census number - -rm(data.bci1, data.bci2) -big.bci <- big.bci[order(big.bci$TreeID), ] -data.bci.extra <- big.bci -rm(big.bci) -### read species names -species.clean <- read.table("data/raw/BCI/TaxonomyDataReport.txt", stringsAsFactors = FALSE, - header = T, sep = "\t") -species.clean$Latin_name <- paste(species.clean[["Genus"]],species.clean[["species"]],sep=" ") -## ## Try to relate SpeciesID in species.clean species names in data.bci -#unique(data.bci$Latin) %in% species.clean$Latin_name -data.bci.extra <- merge(data.bci.extra, data.frame(Latin = species.clean$Latin_name, sp = species.clean$SpeciesID, genus = species.clean$Genus, stringsAsFactors = F),by = "Latin", sort = F) - -## format sp as sp. -data.bci.extra$sp <- paste("sp",data.bci.extra$sp,sep=".") -######## DEFINE THE ADDITIONAL DATASET AS AN DIFFERENT CLUSTER -data.bci$cluster <- rep("A",nrow(data.bci)) -data.bci.extra$cluster <- rep("B",nrow(data.bci.extra)) -data.bci <- rbind(data.bci, data.bci.extra) +######## DEFINE THE PLOTS AS DIFFERENT CLUSTER +data.bci$cluster <- data.bci$Plot ########################################## FORMAT INDIVIDUAL TREE DATA data.bci <- data.bci[order(data.bci[["TreeID"]]),] @@ -105,13 +89,17 @@ data.bci$Date2 <- as.Date(data.bci$Date2) # data.bci$yr2 <- format(strptime(data.bci$Date2, format = '%Y-%m-%d'),'%Y') data.bci$year <- as.numeric(difftime(data.bci$Date2, data.bci$Date1, units = "weeks")/52) ## Not rounded data.bci$obs.id <- 1:nrow(data.bci) #apply(data.bci[,c("TreeID","Census")],1,paste,collapse="_") -data.bci$tree.id <- data.bci$TreeID; data.bci$TreeID <- NULL +data.bci$tree.id <- data.bci$ID; data.bci$x <- data.bci$gx; data.bci$gx <- NULL data.bci$y <- data.bci$gy; data.bci$gy <- NULL data.bci$G <- (data.bci$DBH2 - data.bci$DBH1)/data.bci$year ## diameter growth in mm per year - BASED ON UNROUNDED YEARS data.bci$BA.G <- (pi*(data.bci$DBH2/20)^2 -pi*(data.bci$DBH1/20)^2)/data.bci$year ## BA growth in cm2/yr data.bci$G[data.bci$Status1=="missing" | data.bci$Status2=="missing"] <- NA +data.bci$G[abs(data.bci$HOM2-data.bci$HOM1)>0.5 & !is.na(data.bci$HOM2-data.bci$HOM1)] <- NA +data.bci$G[data.bci$Codes2 %in% c('RF','RP','RT','R')] <- NA data.bci$BA.G[data.bci$Status1=="missing" | data.bci$Status2=="missing"] <- NA +data.bci$BA.G[abs(data.bci$HOM2-data.bci$HOM1)>0.5 & !is.na(data.bci$HOM2-data.bci$HOM1)] <- NA +data.bci$BA.G[data.bci$Codes2 %in% c('RF','RP','RT','R')] <- NA data.bci$G[data.bci$Latin %in% sp.palm.fern] <- NA ## remove fern and palm data.bci$BA.G[data.bci$Latin %in% sp.palm.fern] <- NA ## remove fern and palm data.bci$dead[data.bci$Status1=="missing" | data.bci$Status2=="missing"] <- NA diff --git a/R/process.data/test.tree.CWM.R b/R/process.data/test.tree.CWM.R index 6e23fdd03d484fe23bf986104a716eb94d159df7..6453726ea22b6468b4b4310f995e87f14b787161 100644 --- a/R/process.data/test.tree.CWM.R +++ b/R/process.data/test.tree.CWM.R @@ -85,4 +85,3 @@ dir.create("figs/test.processed", recursive=TRUE,showWarnings=FALSE) lapply(sets,fun.test.set.tree.CWM,filedir=filedir) -#### TODO NEED TO CHECK WARNINGS ROWNAMES DUPLICATED WHY ?