rmaia

require(RColorBrewer)

# RColorBrewer provides the palettes developped by Cynthia Brewer, you can check them
# out at http://colorbrewer2.org/ .

# make some mock data with interaction

x1 <- factor(rep(c('a','b','c'),each=33))
x2 <- rnorm(99,10,2)

rmaia

data(sicalis)

# plot the desired spectra, leaving room horizontally to add labels

plot(sicalis, select=2:4, xlim=c(300,800))

# find the y position to add the labels right where they finish (wavelength=700)
# using nrow to select the last row (for the last wavelength value)

# subset the last wavelength value (column 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

# 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

# (very convoluted) function to plot tetrahedral colorspace

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

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

rmaia

# load required packages

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

data(bird.families) 

# CASE 1: very high lambda

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

# 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(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

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))