Commit 914b9131 authored by Georges Kunstler's avatar Georges Kunstler
Browse files

progress on analysis lmer fixed and include random biomes in stan working

parent cc6aa5d3
......@@ -987,6 +987,9 @@ for (i in traits){
}
}
plot.param <- function(list.res,
model = 'lmer.LOGLIN.ER.AD.Tf',
traits = c('Wood.density', 'SLA', 'Max.height'),
......@@ -1003,7 +1006,7 @@ plot.param <- function(list.res,
pch.vec,
names.bio,
...){
m <- matrix(c(1:3), 1, 3)
m <- matrix(c(1:3), 1, 3)
layout(m, widths=c(5.3, 2.39, 2.39 )/10,
heights=c(4)/10)
for (i in traits){
......@@ -1186,13 +1189,107 @@ legend("center",legend = biomes, col = col.vec[biomes], lty =1, lwd = 2,
bty = 'n', cex = 1.2)
}
## get fixed biomes
fun.get.fixed.biomes <- function(var, list,
biomes.vec = c("Boreal forest",
"Subtropical desert",
"Temperate forest",
"Temperate grassland desert",
"Temperate rain forest",
"Tropical forest savanna",
"Tropical rain forest",
"Tundra",
"Woodland shrubland")){
param <- list$lmer.summary$fixed.coeff.E
remaining <- grep(paste0(':',var), names(param))
boreal <- seq_len(length(param))[names(param) == var]
select <- c(boreal, remaining)
## create design matrix
ff <- ~ biomes.id
mm <- model.matrix(ff, data.frame(biomes.id = biomes.vec))
# compute mean and std based on Bolker post
param.biomes <- drop(mm %*% param[select])
names(param.biomes) <- sort(unique(data$biomes.id))
std.biomes <- sqrt(diag(mm %*% tcrossprod(list$vcov[select, select],mm)))
return(list(fixed.biomes = param.biomes,
fixed.biomes.std = std.biomes))
}
##
plot.param.biomes.fixed <- function(list.res,
model,
biomes = c("Boreal forest",
"Subtropical desert",
"Temperate forest",
"Temperate grassland desert",
"Temperate rain forest",
"Tropical forest savanna",
"Tropical rain forest",
"Tundra",
"Woodland shrubland") ,
traits = c('Wood.density', 'SLA', 'Max.height'),
param.vec = c("Tf", "sumTnBn",
"sumTfBn", "sumTnTfBn.abs"),
param.names = c('Direct trait',
'Compet effect x trait',
'Compet response x trait',
'Compet x trait dissimilarity'),
param.print = 1:5,
col.names = c('#e41a1c', '#377eb8', '#4daf4a',
'#984ea3', '#ff7f00') ,
col.vec,
pch.vec,
names.bio,
...){
m <- matrix(c(1:4), 1, 4)
layout(m, widths=c(5.3, 2, 2 ,2.5)/10,
heights=c(4)/10)
for (i in traits){
list.temp <- list.res[[paste0("all.no.log_", i ,
"_species_", model)]]
## NEED VAR
for (n.vars in seq_len(length(param.vec))){
list.fixed <- fun.get.fixed.biomes(param.vec[n.vars], list.temp,
biomes.vec = biomes)
param.mean <- list.fixed$fixed.biomes
param.std <- list.fixed$fixed.biomes.std
if(i == traits[1]) {
par(mar=c(6,22,3,0))
}else{
par(mar=c(6,0,3,0))
}
seq.jitter <- seq(-10,10, length.out = length(biomes))/120
if(n.vars == 1){
plot(param.mean, seq.jitter+n.vars,
yaxt = 'n', xlab = 'Effect size ', ylab = NA,
, pch = pch.vec[biomes] , cex = 2, cex.lab = 1.5,
ylim = range(seq_len(length(param.vec))),
col = col.vec[biomes], ...)}
if(n.vars != 1){
points(param.mean, seq.jitter+n.vars,
col = col.vec[biomes], pch = pch.vec[biomes], cex = 2)}
mtext(i, side=3, cex =1.2)
box(lwd= 2)
abline(v = 0)
if(i == traits[1]) {lapply(1:length(param.vec),
fun.axis.one.by.one,
side = 2,
labels = param.names,
cols.vec = col.names)
}
fun.plot.error.bar.horiz(param.mean,
seq.jitter+n.vars,
param.std, col = col.vec[biomes])
}
}
par(mar=c(0, 0, 0, 0))
plot(1,1,type="n", axes=F, xlab="", ylab="")
legend("center",legend = biomes, col = col.vec[biomes],
pch = pch.vec[biomes], lty =1, lwd = 2,
bty = 'n', cex = 1)
}
##
......
......@@ -72,7 +72,7 @@ parameters{
real<lower=0.0001,upper=5> sigma_sumBn_biomes;
real<lower=0.0001,upper=5> sigma_sumTfBn_biomes;
real<lower=0.0001,upper=5> sigma_sumTnBn_biomes;
real<lower=0.0001,upper=5> sigma_sumTnTfBn_biomes;
real<lower=0.0001,upper=5> sigma_sumTnTfBn_abs_biomes;
real<lower=0.0001,upper=5> sigma;
real<lower=-5,upper=5> param_Tf_biomes[N_biomes];
real<lower=-5,upper=5> param_sumBn_biomes[N_biomes];
......
......@@ -11,6 +11,12 @@
source('R/analysis/stan.run.R')
run.multiple.model.for.set.one.trait(model.files.stan.Tf.2, run.stan, 'SLA',
type.filling='species', sample.size = 2000, var.sample = 'ecocode',
iter = 1000, warmup = 500, chains = 3,
init.TF = FALSE, parallel.TF = FALSE)
### TEST simple model on France only
df.lmer <- load.data.for.lmer(trait = 'SLA',
......@@ -22,37 +28,10 @@ source('R/analysis/stan.run.R')
stan.list <- fun.turn.in.list.for.jags.stan(df.lmer,
cat.TF = FALSE)
source('R/analysis/model.stan/model.stan.LOGLIN.size.fixed.R', local = TRUE)
rm(df.lmer)
source('R/analysis/model.stan/model.stan.LOGLIN.ER.AD.Tf.r.biomes.R', local = TRUE)
model <- load.model()
fun.init.stan <- function(chain_id= 1, stan.list){
init <- list(intercept = 0+(chain_id-1)/10,
intercept_species = rnorm(stan.list$N_species,
0,
0.5+(chain_id-1)/10),
intercept_plot = rnorm(stan.list$N_plot,
0,
0.5 +(chain_id-1)/10),
mean_logD = 2/3+(chain_id-1)/10,
sigma_inter_species = 0.5 +(chain_id-1)/10,
sigma_inter_plot = 0.5 + +(chain_id-1)/10,
sigma = 0.5 +(chain_id-1)/10)
}
inits <- list(model$init.fun(1, stan.list),
model$init.fun(2, stan.list),
model$init.fun(3, stan.list))
library(rstan)
system.time( stan.output <- stan(model_code = model$stan,
data = stan.list,
pars = model$pars,
init = inits,
iter = 1000,
warmup = 500,
chains = 3,
verbose = FALSE))
source('R/analysis/model.stan/model.stan.LOGLIN.ER.AD.Tf.fixed.biomes.R', local = TRUE)
model <- load.model()
......
......@@ -199,7 +199,7 @@ pandoc.table(mat.param, caption = "Full effect size")
##+ Fig3, fig.cap='**effect size parameters with data re-sampled per species.**', fig.width=10, fig.height=5, echo = FALSE
plot.param(list.all.results.species.id ,
model = 'lmer.LOGLIN.ER.AD.Tf.fixed.biomes.species',
model = 'lmer.LOGLIN.ER.AD.Tf.fixed.species',
traits = c('Wood.density', 'SLA', 'Max.height'),
param.vec = c("Tf","sumBn", "sumTfBn",
"sumTnBn", "sumTnTfBn.abs"),
......@@ -219,7 +219,7 @@ plot.param(list.all.results.species.id ,
##+ Fig3b, fig.cap='**effect size parameters with no subsampling.**', fig.width=10, fig.height=5, echo = FALSE
plot.param(list.all.results ,
model = 'lmer.LOGLIN.ER.AD.Tf.fixed.biomes.species',
model = 'lmer.LOGLIN.ER.AD.Tf.fixed.species',
traits = c('Wood.density', 'SLA', 'Max.height'),
param.vec = c("Tf","sumBn", "sumTfBn",
"sumTnBn", "sumTnTfBn.abs"),
......@@ -313,6 +313,23 @@ plot.param.biomes(list.all.results.ecocode.id.biomes,
col.vec = fun.col.pch.biomes()$col.vec,
pch.vec = fun.col.pch.biomes()$pch.vec,
xlim = c(-0.45, 0.35) )
x11()
plot.param.biomes.fixed(list.all.results.species.id,
model = 'lmer.LOGLIN.ER.AD.Tf.fixed.biomes.species',
param.vec = c("Tf","sumBn", "sumTfBn",
"sumTnBn", "sumTnTfBn.abs"),
param.names = c(expression(paste("Direct trait effect ",
(m[1] %*% t[f]))),
expression(paste('Comp mean ', (c[0]))),
expression(paste('Compet response ',
(c[r] %*% t[f]))),
expression(paste('Compet effect ',
(c[e] %*% t[n]))),
expression(paste('Limit simil ',
(c[l] %*% abs(t[n] - t[f]))))) ,
col.vec = fun.col.pch.biomes()$col.vec,
pch.vec = fun.col.pch.biomes()$pch.vec,
xlim = c(-0.5, 0.5) )
## The Figure 6 represents only the limiting similarity parameters estimated separately per biomes (re-sampling with weight per ecocode).
......
......@@ -7,7 +7,12 @@ dir.create("output/formatted/Germany", recursive = TRUE, showWarnings = FALSE)
#############
## DOWNLOAD DATA FROM WEBSITE
data_path <- "data/raw/Germany"
if(!"wzp_87.rds" %in% list.files(data_path)){
gdi.web1 <- "https://gdi.ti.bund.de/geoserver/bwi_2009dt/ows?typeName=bwi_2009dt:"
gdi.web2 <- "&service=WFS&version=1.0.0&request=GetFeature&outputFormat=SHAPE-ZIP"
files <- c("bwi2002dt_wzp4", "bwi2002dt_trakte",
......@@ -45,31 +50,34 @@ library(sp)
## "bwi2002dt_wzp4"),
## proj4string=CRS("+init=epsg:5676"))
wzp_02 <- readOGR(file.path(data_path, "bwi2002dt_wzp4"), "bwi2002dt_wzp4",
p4s = "+init=epsg:31466")
wzp_02 <- readOGR(file.path(data_path, "bwi2002dt_wzp4"), "bwi2002dt_wzp4")
wzp_87 <- readOGR(file.path(data_path, "bwi1987dt_wzp4"), "bwi1987dt_wzp4",
p4s = "+init=epsg:31466")
wzp_87 <- readOGR(file.path(data_path, "bwi1987dt_wzp4"), "bwi1987dt_wzp4")
class(wzp_02)
## save as rds to load faster later on
saveRDS(wzp_87, file.path(data_path,'wzp_87.rds'))
saveRDS(wzp_02, file.path(data_path,'wzp_02.rds'))
}else{
wzp_87 <- readRDS(file.path(data_path,'wzp_87.rds'))
wzp_02 <- readRDS(file.path(data_path,'wzp_02.rds'))
}
# get proj
p4.bwi <- proj4string(wzp_87)
# get centroids
centroids_02<- getSpPPolygonsLabptSlots(wzp_02)
centroids_87<- getSpPPolygonsLabptSlots(wzp_87)
colnames(centroids_87) <- colnames(centroids_02) <- c('x', 'y')
wzp_02_df<- cbind(as.data.frame(wzp_02), centroids_02)
wzp_87_df<- cbind(as.data.frame(wzp_87), centroids_87)
####### change coordinates system of x y to be in lat long WGS84
library(sp)
library(dismo)
library(rgdal)
data.sp <- wzp_87_df[, c("x", "y")]
coordinates(data.sp) <- c("x", "y")
proj4string(data.sp) <- CRS("+init=epsg:2166")
proj4string(data.sp) <- p4.bwi
# define projection system of our data
summary(data.sp)
data.sp2 <- spTransform(data.sp, CRS("+init=epsg:4326"))
......@@ -82,7 +90,6 @@ newmap <- getMap(resolution = 'high')
## different resolutions available
plot(newmap, xlim = c(7, 16), ylim = c(47, 55))
points(data.sp2, cex = 0.2, col = 'red')
plot(wzp_87_df$x, wzp_87_df$y, cex = 0.2, col = 'red')
## Need to ask Bjoern if good EPSG seems ok
......@@ -123,4 +130,4 @@ data_germany<- merge(wzp_87_df[, vars.select1], wzp_02_df[, vars.select2],
## check growth
plot(data_germany$bhd_130.x, data_germany$bhd_130.y - data_germany$bhd_130.x)
plot(data_germany$bhd_130.x, (data_germany$bhd_130.y - data_germany$bhd_130.x)/15)
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