Skip to content

GitLab

Commit 5d4ce8a5 authored by Georges Kunstler's avatar Georges Kunstler
Browse files

start rstan model for global data

parent 8ec9cedd
Hide whitespace changes
Inline Side-by-side
R/analysis/lmer.run.nolog.all.R
View file @ 5d4ce8a5
......@@ -49,8 +49,8 @@ fun.call.lmer <- function(formula,df.lmer){
fun.call.lmer.and.save <- function(formula,df.lmer,path.out){
browser()
lmer.output <- lmer(formula=formula,data=df.lmer,REML = FALSE)
lmer.output <- lmer(formula=formula,data=df.lmer,REML = FALSE,
control=lmerControl(optCtrl=list(maxfun=20000) ) )
print(summary(lmer.output))
saveRDS(lmer.output,file=file.path(path.out, "results.nolog.all.rds"))
}
......@@ -106,17 +106,20 @@ fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,perc.neigh,BATOT
data.tree <- subset(data.tree,subset=(!is.na(data.tree[["BA.G"]])) &
(!is.na(data.tree[["D"]])) )
## remove tree with zero
data.tree <- subset(data.tree,subset=data.tree[["BA.G"]]>min.BA.G & data.tree[["BA.G"]]<max.BA.G &
data.tree[["D"]]>9.9)
data.tree <- subset(data.tree,subset=data.tree[["BA.G"]]>min.BA.G &
data.tree[["BA.G"]]<max.BA.G &
data.tree[["D"]]>9.9)
## select species with no missing traits
select.temp <- (!is.na(data.tree[[abs.CWM.tntf]]) &
!is.na(data.tree[[BATOT]]))
data.tree <- data.tree[select.temp,]
# select obs abov min perc.neigh
data.tree <- subset(data.tree,subset=data.tree[[perc.neigh]] > min.perc.neigh &
!is.na(data.tree[[perc.neigh]]))
# 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.neigh
data.tree <- subset(data.tree,subset=data.tree[[perc.neigh]] > min.perc.neigh & !is.na(data.tree[[perc.neigh]]))
names(table(factor(data.tree[["sp"]])))[
table(factor(data.tree[["sp"]]))>min.obs])
return(data.tree)
}
......@@ -124,37 +127,33 @@ fun.center.and.standardized.var <- function(x){
return((x-mean(x))/sd(x))
}
### TODO NEED TO CONTIMNUE TO CHANGE FROM HERE
### get variables for lmer
fun.get.the.variables.for.lmer.tree.id <- function(data.tree,BATOT, CWM.tn,
abs.CWM.tntf, tf,
min.BA.G=50){
logG <- fun.center.and.standardized.var(log(data.tree[["BA.G"]]+min.BA.G))
logD <- log(data.tree[["D"]])
species.id <- unclass(factor(paste(data.tree[['set']],
data.tree[["sp"]],sep='.')))
tree.id <- unclass(factor(paste(data.tree[['set']],data.tree[["tree.id"]])))
plot.species.id <- unclass(factor(paste(data.tree[["set"]],
data.tree[["plot"]],
data.tree[["sp"]],sep="")))
set.id <- unclass(factor(data.tree[["set"]]))
ecocode.id <- unclass(factor(paste(data.tree[["set"]],data.tree[['ecocode']])))
#= multiply CWMs by BATOT
sumTnTfBn.abs <- data.tree[[abs.CWM.tntf]]*data.tree[[BATOT]]
sumTnBn <- data.tree[[CWM.tn]]*data.tree[[BATOT]]
sumTfBn <- data.tree[[tf]]*data.tree[[BATOT]]
logD <- fun.center.and.standardized.var(log(data.tree[["D"]]))
species.id <- factor(data.tree[["sp"]])
tree.id <- factor(data.tree[["tree.id"]])
plot.id <- factor(data.tree[["plot"]])
set.id <- factor(data.tree[["set"]])
ecocode.id <- factor(data.tree[['ecocode']])
#= multiply CWMs by BATOT
sumTnTfBn.abs <- data.tree[[abs.CWM.tntf]]*
data.tree[[BATOT]]/max(data.tree[[BATOT]])
sumTnBn <- data.tree[[CWM.tn]]*data.tree[[BATOT]]/max(data.tree[[BATOT]])
sumTfBn <- data.tree[[tf]]*data.tree[[BATOT]]/max(data.tree[[BATOT]])
sumTnTfBn.diff <- sumTnBn-sumTfBn
sumBn <- data.tree[[BATOT]]
sumBn <- data.tree[[BATOT]]/max(data.tree[[BATOT]])
return(data.frame(logG= logG,
logD= logD,
species.id= species.id,
tree.id= tree.id,
set.id= set.id,
ecocode.id= ecocode.id,
plot.species.id= plot.species.id,
plot.id= plot.id,
sumTnTfBn.diff= sumTnTfBn.diff,
sumTnTfBn.abs= sumTnTfBn.abs,
Tf= data.tree[[tf]],
......
R/analysis/run.local.R
View file @ 5d4ce8a5
##### SCRIPT TO TEST BEFORE TO SEND ON CLUSTER
source("R/analysis/lmer.run.nolog.all.R")
run.lmer(model.files.lmer.Tf.1[3],'SLA',type.filling='species')
run.models.for.set.all.traits(model.files.lmer.Tf.1[3], run.lmer,type.filling='species')
trait <- "SLA"
type.filling <- 'species'
source(model.files.lmer.Tf.1[3], local = TRUE)
model <- load.model()
#= Path for output
path.out <- output.dir.lmer(model$name, trait, 'all',
type.filling=type.filling)
df.lmer <- load.and.prepare.data.for.lmer(trait,
min.obs, sample.size,
type.filling=type.filling) # return a DF
lmer.output.s <- lmer(formula=model$lmer.formula.tree.id,data=df.lmer[sample(1:nrow(df.lmer),200000), ],
REML = FALSE,
control=lmerControl(optCtrl=list(maxfun=10000) ) )
lmer.output <- lmer(formula=model$lmer.formula.tree.id,data=df.lmer,
REML = FALSE,
control=lmerControl(optCtrl=list(maxfun=10000) ) )
lmer.output.s@optinfo$conv
lmer.output@optinfo$conv
##### SCRIPT TO TEST BEFORE TO SEND ON CLUSTER
source("R/analysis/lmer.run.nolog.R")
......
R/analysis/stan-fun.R
View file @ 5d4ce8a5
......@@ -3,32 +3,41 @@ library(rstan) #http://www.mc-stan.org/
stan.model <- "
data {
int<lower=0> N_indiv;
int<lower=0> N_species;
int<lower=0> N_plot_species;
int<lower=0> Ncompet_sp;
int<lower=0> N.indiv;
int<lower=0> N.species;
int<lower=0> N.plot;
int<lower=0> N.ecocode;
int[N_species] species_id;
int[N_plot_species] plot_species_id;
int[N.indiv] species.id;
int[N.indiv] plot.id;
int[N.indiv] ecocode.id;
real[N_indiv] logG;
real[N_indiv] logD;
real[N_indiv] logA1;
real[N_indiv] logA2;
matrix[N_indiv,Ncompet_sp] data_compet;
#vector[] trait_compet;
#vector[] trait_target;
real[N.indiv] logG;
real[N.indiv] logD;
real[N.indiv] sumBn;
real[N.indiv] Tf;
real[N.indiv] sumTfBn;
real[N.indiv] sumTnBn;
}
parameters{
real<lower=-1,upper=1> intercept_trait;
real<lower=-10,upper=10>[3] mean_B;
real<lower=0.0001,upper=5>[3] sigma_B;
real<lower=-1,upper=1> intercept;
real<lower=-10,upper=10> mean.logD;
real<lower=-10,upper=10> mean.Tf;
real<lower=-10,upper=10> mean.sumBn;
real<lower=-10,upper=10> mean.sumTfBn;
real<lower=-10,upper=10> mean.sumTnBn;
real<lower=0.0001,upper=5> sigma.inter.sp;
real<lower=0.0001,upper=5> sigma.inter.ecocode;
real<lower=0.0001,upper=5> sigma.inter.plot;
real<lower=0.0001,upper=5> sigma.logD.sp;
real<lower=0.0001,upper=5> sigma.Tf.ecocode;
real<lower=0.0001,upper=5> sigma.sumBn.ecocode;
real<lower=0.0001,upper=5> sigma.sumTfBn.ecocode;
real<lower=0.0001,upper=5> sigma.sumTnBn.ecocode;
real<lower=0.0001,upper=5> sigma;
real<lower=0.0001,upper=5> sigma_plot;
real<lower=0.0001,upper=5> sigma_alpha;
}
......@@ -39,88 +48,123 @@ model {
### species random param
for (n in 1:N_species)
for (n in 1:N.species)
{
B1[n] ~ dnorm(mean_B[1],sigma_B[1])
B2[n] ~ dnorm(mean_B[2],sigma_B[2])
B3[n] ~ dnorm(mean_B[3],sigma_B[3])
alpha[n] ~ dnorm(0,sigma_alpha)
for (j in 1:Ncompet_sp)
{
param_compet[n,j] <- intercept_trait
}
param.logD.sp[n] ~ dnorm(mean.logD,sigma.logD.sp);
intercept.sp[n] ~ dnorm(intercept,sigma.inter.sp);
}
############ Likelihood ###################
for (i in 1:N_indiv) {
### plot random param
for (n in 1:N.plot)
{
intercept.plot[n] ~ dnorm(0,sigma.inter.plot)
}
IC[i] <- dot.product(
log(1+data_compet[i,1:Ncompet_sp]),
param_compet[species.id[i],1:Ncompet_sp]
);
### ecocode random param
for (n in 1:N.ecocode)
{
param.Tf.ecocode[n] ~ dnorm(mean.Tf,sigma.Tf.ecocode);
param.sumBn.ecocode[n] ~ dnorm(mean.sumBn,sigma.sumBn.ecocode);
param.sumTfBn.ecocode[n] ~ dnorm(mean.sumTfBn,sigma.sumTfBn.ecocode);
param.sunTnBn.ecocode[n] ~ dnorm(meanTnBn,sigma.sumTnBn.ecocode);
intercept.ecocode[n] ~ dnorm(0,sigma.inter.ecocode);
}
real theo_g;
theo_g <- alpha[species_id[i]]
+ B1[species_id[i]]*logD[i]
+ B2[species_id[i]]*logA1[i]
+ B3[species_id[i]]*logA2[i]
+IC[i]
+random_plot[plot_species_id[i]]
############ Likelihood ###################
for (i in 1:N.indiv) {
real theo.g;
theo.g <- intercept.sp[species.id[i]]
+ intercept.plot[plot.id[i]]
+ intercept.ecocode[ecocode.id[i]]
+ param.logD.sp[species.id[i]]*logD[i]
+ param.Tf.ecocode[ecocode.id[i]]*Tf[i]
+ param.sumBn.ecocode[ecocode.id[i]]*sumBn[i]
+ param.sumTfBn.ecocode[ecocode.id[i]]*sumTfBn[i]
+ param.sumTnBn.ecocode[ecocode.id[i]]*sumTnBn[i]
;
logG ~ dnorm(theo_g,sigma);
logG ~ dnorm(theo.g,sigma);
}
################################################
########### Hierarchical parameters ########
## plot effect
for (j in 1:N_plot_species)
{
random_plot[j] ~ dnorm(0,sigma_plot)
}
###############################################
########### Non-hierarchical parameters ########
# constants for prior
real sigma0;
sigma0 <- pow(1.0E-4,-2)
sigma0 <- pow(1.0E-2,-2);
# slope and intercept
intercept ~ dnorm(0,sigma0);
mean.logD ~ dnorm(0,sigma0);
mean.Tf ~ dnorm(0,sigma0);
mean.sumBn ~ dnorm(0,sigma0);
mean.sumTfBn ~ dnorm(0,sigma0);
mean.sumTnBn ~ dnorm(0,sigma0);
# sigma
sigma.inter.sp~ uniform(0,6);
sigma.inter.ecocode~ uniform(0,6);
sigma.inter.plot~ uniform(0,6);
sigma.logD.sp~ uniform(0,6);
sigma.Tf.ecocode~ uniform(0,6);
sigma.sumBn.ecocode~ uniform(0,6);
sigma.sumTfBn.ecocode~ uniform(0,6);
sigma.sumTnBn.ecocode~ uniform(0,6);
sigma~ uniform(0,6);
## traits parameters
intercept_trait ~ dnorm(0,sigma0)
## prior for hyper parameter mean and sd
for (l in 1:3)
{
mean_B[l]~ dnorm(0,sigma0)
}
} # End of the jags model
} # End of the stan model
}
}
"
stan.dat <- jags.data
stan.inits.LOGLIN.NULL <- function(){
list("intercept_trait"=runif(1,min=-1,max=1),
"mean_B"=runif(3,min=-0.5,max=1),
"sigma_B"=runif(3,min=0.0001,max=1.5),
list("intercept"=runif(1,min=-1,max=1),
"mean.logD"=runif(1,min=-0.5,max=1),
"mean.Tf"=runif(1,min=-0.5,max=1),
"mean.sumBn"=runif(1,min=-0.5,max=1),
"mean.sumTfBn"=runif(1,min=-0.5,max=1),
"mean.sumTnBn"=runif(1,min=-0.5,max=1),
"sigma.inter.sp"=runif(1,min=0.0001,max=1.5),
"sigma.logD.sp"=runif(1,min=0.0001,max=1.5),
"sigma"=runif(1,min=0.0001,max=4),
"sigma_alpha"=runif(1,min=0.0001,max=4),
"sigma_plot"=runif(1,min=0.0001,max=4))
"sigma.inter.ecocode"=runif(1,min=0.0001,max=4),
"sigma.Tf.ecocode"=runif(1,min=0.0001,max=4),
"sigma.sumBn.ecocode"=runif(1,min=0.0001,max=4),
"sigma.sumTfBn.ecocode"=runif(1,min=0.0001,max=4),
"sigma.sumTnBn.ecocode"=runif(1,min=0.0001,max=4),
"sigma.inter.plot"=runif(1,min=0.0001,max=4))
}
# load data
source("R/analysis/lmer.run.nolog.all.R")
trait <- "SLA"
type.filling <- 'species'
df.lmer <- load.and.prepare.data.for.lmer(trait,
min.obs, sample.size,
type.filling=type.filling) # return a DF
stan.dat <- list(N.indiv=nrow(df.lmer),
N.species=nlevels(df.lmer$species.id),
N.plot=nlevels(df.lmer$plot.id),
N.ecocode=nlevels(df.lmer$ecocode.id),
species.id=unclass(factor(df.lmer$species.id)),
plot.id=unclass(factor(df.lmer$plot.id)),
species.id=unclass(factor(df.lmer$ecocode.id)),
logG=df.lmer$logG,
logD=df.lmer$logD,
Tf=df.lmer$Tf,
sumBn=df.lmer$sumBn,
sumTfBn=df.lmer$sumTfBn,
sumTnBn=df.lmer$sumTnBn)
#warmup = burn-in
stan.out <- stan(
model_code=stan.code,
......
R/format.data/Japan.R
View file @ 5d4ce8a5
......@@ -105,7 +105,7 @@ data.japan$cluster <- data.japan$plot ## cluster number
data.japan$plot <- NA ## no plot
data.japan$htot <- rep(NA,nrow(data.japan)) ## height of tree in m - missing
data.japan$tree.id <- data.japan$id; ## tree unique id
data.japan$tree.id <- paste(data.japan$cluster,data.japan$id) ## tree unique id
data.japan$obs.id <- 1:nrow(data.japan) ## obs uniquue id
data.japan <- subset(data.japan,subset= (data.japan$dead1 == 0))
......@@ -113,8 +113,33 @@ data.japan$dead1 <- NULL
data.japan$dead <- data.japan$dead2
data.japan$dead2 <- NULL
### TODO CHANGE PLOT SIZE AS other
#remove tree with no x y
data.japan <- data.japan[!(is.na(data.japan$x) | is.na(data.japan$y)),]
## create small plot
square.s.t <- 20
make.quad <- do.call("rbind",lapply(unique(data.japan$cluster),
FUN=fun.quadrat.cluster,
cluster=data.japan$cluster,
tree.id=data.japan$obs.id,
x=data.japan$x,
y=data.japan$y,
square.s=square.s.t))
data.japan <- merge(data.japan, make.quad,by.x='obs.id',by.y="tree.id")
## #plot each plot
## pdf("figs/plots.japan.pdf")
## lapply(unique(data.japan[["cluster"]]),FUN=fun.circles.plot,
## data.japan[['x']],data.japan[['y']],
## data.japan[["cluster"]],data.japan[["D"]],
## inches=0.2,fg.l=data.japan$make.quad)
## dev.off()
#
data.japan$plot <- data.japan$make.quad
## LOAD WORLDCLIM
#### get wc climate
......
R/format.data/Paracou.R
View file @ 5d4ce8a5
......@@ -9,23 +9,30 @@ library(reshape,quietly=TRUE)
############################ read individuals tree data
data.paracou <- read.table("data/raw/Paracou/20130717_paracou_1984_2012.csv",
header=TRUE,stringsAsFactors=FALSE,sep = ";", na.strings = "NULL")
#barplot(apply(!is.na(data.paracou[,paste("circ_",1984:2012,sep="")]),MARGIN=2,FUN=sum),las=3)
header=TRUE,stringsAsFactors=FALSE,sep = ";",
na.strings = "NULL")
# select good columns
data.paracou <- data.paracou[,c("foret","parcelle","carre","arbre","vernaculaire","idtaxon",
"x","y","circ_2001","code_2001","circ_2005","code_2005",
"circ_2009","code_2009","campagne_mort","type_mort")]
colnames(data.paracou) <- c("forest","cluster","plot","tree","vernacular","taxonid","x","y",
"circum2001","code2001","circum2005","code2005","circum2009",
data.paracou <- data.paracou[,c("foret","parcelle","carre","arbre",
"vernaculaire","idtaxon",
"x","y","circ_2001","code_2001",
"circ_2005","code_2005",
"circ_2009","code_2009","campagne_mort",
"type_mort")]
colnames(data.paracou) <- c("forest","cluster","plot","tree","vernacular",
"taxonid","x","y",
"circum2001","code2001","circum2005","code2005",
"circum2009",
"code2009","yeardied","typedeath")
### change numeric separator
numeric.col.name <- c("x","y","circum2001","code2001","circum2005",
"code2005","circum2009","code2009")
for(k in numeric.col.name){
data.paracou[,k] <- gsub(",",".",data.paracou[,k]); data.paracou[,k] <- as.numeric(data.paracou[,k])
data.paracou[,k] <- gsub(",",".",data.paracou[,k])
data.paracou[,k] <- as.numeric(data.paracou[,k])
} ## Replace all , in decimals with .
data.paracou$tree.id <- apply(data.paracou[,c("cluster","plot","tree")],1,paste,collapse="_");
data.paracou$tree.id <- apply(data.paracou[,c("cluster","plot","tree")],
1,paste,collapse="_");
data.paracou$sp <- paste("sp",data.paracou[["taxonid"]],sep=".")
data.paracou <- data.paracou[,c(ncol(data.paracou),1:(ncol(data.paracou)-1))]
......@@ -36,14 +43,19 @@ data.paracou <- subset(data.paracou,subset=(!is.na(data.paracou[["y"]])) &
data.paracou[["x"]]<251 &
data.paracou[["y"]]<251)
#### REMOVE PLOTs 16 17 18 ACCORDING TO GHSILAIN
data.paracou <- subset(data.paracou,subset=! data.paracou[["cluster"]] %in% 16:18)
data.paracou <- subset(data.paracou,
subset=! data.paracou[["cluster"]] %in% 16:18)
## keep only tree alive in 2001
data.paracou <- subset(data.paracou,subset=!(as.numeric(data.paracou[["yeardied"]])<=2001 & !is.na(data.paracou[["yeardied"]])))
data.paracou <- subset(data.paracou,
subset=!(as.numeric(data.paracou[["yeardied"]])<=2001 &
!is.na(data.paracou[["yeardied"]])))
## plot each plot
## pdf("figs/plots.paracou.pdf")
## lapply(unique(data.paracou[["cluster"]]),FUN=fun.circles.plot,data.paracou[['x']],data.paracou[['y']],
## data.paracou[["cluster"]],data.paracou[["circum2009"]],inches=0.2,fg.l=data.paracou$plot)
## lapply(unique(data.paracou[["cluster"]]),FUN=fun.circles.plot,
## data.paracou[['x']],data.paracou[['y']],
## data.paracou[["cluster"]],data.paracou[["circum2009"]],
## inches=0.2,fg.l=data.paracou$plot)
## dev.off()
## #
......@@ -70,54 +82,82 @@ data.paracou <- merge(data.paracou, make.quad,by='tree.id')
######################################## MASSAGE TRAIT DATA
######################################## MASSAGE TREE DATA
### read species names
species.clean <- read.csv("data/raw/Paracou/20130717_paracou_taxonomie.csv",stringsAsFactors=FALSE, header = T, sep = ";")
species.clean <- read.csv("data/raw/Paracou/20130717_paracou_taxonomie.csv",
stringsAsFactors=FALSE, header = T, sep = ";")
species.clean$sp <- paste("sp.",species.clean[["idTaxon"]],sep="")
species.clean$Latin_name <- paste(species.clean[["Genre"]],species.clean[["Espece"]],sep=" ")
species.clean$Latin_name <- paste(species.clean[["Genre"]],
species.clean[["Espece"]],sep=" ")
## keep only one row pers idTaxon
species.clean <- subset(species.clean,subset=!duplicated(species.clean[["sp"]]),select=c("sp","Latin_name","Genre","Espece","Famille","idCIRAD"))
species.clean <- subset(species.clean,subset=!duplicated(species.clean[["sp"]]),
select=c("sp","Latin_name","Genre","Espece",
"Famille","idCIRAD"))
## select only species present in data base
species.clean <- subset(species.clean,subset=species.clean[["sp"]] %in% data.paracou[["sp"]])
species.clean <- subset(species.clean,subset=species.clean[["sp"]] %in%
data.paracou[["sp"]])
## percentage of species with no taxonomic identification
length(grep("Indet",species.clean[["Latin_name"]]))/nrow(species.clean) ## 25% agree with Bruno ?
length(grep("Indet",species.clean[["Latin_name"]]))/nrow(species.clean)
## 25% agree with Bruno ?
## GET PALM AND TREE FERN TO REMOVE THEIR GROWTH
sp.palm.fern <- species.clean$sp[species.clean$Famille %in% c('Arecaceae','Cyatheaceae','Dicksoniaceae',
'Metaxyaceae', 'Cibotiaceae',
'Loxomataceae', 'Culcitaceae',
'Plagiogyriaceae', 'Thyrsopteridaceae')]
sp.palm.fern <- species.clean$sp[species.clean$Famille %in%
c('Arecaceae','Cyatheaceae','Dicksoniaceae',
'Metaxyaceae', 'Cibotiaceae',
'Loxomataceae', 'Culcitaceae',
'Plagiogyriaceae', 'Thyrsopteridaceae')]
############################################ FORMAT INDIVIDUAL TREE DATA
data.paracou2 <- data.paracou[rep(1:nrow(data.paracou),each=2),c(1:10,(ncol(data.paracou)-2):ncol(data.paracou))]
rownames(data.paracou2) <- 1:nrow(data.paracou2); data.paracou2 <- as.data.frame(data.paracou2)
data.paracou2$yr1 <- rep(c(2001,2001+4),nrow(data.paracou)); data.paracou2$yr2 <- rep(c(2005,2005+4),nrow(data.paracou))
data.paracou2 <- data.paracou[rep(1:nrow(data.paracou),each=2),
c(1:10,(ncol(data.paracou)-2):ncol(data.paracou))]
rownames(data.paracou2) <- 1:nrow(data.paracou2)
data.paracou2 <- as.data.frame(data.paracou2)
data.paracou2$yr1 <- rep(c(2001,2001+4),nrow(data.paracou))
data.paracou2$yr2 <- rep(c(2005,2005+4),nrow(data.paracou))
data.paracou2$year <- rep(c(4,4),nrow(data.paracou))
data.paracou2$dbh1 <- c(rbind(data.paracou$circum2001/pi,data.paracou$circum2005/pi))
data.paracou2$dbh2 <- c(rbind(data.paracou$circum2005/pi,data.paracou$circum2009/pi))
data.paracou2$code1 <- c(as.numeric(rbind(data.paracou$code2001,data.paracou$code2005)))
data.paracou2$code2 <- c(as.numeric(rbind(data.paracou$code2005,data.paracou$code2009)))
data.paracou2$dbh1 <- c(rbind(data.paracou$circum2001/pi,
data.paracou$circum2005/pi))
data.paracou2$dbh2 <- c(rbind(data.paracou$circum2005/pi,
data.paracou$circum2009/pi))
data.paracou2$code1 <- c(as.numeric(rbind(data.paracou$code2001,
data.paracou$code2005)))
data.paracou2$code2 <- c(as.numeric(rbind(data.paracou$code2005,
data.paracou$code2009)))
data.paracou2$dead <- rep(0,nrow(data.paracou)*2)
data.paracou2$dead[c(as.numeric(data.paracou[["yeardied"]]) %in% 2002:2005 &
(!is.na(data.paracou[["yeardied"]])),
as.numeric(data.paracou[["yeardied"]]) %in% 2006:2009 &
(!is.na(data.paracou[["yeardied"]])))] <- 1
## remove tree dead at first census for both date (census 2001-2005 2005-2009)
data.paracou <- subset(data.paracou2,subset=!(data.paracou2[['yr1']] ==2005 & (as.numeric(data.paracou[["yeardied"]]) %in% 2002:2005 & (!is.na(data.paracou[["yeardied"]])))))
data.paracou <- subset(data.paracou2,
subset=!(data.paracou2[['yr1']] ==2005 &
(as.numeric(data.paracou[["yeardied"]]) %in%
2002:2005 & (!is.na(data.paracou[["yeardied"]])))))
## change unit and names of variables to be the same in all data for the tree
data.paracou$G <- 10*(data.paracou$dbh2-data.paracou$dbh1)/data.paracou$year ## diameter growth in mm per year
data.paracou$G[data.paracou$code1>0] <- NA ## indivs with code indicating problem in dbh measurment at dbh1
data.paracou$G[data.paracou$code2>0] <- NA ## indivs with code indicating problem in dbh measurment at dbh2
data.paracou$BA.G <- (pi*(data.paracou$dbh2/2)^2-pi*(data.paracou$dbh1/2)^2)/data.paracou$year ## BA growth in cm2 per year
data.paracou$BA.G[data.paracou$code1>0] <- NA ## indivs with code indicating problem in dbh measurment at dbh1
data.paracou$BA.G[data.paracou$code2>0] <- NA ## indivs with code indicating problem in dbh measurment at dbh2
data.paracou$G[data.paracou$sp %in% sp.palm.fern] <- NA ## remove palm and fern growth
data.paracou$BA.G[data.paracou$sp %in% sp.palm.fern] <- NA ## remove palm and fern growth
data.paracou$G <- 10*(data.paracou$dbh2-data.paracou$dbh1)/data.paracou$year
## diameter growth in mm per year
data.paracou$G[data.paracou$code1>0] <- NA
## indivs with code indicating problem in dbh measurment at dbh1
data.paracou$G[data.paracou$code2>0] <- NA
## indivs with code indicating problem in dbh measurment at dbh2
data.paracou$BA.G <- (pi*(data.paracou$dbh2/2)^2-pi*(data.paracou$dbh1/2)^2)/
data.paracou$year ## BA growth in cm2 per year
data.paracou$BA.G[data.paracou$code1>0] <- NA
## indivs with code indicating problem in dbh measurment at dbh1
data.paracou$BA.G[data.paracou$code2>0] <- NA
## indivs with code indicating problem in dbh measurment at dbh2
data.paracou$G[data.paracou$sp %in% sp.palm.fern] <- NA
## remove palm and fern growth
data.paracou$BA.G[data.paracou$sp %in% sp.palm.fern] <- NA
## remove palm and fern growth
data.paracou$census <- data.paracou$yr1
data.paracou$D <- data.paracou[["dbh1"]]; data.paracou$D[data.paracou$D == 0] <- NA ;## diameter in cm
data.paracou$cluster <- paste("p",data.paracou$cluster,sep=".")#apply(data.paracou[,c("forest","cluster","plot")],1,paste,collapse=".") ## plot code
data.paracou$D <- data.paracou[["dbh1"]]
data.paracou$D[data.paracou$D == 0] <- NA ;## diameter in cm
data.paracou$cluster <- paste("p",data.paracou$cluster,sep=".")
data.paracou$htot <- rep(NA,length(data.paracou[["G"]])) ## height of tree in m
data.paracou$obs.id <- 1:nrow(data.paracou)
data.paracou$plot <- data.paracou$make.quad
......@@ -126,7 +166,8 @@ data.paracou <- subset(data.paracou,subset=!is.na(data.paracou$D))
## minimum circumfer 30 delete all tree with a dbh <30/pi,
data.paracou <- subset(data.paracou,subset= data.paracou[["D"]]>(30/pi))
## add latin name
data.paracou <- merge(data.paracou,subset(species.clean,select=c("sp","Latin_name")),by="sp")
data.paracou <- merge(data.paracou,subset(species.clean,
select=c("sp","Latin_name")),by="sp")
data.paracou$sp.name <- data.paracou$Latin_name
####################
......@@ -142,7 +183,8 @@ data.paracou$MAT <- clim$MAT
data.paracou$MAP <- clim$MAP
###################### PLOT SELECTION FOR THE ANALYSIS - NEEDS REDOING
vec.basic.var <- c("obs.i