Skip to content

Instantly share code, notes, and snippets.

View rmaia's full-sized avatar

Rafael Maia rmaia

View GitHub Profile
@rmaia
rmaia / regplots.R
Created October 28, 2012 23:40
example regression plots with pretty CI shades
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
rmaia / ggplot_chulls.R
Last active November 29, 2023 07:54
Plot convex hulls according to a factor
require(ggplot2)
require(plyr)
# make toy data
set.seed(1)
dataf <- data.frame(
x = rnorm(100),
y = c(rnorm(50,-1), rnorm(50,1)),
cat = rep(c("A","B"), each=50)
)
@rmaia
rmaia / label_specs.R
Created April 8, 2013 12:28
Adding the sample name to the end of the plotted spectra
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
rmaia / radial_phylo_bars.R
Last active December 15, 2015 23:49
plot a radial ("fan") tree with a phenotypic trait as a bar plot
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 / lmer_pbmc.R
Last active December 16, 2015 11:48
Hypothesis testing in mixed-effects linear models based on parametric bootstrap and monte carlo simulations
# 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 / tcs_s3d.R
Last active December 16, 2015 13:09
Use scatterplot3d to plot a tetrahedral color space
# (very convoluted) function to plot tetrahedral colorspace
require(pavo)
require(scatterplot3d)
require(RColorBrewer)
# get the data
data(sicalis)
tcsdat <- tcs(vismodel(sicalis))
@rmaia
rmaia / phylostan_1tree.R
Created July 5, 2013 18:55
Bayesian phylogenetic regression
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 / phylostan_Ktrees.R
Last active May 23, 2017 12:55
Bayesian phylogenetic regression incorporating phylogenetic uncertainty by sampling from multiple trees (work in progress)
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 / mcmcglmm_lambda.R
Last active January 1, 2016 03:28
Does MCMCglmm heritability correspond to lambda estimates?
# load required packages
require(phytools)
require(geiger)
require(MCMCglmm)
data(bird.families)
# CASE 1: very high lambda
@rmaia
rmaia / plot_phylowbars.R
Created May 4, 2016 15:51
Plot a phylogeny with bars, as separate plots, combined
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))