rmaia

ubercoldist <- function(Qi, ni=c(1,2,2,4), w=0.2){

Qi <- Qi[, 1:length(ni)]
Qi <- log(Qi)


e <- setNames(w/sqrt(ni), names(Qi))

# prepare output
pairsid <- t(combn(nrow(Qi),2))

rmaia

require(XML)
scrapeBirdlife <- function(searchterm){

webaddress <- paste("http://www.birdlife.org/datazone/sowbsearchresults.php?a=ns&SearchTerms", searchterm, sep='=')

aa <- readLines(webaddress)
# scrap only web address, common name, species name, IUCN status

aa <- aa[grep('species/factsheet', aa)]
aa <- gsub('<li><a href=\"','',aa)

rmaia

# Plotting eigenvectors from PCA onto scatterplots 

sigmat <- matrix(c(1,-.9,-.9,1), nrow=2)
cent <- c(8,3)

dat <- MASS::mvrnorm(n=100, mu=cent, Sigma=sigmat)

eig <- eigen(cov(dat))

plot(dat, xlim=cent[1]+c(-3,3)*sigmat[1,1], 

rmaia

require(ape)
require(geiger)

# use Geospiza data for example
data(geospiza)

phylo <- drop.tip(geospiza$geospiza.tree, 'olivacea')

bardata <- setNames(geospiza$geospiza.data[,'beakD'], rownames(geospiza$geospiza.data))

rmaia

# load required packages

require(phytools)
require(geiger)
require(MCMCglmm)

data(bird.families) 

# CASE 1: very high lambda

rmaia

require(rstan)
require(geiger)
require(MCMCglmm)

# load data
data(geospiza)
dat <- geospiza$geospiza.data

# create fake sample of trees
tr <- drop.tip(geospiza$geospiza.tree, 'olivacea')

rmaia

require(rstan)
require(geiger)
require(MCMCglmm)

# load data
data(geospiza)

# get the inverse of the phylogenetic variance-covariance matrix
tr <- drop.tip(geospiza$geospiza.tree, 'olivacea')
invA <- solve(vcv.phylo(tr))

rmaia

# (very convoluted) function to plot tetrahedral colorspace

require(pavo)
require(scatterplot3d)
require(RColorBrewer)

# get the data
data(sicalis)
tcsdat <- tcs(vismodel(sicalis))

rmaia

# PARAMETRIC BOOTSTRAP NULL HYPOTHESIS TESTING OF FIXED-EFFECT IN A MIXED MODEL
# Same rationale can be used to test random effects
# adapted from Faraway 2006

require(lme4)

# build the models to be tested
# note that to test fixed effects, one should use ML and not REML estimation

fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy, REML=F)

rmaia

require(ape)
require(geiger)

# use Geospiza data for example
data(geospiza)

phylo <- drop.tip(geospiza$geospiza.tree, 'olivacea')

bardata <- setNames(geospiza$geospiza.data$beakD, rownames(geospiza$geospiza.data))