View ubercoldist
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))
View scrapeBirdlife.R
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)
View eigensinplot.R
# 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],
View plot_phylowbars.R
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))
View mcmcglmm_lambda.R
# load required packages
require(phytools)
require(geiger)
require(MCMCglmm)
data(bird.families)
# CASE 1: very high lambda
View phylostan_Ktrees.R
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')
View phylostan_1tree.R
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))
View tcs_s3d.R
# (very convoluted) function to plot tetrahedral colorspace
require(pavo)
require(scatterplot3d)
require(RColorBrewer)
# get the data
data(sicalis)
tcsdat <- tcs(vismodel(sicalis))
View lmer_pbmc.R
# 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)
View radial_phylo_bars.R
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))