Commit f4bf5edf authored by Georges Kunstler's avatar Georges Kunstler
Browse files

add convergence test from bolker

parent 914b9131
......@@ -41,7 +41,11 @@ fun.call.lmer.and.save <- function(formula, df.lmer, path.out,
var.sample, select.biome){
lmer.output <- lmer(formula = formula, data = df.lmer, REML = FALSE,
control = lmerControl(optCtrl = list(maxfun = 40000) ) )
print(summary(lmer.output))
print(warnigs())
print(summary(lmer.output))
relgrad <- with(lmer.output@optinfo$derivs,solve(Hessian,gradient))
print('test convergence of relative gradient of hessian')
print(max(abs(relgrad)) )
if(is.na(select.biome)){
name.file <- paste(var.sample,
"results.nolog.all.rds", sep ='.')
......
......@@ -93,15 +93,39 @@ cols.vec <- fun.col.pch.set()$col.vec
cols.vec <- rep('black', length(cols.vec))
names(cols.vec) <- names(fun.col.pch.set()$col.vec)
plot.ellips.on.biome.map(MAP = clim$MAP/10,
MAT = clim$MAT,
set = clim$set, cols = cols.vec)
MAT = clim$MAT,
set = clim$set,
cols = cols.vec)
##+ Fig1b, fig.cap="**Number of tree per Whittaker biomes**", fig.width=7, fig.height=7, echo = FALSE, results = 'hide', warning=FALSE, message=FALSE
data.all <- readRDS.root('output/processed/data.all.no.log.rds')
data.all <- fun.merge.biomes(data.all)
par(mar = c(3, 13, 3,3))
barplot(table(data.all$biomes), horiz = TRUE, las = 2)
rm(data.all)
data.SLA.ecocode <- load.data.for.lmer(trait = 'SLA',
cat.TF = FALSE,
fname = 'data.all.no.log.rds',
base.dir = "output/processed",
sample.size. = 1000000,
var.sample. = 'ecocode')
data.SLA.species <- load.data.for.lmer(trait = 'SLA',
cat.TF = FALSE,
fname = 'data.all.no.log.rds',
base.dir = "output/processed",
sample.size. = 1000000,
var.sample. = 'sp.name')
data.SLA.NA <- load.data.for.lmer(trait = 'SLA',
cat.TF = FALSE,
fname = 'data.all.no.log.rds',
base.dir = "output/processed",
sample.size. = 1000000,
var.sample. = NA)
## \newpage
## #Results
......
This diff is collapsed.
......@@ -22,7 +22,7 @@ samplesize=$1
#
for trait in "'SLA'" "'Wood.density'" "'Max.height'" "'Seed.mass'" "'Leaf.N'"; do
echo "/usr/local/R/R-3.1.1/bin/Rscript -e \"source('R/analysis/lmer.run.R'); run.multiple.model.for.set.one.trait(c(model.files.lmer.Tf.1, model.files.lmer.Tf.2), run.lmer,$trait,type.filling='species', sample.size = $samplesize, var.sample = 'ecocode');print('done')\"" > Rscript_temp/species2${trait}.sh
echo "/usr/local/R/R-3.1.1/bin/Rscript -e \"source('R/analysis/lmer.run.R'); run.multiple.model.for.set.one.trait(c(model.files.lmer.Tf.1, model.files.lmer.Tf.2), run.lmer,$trait,type.filling='species', sample.size = $samplesize, var.sample = 'biomes');print('done')\"" > Rscript_temp/species2${trait}.sh
qsub Rscript_temp/species2${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2${trait}" -q opt32G -j oe
echo "/usr/local/R/R-3.1.1/bin/Rscript -e \"source('R/analysis/lmer.run.R'); run.multiple.model.for.set.one.trait(c(model.files.lmer.Tf.1, model.files.lmer.Tf.2), run.lmer,$trait,type.filling='species', sample.size = NA, var.sample = NA);print('done')\"" > Rscript_temp/species2T${trait}.sh
......
......@@ -36,7 +36,10 @@ lapply(files, FUN = download.file2 ,
#unzip data
fun.unzip <- function(file, data_path){
unzip(file.path(data_path, paste0(file, ".zip")), exdir = file.path(data_path, file))
unzip(file.path(data_path, paste0(file, ".zip")), exdir = file.path(data_path, file))
obj <- readOGR(file.path(data_path, file), file, encoding = "latin1")
saveRDS(obj, file.path(data_path, paste0(file, ".rds")))
unlink(file.path(data_path, file), recursive = TRUE)
}
lapply(files, FUN = fun.unzip ,
......@@ -91,8 +94,6 @@ newmap <- getMap(resolution = 'high')
plot(newmap, xlim = c(7, 16), ylim = c(47, 55))
points(data.sp2, cex = 0.2, col = 'red')
## Need to ask Bjoern if good EPSG seems ok
### VARIABLES SELECTION
vars.select1 <- c('bhd_130',# diamter at 130 cm corrected in mm
......@@ -129,5 +130,44 @@ data_germany<- merge(wzp_87_df[, vars.select1], wzp_02_df[, vars.select2],
by = c("tree.id"), all.x = TRUE, all.y = FALSE)
## check growth
data_germany$d_bhd <- data_germany$bhd_130.y - data_germany$bhd_130.x
data_germany$plot.id <- paste(data_germany$tnr.x, data_germany$enr.x, sep = '_')
## check growth
plot(data_germany$bhd_130.x, data_germany$d_bhd)
data_DE <- data_germany[!is.na(data_germany$d_bhd), ]
nrow(data_DE)
## plot with no remeasurment at all
plot.perc.remeasure <- tapply(is.na(data_germany$d_bhd), INDEX = data_germany$plot.id,
FUN = function(x) sum(x)/length(x))
plot.no.remeasure <- names(plot.perc.remeasure)[plot.perc.remeasure==1]
data_remeasure <- data_germany[!data_germany$plot.id %in% plot.no.remeasure, ]
table(is.na(data_remeasure$d_bhd))
## still way to much missing tree to be natural mortality. How to tell apart harvested tree??
# Elevation a.s.l.
# Aufnahmezeitpunkt (Monat/Jahr)
# Nutzungsart (Plantage) und Eingriffe (Störungen/Ernten)
plot(data_germany$bhd_130.x, data_germany$bhd_130.y - data_germany$bhd_130.x)
ecken_02 <- readRDS(file.path(data_path, "bwi2002dt_waldecken.rds"))
ecken_02_df <- as.data.frame(ecken_02)
head(ecken_02_df)
plot(m_hoe ~ m_bhd, data = subset(wzp_02_df, ba_gattung == "Picea" & ba_art == "abies"))
points(m_hoe ~ m_bhd, data = subset(wzp_02_df, ba_gattung == "Pinus" & ba_art == "sylvestris"), col = "red")
points(m_hoe ~ m_bhd, data = subset(wzp_02_df, ba_gattung == "Fagus" & ba_art == "sylvatica"), col = "blue")
points(m_hoe ~ m_bhd, data = subset(wzp_02_df, ba_gattung == "Quercus" & ba_art == "petraea"), col = "green")
points(m_hoe ~ m_bhd, data = subset(wzp_02_df, ba_gattung == "Larix" & ba_art == "decidua"), col = "orange")
plot(data_germany$bhd_130.x, (data_germany$bhd_130.y - data_germany$bhd_130.x)/15)
plot(m_bhd ~ al_ba2002, data = subset(wzp_02_df, ba_gattung == "Picea" & ba_art == "abies"))
plot(m_hoe ~ al_ba2002, data = subset(wzp_02_df, ba_gattung == "Picea" & ba_art == "abies"))
plot(m_hoe ~ al_ba2002, data = subset(wzp_02_df, ba_gattung == "Pinus" & ba_art == "sylvestris"))
points(m_hoe ~ al_ba2002, data = subset(wzp_02_df, ba_gattung == "Pinus" & ba_art == "sylvestris" & bkl_txt == "vorherrschender Baum"), col = "red")
points(m_hoe ~ al_ba2002, data = subset(wzp_02_df, ba_gattung == "Pinus" & ba_art == "sylvestris" & bkl_txt == "beherrschter Baum"), col = "blue", pch = 19)
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