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

first version of quick report in knitr

parent 82a1a01f
......@@ -77,7 +77,7 @@ logexp <- function(exposure = 1)
glmer.output <- glmer(formula = formula, data = df.glmer,
family = binomial(link = logexp(exp(df.glmer$logyear))),
control=glmerControl(optimizer = "bobyqa",
tolPwrss = 1e-3,
tolPwrss = 1e-1,
optCtrl = list(maxfun=500)))
end <- Sys.time()
print(end - Start)
......
......@@ -817,7 +817,7 @@ plot.param.BLUP <- function(list.res,
'Compet x trait dissimilarity'),
col.vec,
pch.vec,
names.bio){
names.bio, ...){
m <- matrix(c(1:4), 1, 4)
layout(m, widths=c(7, 4, 4, 3 )/10,
heights=c(4)/10)
......@@ -838,7 +838,7 @@ for (i in traits){
plot(param.mean, 1:length(param.vec),
yaxt = 'n', xlab = paste('Effect size '), ylab = NA,
xlim = c(-0.45, 0.25) , pch = 16 , cex = 2, cex.lab = 1.5)
, pch = 16 , cex = 2, cex.lab = 1.5, ...)
mtext(i, side=3, cex =1.2)
box(lwd= 2)
abline(v = 0)
......@@ -878,7 +878,8 @@ plot.param.BLUP <- function(list.res,
col.names = c('#e41a1c', '#377eb8', '#4daf4a', '#984ea3', '#ff7f00') ,
col.vec,
pch.vec,
names.bio){
names.bio,
...){
m <- matrix(c(1:4), 1, 4)
layout(m, widths=c(5.3, 2, 2 ,2)/10,
heights=c(4)/10)
......@@ -898,8 +899,8 @@ for (i in traits){
}
plot(param.mean[param.print], (1:length(param.vec))[param.print],
yaxt = 'n', xlab = paste('Effect size '), ylab = NA,
xlim = c(-0.45, 0.25) , pch = 16 , cex = 2, cex.lab = 1.5,
ylim = range(1:length(param.vec)))
pch = 16 , cex = 2, cex.lab = 1.5,
ylim = range(1:length(param.vec)), ...)
mtext(i, side=3, cex =1.2)
box(lwd= 2)
abline(v = 0)
......@@ -1000,7 +1001,8 @@ plot.param <- function(list.res,
col.names = c('#e41a1c', '#377eb8', '#4daf4a', '#984ea3', '#ff7f00') ,
col.vec,
pch.vec,
names.bio){
names.bio,
...){
m <- matrix(c(1:3), 1, 3)
layout(m, widths=c(5.3, 2.39, 2.39 )/10,
heights=c(4)/10)
......@@ -1020,8 +1022,8 @@ for (i in traits){
}
plot(param.mean[param.print], (1:length(param.vec))[param.print],
yaxt = 'n', xlab = paste('Effect size '), ylab = NA,
xlim = c(-0.45, 0.25) , pch = 16 , cex = 2, cex.lab = 1.5,
ylim = range(1:length(param.vec)))
pch = 16, cex = 2, cex.lab = 1.5,
ylim = range(1:length(param.vec)), ...)
mtext(i, side=3, cex =1.2)
box(lwd= 2)
abline(v = 0)
......@@ -1039,11 +1041,26 @@ for (i in traits){
extract.param <- function(trait, list.res,
model = 'lmer.LOGLIN.ER.AD.Tf.r.biomes.species',
param.vec = c("Tf","sumBn", "sumTnBn",
"sumTfBn", "sumTnTfBn.abs")){
list.temp <- list.res[[paste("all.no.log_", trait ,
"_species_", model,
sep = '')]]$lmer.summary
param.mean <- list.temp$fixed.coeff.E[param.vec]
return(param.mean)
}
plot.param.biomes <- function(list.res,
biomes.mat ,
biomes = c("Temperate rain forest",
"Tropical rain forest", "Tropical forest savanna",
"Boreal forest", "Woodland shrubland",
"Tropical rain forest",
"Tropical forest savanna",
"Boreal forest",
"Woodland shrubland",
"Temperate forest") ,
traits = c('Wood.density', 'SLA', 'Max.height'),
param.vec = c("Tf", "sumTnBn",
......@@ -1056,7 +1073,8 @@ plot.param.biomes <- function(list.res,
col.names = c('#e41a1c', '#377eb8', '#4daf4a', '#984ea3', '#ff7f00') ,
col.vec,
pch.vec,
names.bio){
names.bio,
...){
m <- matrix(c(1:4), 1, 4)
layout(m, widths=c(5.3, 2, 2 ,2)/10,
heights=c(4)/10)
......@@ -1078,12 +1096,12 @@ for (i in traits){
if(biome == biomes[1]){
plot(param.mean[param.print], seq.jitter[biomes == biome]+(1:length(param.vec))[param.print],
yaxt = 'n', xlab = paste('Effect size '), ylab = NA,
xlim = c(-0.45, 0.25) , pch = 16 , cex = 2, cex.lab = 1.5,
, pch = pch.vec[biome] , cex = 2, cex.lab = 1.5,
ylim = range(1:length(param.vec)),
col = col.vec[biome])}
col = col.vec[biome], ...)}
if(!biome == biomes[1]){
points(param.mean[param.print], seq.jitter[biomes == biome] + (1:length(param.vec))[param.print],
col = col.vec[biome], , pch = 16 , cex = 2)}
col = col.vec[biome], , pch = pch.vec[biome] , cex = 2)}
mtext(i, side=3, cex =1.2)
box(lwd= 2)
......@@ -1101,7 +1119,7 @@ for (i in traits){
}
par(mar=c(0, 0, 0, 0))
plot(1,1,type="n", axes=F, xlab="", ylab="")
legend("center",legend = biomes, col = col.vec[biomes], lty =1, lwd = 2,
legend("center",legend = biomes, col = col.vec[biomes],pch = pch.vec[biomes], lty =1, lwd = 2,
bty = 'n', cex = 1)
}
......@@ -1141,12 +1159,12 @@ for (i in traits){
if(biome == biomes[1]){
plot(param.mean[param.print], (1:length(param.vec))[param.print]+seq.jitter[biomes == biome],
yaxt = 'n', xlab = paste('Effect size '), ylab = NA,
xlim = c(-0.02, 0.15) , pch = 16 , cex = 2, cex.lab = 1.5,
xlim = c(-0.02, 0.15) , pch = pch.vec[biome] , cex = 2, cex.lab = 1.5,
ylim = range(1:length(param.vec)),
col = col.vec[biome])}
if(!biome == biomes[1]){
points(param.mean[param.print], (1:length(param.vec))[param.print]+seq.jitter[biomes == biome],
col = col.vec[biome], , pch = 16 , cex = 2)}
col = col.vec[biome], , pch = pch.vec[biome] , cex = 2)}
mtext(i, side=3, cex =1.2)
box(lwd= 2)
......
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.CAT.r.biomes.species",
var.BLUP = 'biomes.id',
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|plot.id)+Tf.A_EV+Tf.A_D+Tf.C+logD+sumBn+sumTfBn.A_EV+sumTfBn.A_D+sumTfBn.C+sumTnBn.A_EV+sumTnBn.A_D+sumTnBn.C+sumTnTfBn.abs+(1|species.id) +(logD-1|species.id) +(sumBn-1|species.id)+(Tf-1|biomes.id)+(sumBn-1|biomes.id)+(sumTfBn.A_EV-1|biomes.id)+(sumTfBn.A_D-1|biomes.id)+(sumTfBn.C-1|biomes.id)+(sumTnBn.A_EV-1|biomes.id)+(sumTnBn.A_D-1|biomes.id)+(sumTnBn.C-1|biomes.id)+(sumTnTfBn.abs-1|biomes.id)"))
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|plot.id)+Tf.A_EV+Tf.A_D+Tf.C+logD+sumBn+sumTfBn.A_EV+sumTfBn.A_D+sumTfBn.C+sumTnBn.A_EV+sumTnBn.A_D+sumTnBn.C+sumTnTfBn.abs+(1|species.id) +(logD-1|species.id) +(sumBn-1|species.id)+(Tf.A_EV-1|biomes.id)+(Tf.A_D-1|biomes.id)+(Tf.C-1|biomes.id)+(sumBn-1|biomes.id)+(sumTfBn.A_EV-1|biomes.id)+(sumTfBn.A_D-1|biomes.id)+(sumTfBn.C-1|biomes.id)+(sumTnBn.A_EV-1|biomes.id)+(sumTnBn.A_D-1|biomes.id)+(sumTnBn.C-1|biomes.id)+(sumTnTfBn.abs-1|biomes.id)"))
}
......
......@@ -172,73 +172,25 @@ boxplot(y~x.cut,at=mid.points,outline=FALSE,names=round(mid.points,0),
}
################
### PLOTS function for biomes
fun.poly.obj <- function(){
library(sp)
biomes.data <- read.csv('R/utils/biomes.csv', stringsAsFactors = FALSE)
list.coord <- lapply(unique(biomes.data$biome),
function(biome.id,data)
cbind(data$y[data$biome == biome.id] * 10,
data$x[data$biome == biome.id]),
data=biomes.data)
names(list.coord) <- unique(biomes.data$biome)
list.Polygon <- lapply(list.coord, Polygon)
list.Polygons <- lapply(1:length(list.Polygon),
function(i, x) {Polygons(list(x[[i]]), names(x)[i])},
x = list.Polygon)
poly.biomes <- SpatialPolygons(list.Polygons)
DF <- data.frame(biomes= names(poly.biomes), row.names = names(poly.biomes))
poly.DF <- SpatialPolygonsDataFrame(poly.biomes, DF)
return(list(poly.DF = poly.DF))
}
## function to overlay points on biomes
fun.overly.plot.on.biomes <- function(MAP, MAT, names.vec,
poly = fun.poly.obj()$poly.DF){
require(BIOMEplot)
xy = cbind(MAP, MAT)
dimnames(xy)[[1]] = names.vec
pts = SpatialPoints(xy)
return(over(pts, poly))
}
plot.biome.map <- function(poly = fun.poly.obj()$poly.DF, add.legend = TRUE){
require(sp)
# MAKE EMPTY PLOT
par(mar = c(5.1, 3.5, 0.1, 0.1),
mgp
= c(1.8, 0.6,0))
plot(0,0, type = 'n', xlim = c(0, 450), ylim = c(30, -15),
xlab = "Mean annual precipitation (cm)",
ylab = expression(paste('Mean annual temperature ', (~degree~C))))
# MAKE POLYGONS
plot(poly,col = make.transparent(c("navajowhite3", "darkgoldenrod1",
"sienna","darkolivegreen4",
"darkseagreen3", "forestgreen",
"darkgreen","olivedrab","gray"),
0.9),
border = NA,
add = TRUE)
if(add.legend) {
legend("topright", legend = poly$biomes,
fill =make.transparent(c("navajowhite3", "darkgoldenrod1",
"sienna","darkolivegreen4",
"darkseagreen3", "forestgreen",
"darkgreen","olivedrab","gray"),
0.9),
bty = 'n')}
}
plot.points.on.biome.map.diff.tropical <- function(MAP, MAT, set,ecocode,
cols = niceColors(),
add.legend=TRUE){
plot.biome.map()
require(BIOMEplot)
plot_biome()
col <- c()
pch <- c()
cex <- c()
......@@ -264,7 +216,8 @@ plot.points.on.biome.map.diff.tropical <- function(MAP, MAT, set,ecocode,
plot.points.on.biome.map.ecocode <- function(MAP, MAT, MAP.sd, MAT.sd, set,
ecocode, cols = niceColors(),
add.legend = TRUE){
plot.biome.map()
require(BIOMEplot)
plot_biome()
col <- c()
pch <- c()
cex <- c()
......@@ -286,7 +239,8 @@ plot.points.on.biome.map.ecocode <- function(MAP, MAT, MAP.sd, MAT.sd, set,
plot.points.on.biome.map <- function(MAP, MAT, set,
cols = fun.col.pch.set()$col.vec){
plot.biome.map()
require(BIOMEplot)
plot_biome()
points(MAP, MAT, col = cols[set])
}
......@@ -309,7 +263,8 @@ plot.points.on.biome.map(MAP[j], MAT[j], set[j], cols = cols[sets[i]])
plot.ellips.on.biome.map <- function(MAP, MAT, set,
cols = fun.col.pch.set()$col.vec,
lty.vec = fun.col.pch.set()$pch.vec){
plot.biome.map()
require(BIOMEplot)
plot_biome()
lapply(unique(set), fun.ellipsoid.hull,
x = MAP, y = MAT, sets.v = set,
col.vec = cols, lty.vec = lty.vec)
......
D1 := ../../figs
figs:= biome_ellipse.pdf
all: quick.report.pdf quick.report.html figs
all: quick.report.pdf quick.report.html
quick.report.md: quick.report.R
Rscript -e "library(sowsear); sowsear('quick.report.R', 'Rmd')"
Rscript -e "library(knitr); knit('quick.report.Rmd', output = 'quick.report.md')"
figure/biome_ellipse.pdf: $(D1)/biome_ellipse.pdf
cp $< $@
quick.report.pdf: quick.report.md figure/biome_ellipse.pdf
pandoc $< --standalone --template=include.tex --variable mainfont="Times New Roman" --variable sansfont=Arial --variable fontsize=12pt --latex-engine=xelatex -o $@
......
(TeX-add-style-hook "include"
(lambda ()
(LaTeX-add-bibliographies
"$biblio-files$")
(TeX-add-symbols
"euro"
"maxwidth"
"Oldincludegraphics")
(TeX-run-style-hooks
"babel"
"$lang$"
"polyglossia"
"ulem"
"normalem"
"unicode=true"
"hyperref"
"xetex"
"unicode=false"
"setpagesize=false"
"graphicx"
"booktabs"
"longtable"
"fancyvrb"
"listings"
"biblatex"
"natbib"
"$endfor$"
"$for(geometry)$$geometry$$sep$"
"fancyhdr"
"xunicode"
"xltxtra"
"fontspec"
"eurosym"
"inputenc"
"utf8"
"upquote"
"microtype"
"fixltx2e"
"ifluatex"
"ifxetex"
"amsmath"
"amssymb"
"lmodern"
""
"fontenc"
"T1"
"geometry"
"vmargin=2cm"
"hmargin=2cm"
"latex2e"
"$documentclass$10"
"$documentclass$"
"a4paper"
"$endif$"
"$endif$$if(lang)$$lang$"
"$if(fontsize)$$fontsize$")))
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -77,7 +77,7 @@ $if(verbatim-in-note)$
\usepackage{fancyvrb}
$endif$
$if(tables)$
\usepackage{longtable}
\usepackage{longtable,booktabs}
$endif$
$if(graphics)$
\usepackage{graphicx}
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -22,20 +22,20 @@ for trait in "'SLA'" "'Wood.density'" "'Max.height'"; do
# echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.multiple.model.for.set.one.trait(model.files.glmer.Tf.1, run.glmer,$trait,type.filling='species', sample.size = $samplesize, var.sample = 'ecocode.id');print('done')\"" > trait.workshop/glmerspecies1${trait}.sh
# qsub trait.workshop/glmerspecies1${trait}.sh -l nodes=1:ppn=1,mem=8gb -N "glmerall1${trait}" -q opt32G -j oe
# echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.multiple.model.for.set.one.trait(c(model.files.glmer.Tf.2), run.glmer,$trait,type.filling='species', sample.size = $samplesize, var.sample = 'ecocode.id');print('done')\"" > trait.workshop/glmerspecies2${trait}.sh
# qsub trait.workshop/glmerspecies2${trait}.sh -l nodes=1:ppn=1,mem=8gb -N "glmerall2${trait}" -q opt32G -j oe
echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.multiple.model.for.set.one.trait(c(model.files.glmer.Tf.2), run.glmer,$trait,type.filling='species', sample.size = $samplesize, var.sample = 'ecocode.id');print('done')\"" > trait.workshop/glmerspecies2${trait}.sh
qsub trait.workshop/glmerspecies2${trait}.sh -l nodes=1:ppn=1,mem=8gb -N "glmerall2${trait}" -q opt32G -j oe
echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.multiple.model.for.set.one.trait(c(model.files.glmer.Tf.3), run.glmer,$trait,type.filling='species', sample.size = $samplesize, var.sample = 'ecocode.id');print('done')\"" > trait.workshop/glmerspecies3${trait}.sh
qsub trait.workshop/glmerspecies3${trait}.sh -l nodes=1:ppn=1,mem=8gb -N "glmerall3${trait}" -q opt32G -j oe
# echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.multiple.model.for.set.one.trait(c(model.files.glmer.Tf.3), run.glmer,$trait,type.filling='species', sample.size = $samplesize, var.sample = 'ecocode.id');print('done')\"" > trait.workshop/glmerspecies3${trait}.sh
# qsub trait.workshop/glmerspecies3${trait}.sh -l nodes=1:ppn=1,mem=8gb -N "glmerall3${trait}" -q opt32G -j oe
# echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.multiple.model.for.set.one.trait(model.files.glmer.Tf.1, run.glmer,$trait,type.filling='species', sample.size = $samplesize, var.sample = 'species.id');print('done')\"" > trait.workshop/glmerspeciesS1${trait}.sh
# qsub trait.workshop/glmerspeciesS1${trait}.sh -l nodes=1:ppn=1,mem=8gb -N "glmerallS1${trait}" -q opt32G -j oe
# echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.multiple.model.for.set.one.trait(c(model.files.glmer.Tf.2), run.glmer,$trait,type.filling='species', sample.size = $samplesize, var.sample = 'species.id');print('done')\"" > trait.workshop/glmerspeciesS2${trait}.sh
# qsub trait.workshop/glmerspeciesS2${trait}.sh -l nodes=1:ppn=1,mem=8gb -N "glmerallS2${trait}" -q opt32G -j oe
echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.multiple.model.for.set.one.trait(c(model.files.glmer.Tf.2), run.glmer,$trait,type.filling='species', sample.size = $samplesize, var.sample = 'species.id');print('done')\"" > trait.workshop/glmerspeciesS2${trait}.sh
qsub trait.workshop/glmerspeciesS2${trait}.sh -l nodes=1:ppn=1,mem=8gb -N "glmerallS2${trait}" -q opt32G -j oe
echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.multiple.model.for.set.one.trait(c(model.files.glmer.Tf.3), run.glmer,$trait,type.filling='species', sample.size = $samplesize, var.sample = 'species.id');print('done')\"" > trait.workshop/glmerspeciesS3${trait}.sh
qsub trait.workshop/glmerspeciesS3${trait}.sh -l nodes=1:ppn=1,mem=8gb -N "glmerallS3${trait}" -q opt32G -j oe
# echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.multiple.model.for.set.one.trait(c(model.files.glmer.Tf.3), run.glmer,$trait,type.filling='species', sample.size = $samplesize, var.sample = 'species.id');print('done')\"" > trait.workshop/glmerspeciesS3${trait}.sh
# qsub trait.workshop/glmerspeciesS3${trait}.sh -l nodes=1:ppn=1,mem=8gb -N "glmerallS3${trait}" -q opt32G -j oe
......
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