Skip to content

Instantly share code, notes, and snippets.

Rafael Maia rmaia

Block or report user

Report or block rmaia

Hide content and notifications from this user.

Learn more about blocking users

Contact Support about this user’s behavior.

Learn more about reporting abuse

Report abuse
View GitHub Profile
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],
@rmaia
rmaia / plot_phylowbars.R
Created May 4, 2016
Plot a phylogeny with bars, as separate plots, combined
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))
@rmaia
rmaia / mcmcglmm_lambda.R
Last active Jan 1, 2016
Does MCMCglmm heritability correspond to lambda estimates?
View mcmcglmm_lambda.R
# load required packages
require(phytools)
require(geiger)
require(MCMCglmm)
data(bird.families)
# CASE 1: very high lambda
@rmaia
rmaia / phylostan_Ktrees.R
Last active May 23, 2017
Bayesian phylogenetic regression incorporating phylogenetic uncertainty by sampling from multiple trees (work in progress)
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')
@rmaia
rmaia / phylostan_1tree.R
Created Jul 5, 2013
Bayesian phylogenetic regression
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))
@rmaia
rmaia / tcs_s3d.R
Last active Dec 16, 2015
Use scatterplot3d to plot a tetrahedral color space
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))
@rmaia
rmaia / lmer_pbmc.R
Last active Dec 16, 2015
Hypothesis testing in mixed-effects linear models based on parametric bootstrap and monte carlo simulations
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)
@rmaia
rmaia / radial_phylo_bars.R
Last active Dec 15, 2015
plot a radial ("fan") tree with a phenotypic trait as a bar plot
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))
You can’t perform that action at this time.