Commit 3b4e3196 authored by Georges Kunstler's avatar Georges Kunstler
Browse files

start write lmer all data model

parent 87744325
###########################
###########################
### FUNCTION TO RUN LMER ESTIMATION
### FUNCTION TO RUN LMER ESTIMATION WITH NO logBA
library(lme4)
......
......@@ -80,10 +80,28 @@ data.cat.extract <- fun.change.factor.pheno.try(data.cat.extract)
data.cat.extract <- fun.change.factor.angio.try(data.cat.extract)
data.cat.extract <- fun.fill.pheno.try.with.zanne(data.cat.extract)
## READ pheno data from BCI
pheno.bci <- read.csv("data/raw/BCI/BCI.pheno.csv", stringsAsFactors = FALSE)
pheno.bci$Latin_name <- paste(pheno.bci$genus,pheno.bci$species,sep=" ")
## merge bci data
data.cat.extract <- merge(data.cat.extract,
subset(pheno.bci,select=c("Latin_name","deciduous")),
by="Latin_name",all.x=TRUE,all.y=FALSE)
write.csv(data.cat.extract,file="output/pheno.comparison.bci.joe.csv",row.names = FALSE)
data.cat.extract$deciduous[data.cat.extract$deciduous %in% c('DB','DF','D-NOVEMBER') &
!is.na(data.cat.extract$deciduous)] <- 'D_EV'
data.cat.extract$deciduous[data.cat.extract$deciduous %in% c('DO') &
!is.na(data.cat.extract$deciduous)] <- 'D'
data.cat.extract$deciduous[data.cat.extract$deciduous %in% c('E') &
!is.na(data.cat.extract$deciduous)] <- 'EV'
data.cat.extract$Pheno.T[!is.na(data.cat.extract$deciduous)] <-
data.cat.extract$deciduous[!is.na(data.cat.extract$deciduous)]
data.traits <- merge(data.traits,data.cat.extract[,c("sp","Phylo.group","Pheno.T")],by="sp")
###
write.csv(data.traits,file="output/formatted/BCI/traits.csv",row.names = FALSE)
......
......@@ -60,7 +60,6 @@ data.cat.extract[data.cat.extract$Latin_name %in% c('Machilus zuihoensis',
'Osmanthus matsumuranus','Symplocos wikstroemiifolia','Itea parviflora',
'Symplocos sonoharae','Cyclobalanopsis sessilifolia',
'Eriobotrya deflexa','Neolitsea konishii'),"Pheno.T"] <- 'EV'
data.cat.extract[is.na(data.cat.extract$Pheno.T),]
data.traits <- merge(data.traits,data.cat.extract[,c("sp","Phylo.group","Pheno.T")],by="sp")
......
......@@ -166,7 +166,7 @@ data.bci[["ecocode"]] <- "tropical"
###################### PLOT SELECTION FOR THE ANALYSIS
vec.basic.var <- c("obs.id","tree.id", "sp", "sp.name",
"cluster","plot", "ecocode", "D", "G","BA.G", "dead",
"cluster","plot", "ecocode", "D", "G","BA.G","year", "dead",
'Lon','Lat',"x", "y", "census",'MAT','MAP')
data.tree <- subset(data.bci, select = c(vec.basic.var))
......
......@@ -91,7 +91,7 @@ data.canada <- merge(data.canada, data.frame(plot = names(perc.dead), perc.dead
table(data.canada$dead)
## data.canada <- data.canada[data.canada$dead == 0,]
vec.abio.var.names <- c("MAT", "MAP")
vec.basic.var <- c("obs.id","tree.id", "sp","sp.name", "cluster", "plot", "ecocode", "D", "G","BA.G", "dead",
vec.basic.var <- c("obs.id","tree.id", "sp","sp.name", "cluster", "plot", "ecocode", "D", "G","BA.G","year", "dead",
"Lon", "Lat","weights","census")
data.tree <- subset(data.canada, select = c(vec.basic.var, vec.abio.var.names))
......
......@@ -82,7 +82,7 @@ data.france$MAP <- data.france$SAP
## names of variables for abiotic conditions
vec.abio.var.names <- c("MAT", "MAP")
## other var
vec.basic.var <- c("obs.id","tree.id", "sp","sp.name", "cluster", "plot", "ecocode", "D", "G","BA.G", "dead",
vec.basic.var <- c("obs.id","tree.id", "sp","sp.name", "cluster", "plot", "ecocode", "D", "G","BA.G","year", "dead",
"Lon", "Lat","weights","census")
data.tree <- subset(data.france, select = c(vec.basic.var, vec.abio.var.names))
## select only tree above 10cm DBH
......
......@@ -64,7 +64,7 @@ data.fushan$MAP <- clim$MAP
###################### PLOT SELECTION FOR THE ANALYSIS - NEEDS REDOING
vec.basic.var <- c("obs.id","tree.id", "sp", "sp.name", "cluster","plot", "ecocode", "D", "G","BA.G", "dead",
vec.basic.var <- c("obs.id","tree.id", "sp", "sp.name", "cluster","plot", "ecocode", "D", "G","BA.G","year", "dead",
'Lon','Lat',"x", "y", "census",'MAT','MAP')
data.tree <- subset(data.fushan, select = c(vec.basic.var))
......
......@@ -169,7 +169,7 @@ data.japan$MAP <- clim$MAP
###################### PLOT SELECTION FOR THE ANALYSIS
vec.basic.var <- c("obs.id","tree.id", "sp", "sp.name", "cluster","plot", "ecocode", "D", "G","BA.G", "dead",
vec.basic.var <- c("obs.id","tree.id", "sp", "sp.name", "cluster","plot", "ecocode", "D", "G","BA.G","year", "dead",
'Lon','Lat',"x", "y", "census",'MAT','MAP')
data.tree <- subset(data.japan, select = c(vec.basic.var))
......
......@@ -116,10 +116,10 @@ data.luq[["ecocode"]] <- "tropical"
###################### PLOT SELECTION FOR THE ANALYSIS - NEEDS REDOING
vec.basic.var <- c("obs.id","tree.id", "sp", "sp.name", "cluster","plot", "ecocode", "D", "G","BA.G", "dead",
vec.basic.var <- c("obs.id","tree.id", "sp", "sp.name", "cluster","plot", "ecocode", "D", "G","BA.G","year", "dead",
'Lon','Lat',"x", "y", "census",'MAT','MAP')
data.tree <- subset(data.luq, select = c(vec.basic.var))
y
data.tree <- subset(data.tree,subset=!is.na(data.tree$x) & !is.na(data.tree$y))
## convert var factor in character or numeric
data.tree <- fun.convert.type.B(data.tree)
......
......@@ -135,7 +135,7 @@ data.mbaiki$MAT <- clim$MAT
data.mbaiki$MAP <- clim$MAP
###################### PLOT SELECTION FOR THE ANALYSIS - NEEDS REDOING
vec.basic.var <- c("obs.id","tree.id", "sp", "sp.name", "cluster","plot", "ecocode", "D", "G","BA.G", "dead",
vec.basic.var <- c("obs.id","tree.id", "sp", "sp.name", "cluster","plot", "ecocode", "D", "G","BA.G","year", "dead",
'Lon','Lat',"x", "y", "census",'MAT','MAP')
data.tree <- subset(data.mbaiki, select = c(vec.basic.var))
......
......@@ -186,7 +186,7 @@ data.nsw$MAT <- clim$MAT
data.nsw$MAP <- clim$MAP
vec.abio.var.names <- c("MAT", "MAP")
vec.basic.var <- c("obs.id","tree.id", "sp","sp.name", "cluster", "plot", "ecocode", "D", "G","BA.G", "dead",
vec.basic.var <- c("obs.id","tree.id", "sp","sp.name", "cluster", "plot", "ecocode", "D", "G","BA.G","year", "dead",
"Lon", "Lat","weights","census")
data.tree <- subset(data.nsw, select = c(vec.basic.var,vec.abio.var.names ))
## select only tree above 10cm DBH
......
......@@ -88,7 +88,7 @@ data.nz <- merge(data.nz, data.frame(plot = names(perc.dead), perc.dead = perc.d
###### REMOVE DEAD TREE at first census
data.nz <- subset(data.nz,subset=!is.na(data.nz[["D"]]))
vec.abio.var.names <- c("MAT", "MAP")
vec.basic.var <- c("obs.id","tree.id", "sp","sp.name", "cluster", "plot", "ecocode", "D", "G","BA.G", "dead",
vec.basic.var <- c("obs.id","tree.id", "sp","sp.name", "cluster", "plot", "ecocode", "D", "G","BA.G","year", "dead",
"Lon", "Lat","weights","census")
data.tree <- subset(data.nz, select = c(vec.basic.var, vec.abio.var.names))
## select only tree above 10cm DBH
......
......@@ -109,7 +109,7 @@ data.paracou$MAT <- clim$MAT
data.paracou$MAP <- clim$MAP
###################### PLOT SELECTION FOR THE ANALYSIS - NEEDS REDOING
vec.basic.var <- c("obs.id","tree.id", "sp", "sp.name", "cluster","plot", "ecocode", "D", "G","BA.G", "dead",
vec.basic.var <- c("obs.id","tree.id", "sp", "sp.name", "cluster","plot", "ecocode", "D", "G","BA.G","year", "dead",
'Lon','Lat',"x", "y", "census",'MAT','MAP')
data.tree <- subset(data.paracou, select = c(vec.basic.var))
......
......@@ -120,7 +120,7 @@ colnames(data.spain)[colnames(data.spain) %in% c("mat", "pp", "PET")] <- c("MAT"
data.spain$ecocode <- NULL
colnames(data.spain)[names(data.spain) == "ecocode.merged"] <- c("ecocode")
vec.abio.var.names <- c("MAT", "MAP")
vec.basic.var <- c("obs.id","tree.id", "sp","sp.name", "cluster", "plot", "ecocode", "D", "G","BA.G", "dead",
vec.basic.var <- c("obs.id","tree.id", "sp","sp.name", "cluster", "plot", "ecocode", "D", "G","BA.G","year", "dead",
"Lon", "Lat","weights","census")
data.tree <- subset(data.spain, select = c(vec.basic.var, vec.abio.var.names))
## select only tree above 10cm DBH
......
......@@ -9,9 +9,9 @@ dir.create("output/formatted/Sweden", recursive=TRUE,showWarnings=FALSE)
######################### READ DATA read individuals tree data
data.swe <- read.table("data/raw/Sweden/Swe_NFI_all.txt",header=T,stringsAsFactors=F,sep="\t")
#head(data.swe)
data.swe$tree.id <- apply(cbind(data.swe[["Year"]],data.swe[["TractNr"]],data.swe[["PlotNr"]],
data.swe$tree.id <- apply(cbind(data.swe[["year"]],data.swe[["TractNr"]],data.swe[["PlotNr"]],
data.swe[["TreeNr"]]),1,paste,collapse="_")
data.swe$plot <- apply(cbind(data.swe[["Year"]],data.swe[["TractNr"]],data.swe[["PlotNr"]]),
data.swe$plot <- apply(cbind(data.swe[["year"]],data.swe[["TractNr"]],data.swe[["PlotNr"]]),
1,paste,collapse="_")
dim(data.swe)
#table(table(data.swe$TreeID))
......@@ -128,7 +128,7 @@ data.swe <- merge(data.swe,data.frame(plot=names(perc.dead),perc.dead=perc.dead)
## vec.abio.var.names <- c("MAT", "MAP") ## TODO NO MAT MAP NEED TO LOAD FROM WORLDCLIM
vec.abio.var.names <- c("MAT", "MAP")
vec.basic.var <- c("obs.id","tree.id", "sp","sp.name", "cluster", "plot", "ecocode", "D", "G","BA.G", "dead",
vec.basic.var <- c("obs.id","tree.id", "sp","sp.name", "cluster", "plot", "ecocode", "D", "G","BA.G","year", "dead",
"Lon", "Lat","weights","census")
data.tree <- subset(data.swe, select = c(vec.basic.var, vec.abio.var.names))
## select only tree above 10cm DBH
......
......@@ -104,7 +104,7 @@ data.swiss <- merge(data.swiss, data.frame(plot = data.climate$CLNR, swb = data.
rm(data.climate)
### select good columns
vec.abio.var.names <- c("MAT", "MAP")
vec.basic.var <- c("obs.id","tree.id", "sp","sp.name", "cluster", "plot", "ecocode", "D", "G","BA.G", "dead",
vec.basic.var <- c("obs.id","tree.id", "sp","sp.name", "cluster", "plot", "ecocode", "D", "G","BA.G","year", "dead",
"Lon", "Lat","weights","census")
data.tree <- subset(data.swiss, select = c(vec.basic.var, vec.abio.var.names))
## select only tree above 10cm DBH
......
......@@ -120,7 +120,8 @@ rm(greco, perc.dead, tab.small.div, sel.small.div)
## variables to keep
vec.abio.var.names <- c("MAT", "MAP")
vec.basic.var <- c("obs.id","tree.id", "sp","sp.name", "cluster", "plot", "ecocode", "D", "G","BA.G", "dead",
vec.basic.var <- c("obs.id","tree.id", "sp","sp.name", "cluster", "plot", "ecocode",
"D", "G","BA.G","year", "dead",
"Lon", "Lat","weights","census")
data.tree <- subset(data.us, select = c(vec.basic.var, vec.abio.var.names))
## select only tree above 10cm DBH
......
......@@ -447,6 +447,25 @@ return(accept)
## function compute FD per plot
## TODO
fun.compute.sd.var.cluster <- function(data,var){
cluster.unique.id <- paste(data[['set'],data[['ecocode']],data[['cluster']])
tapply(data[[var]],INDEX=cluster.unique.id,FUN=sd,na.rm=TRUE)
cluster.unique.id <- paste(data[['set']],data[['ecocode']],data[['cluster']])
return(tapply(data[[var]],INDEX=cluster.unique.id,FUN=sd,na.rm=TRUE))
}
fun.compute.mean.var.cluster <- function(data,var){
cluster.unique.id <- paste(data[['set']],data[['ecocode']],data[['cluster']])
return(tapply(data[[var]],INDEX=cluster.unique.id,FUN=mean,na.rm=TRUE))
}
fun.compute.perc.var.cluster <- function(data,var){
cluster.unique.id <- paste(data[['set']],data[['ecocode']],data[['cluster']])
return(tapply(data[[var]],INDEX=cluster.unique.id,FUN=function(x) sum(x,na.rm=TRUE)/length(x)))
}
fun.compute.all.var.cluster <- function(data){
trait.name <- c("Leaf.N","Seed.mass","SLA","Wood.density","Max.height")
var.for.mean <-c('MAT','MAP',"BATOT",paste(trait.name,"abs.CWM.fill",sep="."))
var.for.sd <-c(paste(trait.name,"focal",sep="."))
var.for.per <- c()
}
### TODO LOOK AT FD R cran package to see diversity index that can be computed
......@@ -35,15 +35,25 @@ mat.perc.I <- do.call("rbind", lmat.perc.I)
lmat.perc.B <-mclapply(sets.B,fun.test.set.tree.CWM.B,filedir=filedir,mc.cores= getOption("mc.cores", 5))
mat.perc.B <- do.call("rbind",lmat.perc.B)
rm(lmat.perc.B,lmat.perc.I)
mat.perc <- data.frame(rbind(mat.perc.I,mat.perc.B),stringsAsFactors=FALSE)
rm(mat.perc.B,mat.perc.I)
mat.perc <- data.frame(lapply(mat.perc, function(x) (unlist(x))))
write.csv(mat.perc,file=file.path(filedir, "all.sites.perc.traits.csv"),
row.names=FALSE)
###############################################################
## format table for report
## genus perc
mat.num.g <- mat.perc[,c(1:3,9:13)]
mat.num.g[,4:8] <- mat.perc[,9:13]/mat.perc[,3]
mat.perc <- read.csv(file.path(filedir, "all.sites.perc.traits.csv"))
mat.num.g <- mat.perc[,c('set','ecocode','num.obs',"Leaf.N.num.genus",
"Seed.mass.num.genus","SLA.num.genus",
"Wood.density.num.genus","Max.height.num.genus")]
mat.num.g[,4:8] <- mat.perc[, c("Leaf.N.num.genus","Seed.mass.num.genus",
"SLA.num.genus","Wood.density.num.genus",
"Max.height.num.genus")]/mat.perc[,"num.obs"]
names(mat.num.g) <- c('set','ecoregion','P obs total','P Leaf N',
'P Seed mass','P SLA','P Wood density','P Max height')
......@@ -54,24 +64,42 @@ pandoc.table(mat.num.g,
split.tables='Inf')
## species perc
mat.num.sp <- mat.perc[,c(1:8)]
mat.num.sp[,4:8] <- mat.perc[,4:8]/mat.perc[,3]
mat.num.sp <- mat.perc[,c('set','ecocode','num.obs',"Leaf.N.num.species",
"Seed.mass.num.species","SLA.num.species",
"Wood.density.num.species","Max.height.num.species")]
mat.num.sp[,4:8] <- mat.perc[,4:8]/mat.perc[,'num.obs']
names(mat.num.sp) <- c('set','ecoregion','P obs total','P Leaf N',
'P Seed mass','P SLA','P Wood density','P Max height')
pandoc.table(mat.num.sp,
caption="Number of tree radial growth observation per data sets and ecoregion.",
split.tables='Inf')
### read all data
data.all <- read.csv(file=file.path(filedir, "data.all.csv"))
data.all.sample <- read.csv(file=file.path(filedir, "data.all.csv"),
stringsAsFactors=FALSE,nrows=1000)
classes <- sapply(data.all.sample,class)
classes[classes=='integer'] <- "numeric"
nrows <- as.numeric(system('wc -l < output/processed/data.all.csv',intern=TRUE))
data.all <- read.csv(file=file.path(filedir, "data.all.csv"),
stringsAsFactors=FALSE,
nrows=nrows,
colClasses=classes)
if(dim(data.all)[1] != sum(mat.perc[['num.obs']])) stop('error not same dimension per ecoregion and total')
## data.all[['G']][!(trim.positive.growth(data.all[['G']]) &
## trim.negative.growth(dbh1=data.all[['D']]*10,
## dbh2=data.all[['D']]*10 +data.all[['year']]*data.all[['G']]))] <- NA
## data.all[['BA.G']][!(trim.positive.growth(data.all[['G']]) &
## trim.negative.growth(data.all[['D']]*10,dbh2=data.all[['D']]*10 +data.all[['year']]*data.all[['G']]))] <- NA
data.all[['G']][!(trim.positive.growth(data.all[['G']]) &
trim.negative.growth(dbh1=data.all[['D']]*10,
dbh2=data.all[['D']]*10 +data.all[['year']]*data.all[['G']]))] <- NA
(data.all[['G']] > -50 & !is.na(data.all[['G']])))] <- NA
data.all[['BA.G']][!(trim.positive.growth(data.all[['G']]) &
trim.negative.growth(data.all[['D']]*10,dbh2=data.all[['D']]*10 +data.all[['year']]*data.all[['G']]))] <- NA
(data.all[['G']] > -50 & !is.na(data.all[['G']])))] <- NA
## plot
......@@ -92,6 +120,8 @@ to.dev(fun.plot.xy.set(data.all,var.x='D',var.y='BA.G',cex=0.6),dev=png,
filename="figs/test.processed/fig.xy.D.BA.G.set.png")
to.dev(fun.plot.xy.set(data.all,var.x='D',var.y='G',cex=0.6),dev=png,
filename="figs/test.processed/fig.xy.D.G.set.png")
to.dev(fun.plot.xy.set(data.all,var.x='BATOT',var.y='G',cex=0.6),dev=png,
filename="figs/test.processed/fig.xy.BATOT.G.set.png")
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment