Commit 13bfdac4 authored by Georges Kunstler's avatar Georges Kunstler
Browse files

progress on quick report for skype meeting

parent e9c709ce
......@@ -11,9 +11,7 @@ files <- c()
for (trait in traits){
for(model in c(model.files.glmer.Tf.1,
model.files.glmer.Tf.2,
model.files.glmer.Tf.3,
model.files.glmer.Tf.4,
model.files.glmer.Tf.5)){
model.files.glmer.Tf.3)){
source(model, local = TRUE)
model.obj <- load.model()
pathout <- output.dir.glmer('all.no.log', model.obj$name, trait,
......@@ -27,3 +25,27 @@ for (trait in traits){
out <- lapply(files, summarise.glmer.output.all.list)
names(out) <- lapply(lapply(files,files.details.all),function(x) paste(as.vector(x[-5]),collapse="_"))
saveRDS(out,file='output/glmer.global.list.out.ecocode.id.rds')
## LOOP OVER FILES AND NOT SCAN ALL FILES
traits <- c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass")
type.filling <- 'species'
var.sample <- 'species.id'
files <- c()
for (trait in traits){
for(model in c(model.files.glmer.Tf.1,
model.files.glmer.Tf.2,
model.files.glmer.Tf.3)){
source(model, local = TRUE)
model.obj <- load.model()
pathout <- output.dir.glmer('all.no.log', model.obj$name, trait,
type.filling=type.filling)
files <- c(files,file.path(pathout,paste(var.sample, "glmer.results.global.rds")))
}
}
out <- lapply(files, summarise.glmer.output.all.list)
names(out) <- lapply(lapply(files,files.details.all),function(x) paste(as.vector(x[-5]),collapse="_"))
saveRDS(out,file='output/glmer.global.list.out.species.id.rds')
......@@ -44,7 +44,7 @@ model.files.glmer.Tf.4 <-
model.files.glmer.Tf.5 <-
c("R/analysis/model.glmer/model.glmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.ecocode.species.R")
### TODO ER QND QD TO CHQNGE NEED TO ADD CAT
### TODO ER AND AD TO CHAbNGE NEED TO ADD CAT
fun.call.glmer.and.save <- function(formula, df.glmer, path.out, var.sample){
library(MASS)
......@@ -75,7 +75,10 @@ logexp <- function(exposure = 1)
## glmer.output <- glmer(formula = formula, data = df.glmer,
## family = binomial(link = 'cloglog'))
glmer.output <- glmer(formula = formula, data = df.glmer,
family = binomial(link = logexp(exp(df.glmer$logyear))))
family = binomial(link = logexp(exp(df.glmer$logyear))),
control=glmerControl(optimizer = "bobyqa",
tolPwrss = 1e-3,
optCtrl = list(maxfun=500)))
end <- Sys.time()
print(end - Start)
print(summary(glmer.output))
......@@ -258,6 +261,7 @@ return(
plot.id = plot.id,
set.id = set.id,
ecocode.id = ecocode.id,
biomes.id = biomes.id,
sumTnTfBn.diff = fun.standardized.var(sumTnTfBn.diff),
sumTnTfBn.abs = fun.standardized.var(sumTnTfBn.abs),
Tf = fun.standardized.var(data.tree[[tf]]),
......
......@@ -3,30 +3,62 @@
source("R/analysis/lmer.output-fun.R")
source("R/analysis/lmer.run.R")
## LOOP OVER FILES AND NOT SCAN ALL FILES
traits <- c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass")
traits <- c("SLA", "Wood.density","Max.height")
type.filling <- 'species'
files <- c()
for (trait in traits){
for(model in c(model.files.lmer.Tf.1[3], model.files.lmer.Tf.2[3],
model.files.lmer.Tf.3[3], model.files.lmer.Tf.4[3],
model.files.lmer.Tf.5,model.files.lmer.Tf.6)){
for(model in c(model.files.lmer.Tf.1,
model.files.lmer.Tf.2,
model.files.lmer.Tf.3,
model.files.lmer.Tf.5,
model.files.lmer.Tf.CAT.1,
model.files.lmer.Tf.CAT.2)){
source(model, local = TRUE)
model.obj <- load.model()
pathout <- output.dir('lmer', model.obj$name, trait, 'all.log',
type.filling=type.filling)
files <- c(files,file.path(pathout,"results.nolog.all.rds"))
files <- c(files,file.path(pathout,"ecocode.id.results.log.all.rds"))
}
}
out <- lapply(files, summarise.lmer.output.all.list, random.name = 'biomes.id')
names(out) <- lapply(lapply(files,files.details.all),
function(x) paste(as.vector(x[names(x) != 'file']),
collapse="_"))
### remove missing
out <- out[!unlist(lapply(out,FUN=function(x) is.null(x$lmer.summary)))]
saveRDS(out,file='output/list.lmer.out.all.ecocode.id.log.rds')
out <- lapply(files, summarise.lmer.output.all.list)
## LOOP OVER FILES AND NOT SCAN ALL FILES
traits <- c("SLA", "Wood.density","Max.height")
type.filling <- 'species'
files <- c()
for (trait in traits){
for(model in c(model.files.lmer.Tf.1,
model.files.lmer.Tf.2,
model.files.lmer.Tf.3,
model.files.lmer.Tf.5,
model.files.lmer.Tf.CAT.1,
model.files.lmer.Tf.CAT.2)){
source(model, local = TRUE)
model.obj <- load.model()
pathout <- output.dir('lmer', model.obj$name, trait, 'all.log',
type.filling=type.filling)
files <- c(files,file.path(pathout,"species.id.results.log.all.rds"))
}
}
out <- lapply(files, summarise.lmer.output.all.list, random.name = 'biomes.id')
names(out) <- lapply(lapply(files,files.details.all),
function(x) paste(as.vector(x[names(x) != 'file']),
collapse="_"))
### remove missing
out <- out[!unlist(lapply(out,FUN=function(x) is.null(x$lmer.summary)))]
saveRDS(out,file='output/list.lmer.out.all.log.rds')
saveRDS(out,file='output/list.lmer.out.all.species.id.log.rds')
......@@ -3,33 +3,29 @@
source("R/analysis/lmer.output-fun.R")
source("R/analysis/lmer.run.R")
## LOOP OVER FILES AND NOT SCAN ALL FILES
traits <- c("SLA", "Wood.density","Max.height")
type.filling <- 'species'
files <- c()
for (trait in traits){
for (biome in c('Temperate rain forest', 'Temperate grassland desert', 'Tropical rain forest', 'Tropical forest savanna', 'Boreal forest', 'Woodland shrubland', 'Temperate forest')){
for(model in c(model.files.lmer.Tf.1,
model.files.lmer.Tf.2,
model.files.lmer.Tf.3,
model.files.lmer.Tf.4,
model.files.lmer.Tf.5,
model.files.lmer.Tf.CAT.1)){
model.files.lmer.Tf.CAT.1,
model.files.lmer.Tf.CAT.2)){
source(model, local = TRUE)
model.obj <- load.model()
pathout <- output.dir('lmer', model.obj$name, trait, 'all.no.log',
type.filling=type.filling)
name.file <- paste('ecocode.id', gsub(' ', '_',biome),
"results.nolog.all.rds", sep ='.')
files <- c(files,file.path(pathout,name.file))
}
files <- c(files,file.path(pathout,"ecocode.id.results.nolog.all.rds"))
}
}
out <- lapply(files, summarise.lmer.output.all.list)
names(out) <- lapply(lapply(files,files.details.all),
function(x) paste(as.vector(x[names(x) != 'file']),
......@@ -38,54 +34,33 @@ names(out) <- lapply(lapply(files,files.details.all),
out <- out[!unlist(lapply(out,FUN=function(x) is.null(x$lmer.summary)))]
saveRDS(out,file='output/list.lmer.out.all.ecocode.id.no.log.rds')
## ## LOOP OVER FILES AND NOT SCAN ALL FILES
## traits <- c("SLA", "Wood.density","Max.height")
## type.filling <- 'species'
## files <- c()
## for (trait in traits){
## for(model in c(model.files.lmer.Tf.1,
## model.files.lmer.Tf.2,
## model.files.lmer.Tf.3)){
## source(model, local = TRUE)
## model.obj <- load.model()
## pathout <- output.dir('lmer', model.obj$name, trait, 'all.no.log',
## type.filling=type.filling)
## files <- c(files,file.path(pathout,"species.id.results.nolog.all.rds"))
## }
## }
## out <- lapply(files, summarise.lmer.output.all.list, random.name = 'biomes.id')
## names(out) <- lapply(lapply(files,files.details.all),
## function(x) paste(as.vector(x[names(x) != 'file']),
## collapse="_"))
## ### remove missing
## out <- out[!unlist(lapply(out,FUN=function(x) is.null(x$lmer.summary)))]
## saveRDS(out,file='output/list.lmer.out.all.species.id.no.log.rds')
## LOOP OVER FILES AND NOT SCAN ALL FILES
traits <- c("SLA", "Wood.density","Max.height")
type.filling <- 'species'
files <- c()
for (trait in traits){
for(model in c(model.files.lmer.Tf.1,
model.files.lmer.Tf.2,
model.files.lmer.Tf.3,
model.files.lmer.Tf.5,
model.files.lmer.Tf.CAT.1,
model.files.lmer.Tf.CAT.2)){
source(model, local = TRUE)
model.obj <- load.model()
pathout <- output.dir('lmer', model.obj$name, trait, 'all.no.log',
type.filling=type.filling)
files <- c(files,file.path(pathout,"species.id.results.nolog.all.rds"))
## ## LOOP OVER FILES AND NOT SCAN ALL FILES
## traits <- c("SLA", "Wood.density","Max.height")
## type.filling <- 'species'
## files <- c()
## for (trait in traits){
## for(model in c(model.files.lmer.Tf.4b,
## model.files.lmer.Tf.3)){
## source(model, local = TRUE)
## model.obj <- load.model()
## pathout <- output.dir('lmer', model.obj$name, trait, 'all.no.log',
## type.filling=type.filling)
}
}
## files <- c(files,file.path(pathout,"obs.id.results.nolog.all.rds"))
out <- lapply(files, summarise.lmer.output.all.list)
names(out) <- lapply(lapply(files,files.details.all),
function(x) paste(as.vector(x[names(x) != 'file']),
collapse="_"))
### remove missing
out <- out[!unlist(lapply(out,FUN=function(x) is.null(x$lmer.summary)))]
saveRDS(out,file='output/list.lmer.out.all.species.id.no.log.rds')
## }
## }
## out <- lapply(files, summarise.lmer.output.all.list, random.name = 'biomes.id')
## names(out) <- lapply(lapply(files,files.details.all),
## function(x) paste(as.vector(x[names(x) != 'file']),
## collapse="_"))
## ### remove missing
## out <- out[!unlist(lapply(out,FUN=function(x) is.null(x$lmer.summary)))]
## saveRDS(out,file='output/list.lmer.out.all.obs.id.no.log.rds')
TARGETS = $(subst md,pdf,$(shell ls *.md))
all: $(TARGETS)
%.pdf: %.md include.tex
cp ../../../figs/biome.ecocode.xy.pdf biome_ecocode_xy.pdf
cp ../../../figs/world_map_all.png world_map_all.png
pandoc $< --csl=journal-of-applied-ecology.csl --filter pandoc-citeproc --bibliography=refs.bib --template=include.tex --variable mainfont="Times New Roman" --variable sansfont=Arial --variable fontsize=12pt --latex-engine=xelatex -o $@
clean:
rm -f *.pdf
#
This diff is collapsed.
D1 := ../../figs
figs:= biome_ellipse.pdf
all: quick.report.pdf quick.report.html figs
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 $@
quick.report.html: quick.report.md figure/biome_ellipse.pdf
pandoc $< --standalone --webtex -o $@
clean:
rm -f *.pdf
rm -f *.Rmd
rm -f *.html
#
% THIS FILE TAKEN FROM PANDOC DEMOS PAGE:http://johnmacfarlane.net/pandoc/demos.html,http://johnmacfarlane.net/pandoc/demo/mytemplate.tex
\documentclass[$if(fontsize)$$fontsize$,$endif$$if(lang)$$lang$,$endif$, a4paper]{$documentclass$}
\RequirePackage[hmargin=2cm,vmargin=2cm]{geometry}
\usepackage[T1]{fontenc}
\usepackage{lmodern}
\usepackage{amssymb,amsmath}
\usepackage{ifxetex,ifluatex}
\usepackage{fixltx2e} % provides \textsubscript
% use microtype if available
\IfFileExists{microtype.sty}{\usepackage{microtype}}{}
% use upquote if available, for straight quotes in verbatim environments
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
\ifnum 0\ifxetex 1\fi\ifluatex 1\fi=0 % if pdftex
\usepackage[utf8]{inputenc}
$if(euro)$
\usepackage{eurosym}
$endif$
\else % if luatex or xelatex
\usepackage{fontspec}
\ifxetex
\usepackage{xltxtra,xunicode}
\fi
\defaultfontfeatures{Mapping=tex-text,Scale=MatchLowercase}
\newcommand{\euro}{}
$if(mainfont)$
\setmainfont{$mainfont$}
$endif$
$if(sansfont)$
\setsansfont{$sansfont$}
$endif$
$if(monofont)$
\setmonofont{$monofont$}
$endif$
$if(mathfont)$
\setmathfont{$mathfont$}
$endif$
\fi
% Fancy HEADER
\usepackage{fancyhdr}
\pagestyle{fancy}
\pagenumbering{arabic}
\lhead{\itshape{\nouppercase{\leftmark}}}
\chead{}
\rhead{\thepage}
%\lfoot{v $version$}
\cfoot{}
%\rfoot{\thepage}
$if(geometry)$
\usepackage[$for(geometry)$$geometry$$sep$,$endfor$]{geometry}
$endif$
$if(natbib)$
\usepackage{natbib}
\bibliographystyle{plainnat}
$endif$
$if(biblatex)$
\usepackage{biblatex}
$if(biblio-files)$
\bibliography{$biblio-files$}
$endif$
$endif$
$if(listings)$
\usepackage{listings}
$endif$
$if(lhs)$
\lstnewenvironment{code}{\lstset{language=Haskell,basicstyle=\small\ttfamily}}{}
$endif$
$if(highlighting-macros)$
$highlighting-macros$
$endif$
$if(verbatim-in-note)$
\usepackage{fancyvrb}
$endif$
$if(tables)$
\usepackage{longtable}
$endif$
$if(graphics)$
\usepackage{graphicx}
% We will generate all images so they have a width \maxwidth. This means
% that they will get their normal width if they fit onto the page, but
% are scaled down if they would overflow the margins.
\makeatletter
\def\maxwidth{\ifdim\Gin@nat@width>\linewidth\linewidth
\else\Gin@nat@width\fi}
\makeatother
\let\Oldincludegraphics\includegraphics
\renewcommand{\includegraphics}[1]{\Oldincludegraphics[width=\maxwidth]{#1}}
$endif$
\ifxetex
\usepackage[setpagesize=false, % page size defined by xetex
unicode=false, % unicode breaks when used with xetex
xetex]{hyperref}
\else
\usepackage[unicode=true]{hyperref}
\fi
\hypersetup{breaklinks=true,
bookmarks=true,
pdfauthor={$author-meta$},
pdftitle={$title-meta$},
colorlinks=true,
urlcolor=$if(urlcolor)$$urlcolor$$else$blue$endif$,
linkcolor=$if(linkcolor)$$linkcolor$$else$magenta$endif$,
pdfborder={0 0 0}}
$if(links-as-notes)$
% Make links footnotes instead of hotlinks:
\renewcommand{\href}[2]{#2\footnote{\url{#1}}}
$endif$
$if(strikeout)$
\usepackage[normalem]{ulem}
% avoid problems with \sout in headers with hyperref:
\pdfstringdefDisableCommands{\renewcommand{\sout}{}}
$endif$
\setlength{\parindent}{0pt}
\setlength{\parskip}{6pt plus 2pt minus 1pt}
\setlength{\emergencystretch}{3em} % prevent overfull lines
$if(numbersections)$
$else$
\setcounter{secnumdepth}{0}
$endif$
$if(verbatim-in-note)$
\VerbatimFootnotes % allows verbatim text in footnotes
$endif$
$if(lang)$
\ifxetex
\usepackage{polyglossia}
\setmainlanguage{$mainlang$}
\else
\usepackage[$lang$]{babel}
\fi
$endif$
$for(header-includes)$
$header-includes$
$endfor$
$if(title)$
\title{$title$}
$endif$
\author{$for(author)$$author$$sep$ \and $endfor$}
\date{$date$}
\begin{document}
$if(title)$
\maketitle
$endif$
$for(include-before)$
$include-before$
$endfor$
$if(toc)$
{
\hypersetup{linkcolor=black}
\setcounter{tocdepth}{$toc-depth$}
\tableofcontents
}
$endif$
$body$
$if(natbib)$
$if(biblio-files)$
$if(biblio-title)$
$if(book-class)$
\renewcommand\bibname{$biblio-title$}
$else$
\renewcommand\refname{$biblio-title$}
$endif$
$endif$
\bibliography{$biblio-files$}
$endif$
$endif$
$if(biblatex)$
\printbibliography$if(biblio-title)$[title=$biblio-title$]$endif$
$endif$
$for(include-after)$
$include-after$
$endfor$
\end{document}
## % Compare output of growth model with lmer for 3 diffferent subsampling
## This deals with some path issues.
git.root <- function() {
system("git rev-parse --show-toplevel", intern=TRUE)
}
source.root <- function(filename) {
source(file.path(git.root(), filename), chdir=TRUE)
}
readRDS.root <- function(filename) {
readRDS(file.path(git.root(), filename))
}
path.root <- git.root()
source.root("R/analysis/lmer.output-fun.R")
source.root("R/analysis/lmer.run.R")
source.root("R/utils/plot.R")
## load results of a fit with 700000 observations
## TODO LOAD RESULTS FOR 1000000 obs
## TODO LOOK AT RESULTS WITH MAT MAP and by biomes or ecoregions!
## LOOK AT SURVIVAL SEEMS VERY LAREG CI
# NEED TO RUN STAN FOR SURVIVAL
## resampling with no weights
## list.all.results.obs.id <-
## readRDS.root('output/list.lmer.out.all.obs.id.no.log.rds')
## resampling with weights equal to inverse focal species abundance
list.all.results.species.id <-
readRDS.root('output/list.lmer.out.all.species.id.no.log.rds')
## resampling with weights equal to inverse focal ecocode abundance
list.all.results.ecocode.id <-
readRDS.root('output/list.lmer.out.all.ecocode.id.no.log.old2.rds')
list.all.results.ecocode.id.old <-
readRDS.root('output/list.lmer.out.all.ecocode.id.no.log.old2.rds')
list.all.results.ecocode.id.biomes <-
readRDS.root("output/list.lmer.out.all.ecocode.id.biomesno.log.rds")
## check convergence
fun.get.conv(list.all.results.obs.id)
fun.get.conv(list.all.results.species.id)
fun.get.conv(list.all.results.ecocode.id)
## there is convergence problem with all models. In general the optimizer has
## a good converge but lme4 report some error. Which gives only warnings.
## Need to try to re-run with control=lmerControl(optimizer="bobyqa")
## Plot param estimates with biomes random effect for resampling with
## weights per ecocode
names.biomes <- c('Subtrop', 'Temp grassland', 'Medit', 'Temp forest', 'Boreal', 'Temp rain forest',
'Trop rain forest', 'Trop forest', 'Tundra')
pdf(file.path(path.root, 'figs/ecocode_res1.pdf'), width = 10, height =5)
plot.param(list.all.results.ecocode.id ,
model = 'lmer.LOGLIN.ER.AD.Tf.r.biomes.species',
traits = c('Wood.density'),
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]))))) ,
param.print = 1,
col.vec = fun.col.pch.biomes()$col.vec,
pch.vec = fun.col.pch.biomes()$pch.vec,
names.bio = names.biomes )
dev.off()
pdf(file.path(path.root, 'figs/ecocode_res2.pdf'), width = 10, height =5)
plot.param(list.all.results.ecocode.id ,
model = 'lmer.LOGLIN.ER.AD.Tf.r.biomes.species',
traits = c('Wood.density'),
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]))))) ,
param.print = 1:2,
col.vec = fun.col.pch.biomes()$col.vec,
pch.vec = fun.col.pch.biomes()$pch.vec,
names.bio = names.biomes )