Weecology can has new mammal dataset

Recology has moved - go to http://recology.info/2011/12/weecology-can-has-new-mammal-dataset

So the Weecology folks have published a large dataset on mammal communities in a data paper in Ecology.  I know nothing about mammal communities, but that doesn't mean one can't play with the data...

Their dataset consists of five csv files:  communities, references, sites, species, and trapping data.

Where are these sites, and by the way, do they vary much in altitude?

Let's zoom in on just 'the states'?



What phylogenies can we get for the species in this dataset?
We can use the rOpenSci package treebase to search the online phylogeny repository TreeBASE.  Limiting to returning a max of 1 tree (to save time), we can see that X species are in at least 1 tree on the TreeBASE database.  Nice. 

So there are 321 species in the database with at least 1 tree in the TreeBASE database.  Of course there could be many more, but we limited results from TreeBASE to just 1 tree per query. 

Here's the code:

	# setwd("/Mac/R_stuff/Blog_etc/Mammal_Dataset")
	# URLs for datasets
	comm <- "http://esapubs.org/archive/ecol/E092/201/data/MCDB_communities.csv"
	refs <- "http://esapubs.org/archive/ecol/E092/201/data/MCDB_references.csv"
	sites <- "http://esapubs.org/archive/ecol/E092/201/data/MCDB_sites.csv"
	spp <- "http://esapubs.org/archive/ecol/E092/201/data/MCDB_species.csv"
	trap <- "http://esapubs.org/archive/ecol/E092/201/data/MCDB_trapping.csv"
	# read them
	require(plyr)
	datasets <- llply(list(comm, refs, sites, spp, trap), read.csv, .progress='text')
	str(datasets[[1]]); head(datasets[[1]]) # cool, worked
	llply(datasets, head) # see head of all data.frame's
	# Map the communities
	require(ggplot2)
	require(maps)
	sitesdata <- datasets[[3]]
	sitesdata <- sitesdata[!sitesdata$Latitude == 'NULL',]
	sitesdata$Latitude <- as.numeric(as.character(sitesdata$Latitude))
	sitesdata$Longitude <- as.numeric(as.character(sitesdata$Longitude))
	sitesdata$Elevation_high <- as.numeric(as.character(sitesdata$Elevation_high))
	sitesdata <- sitesdata[sitesdata$Longitude > -140,]
	world_map <- map_data("world")
	us <- world_map[world_map$region=='USA',]
	# World map
	ggplot(world_map, aes(long, lat)) +
	theme_bw(base_size=16) +
	geom_polygon(aes(group = group), colour="grey60", fill = 'white', size = .3) +
	ylim(-55, 85) +
	geom_jitter(data = sitesdata, aes(Longitude, Latitude, size=Elevation_high), alpha=0.3)
	ggsave("worldmap.png")
	# US only
	ggplot(us, aes(long, lat)) +
	theme_bw(base_size=16) +
	geom_polygon(aes(group = group), colour="grey60", fill = 'white', size = .3) +
	ylim(25, 50) +
	xlim(-130, -60) +
	geom_jitter(data = sitesdata, aes(Longitude, Latitude, size=Elevation_high), alpha=0.3)
	ggsave("usmap.png")
	# What phylogenies can we get for these species?
	install.packages("treebase")
	require(treebase); require(RCurl)
	datasets[[4]]$gensp <- as.factor(paste(datasets[[4]]$Genus, datasets[[4]]$Species))
	trees <- llply(datasets[[4]]$gensp, function(x) search_treebase(paste("\"", x, '\"', sep=''),
	by="taxon", max_trees=1, exact_match=TRUE, curl=getCurlHandle()))
	spwithtrees <- data.frame(datasets[[4]]$gensp,
	laply(trees, function(x) if(length(unlist(x)) > 0){1} else{0}))
	names(spwithtrees) <- c('species', 'trees')
	length(spwithtrees[spwithtrees$trees == 1, 2])

view raw mammaldataset.R hosted with ❤ by GitHub

Recology is 1 yr old...

Recology has moved, go to http://recology.info/2011/12/recology-is-1-yr-old

This blog has lasted a whole year already.  Thanks for reading and commenting.

There are a couple of announcements:

1. Less blogging:  I hope to put in many more years blogging here, but in full disclosure, I am blogging for Journal of Ecology now, so I am going to be (and already have been) blogging less here. 
2. More blogging:  If anyone wants to write guest posts at Recology on the topics of using R for ecology and evolution, or open science, please contact me! 
3. Different blogging:  I was going to roll out the new dynamic views for this blog, but Google doesn't allow javascript, which is how I include code using GitHub gists. Oh well... 

Anywho, here is the breakdown of visits to this blog, visualized using #ggplot2, of course.  There were a total of about 23,000 pageviews in the first year of this blog.      Here is the pie chart code I used:

	piechart <- function(set) {
	ggplot(dat[dat$type==set,], aes(x = "", y = pageviews, fill = name)) +
	theme_bw() +
	geom_bar(width = 1) +
	coord_polar(theta = "y") +
	labs(y="", x='') +
	opts(panel.border = theme_rect(colour = 'white'),
	axis.text.x = theme_text(colour='white'))
	}
	# Where dat is our data.frame, and set is the level of a factor in a column called 'type' that you want to plot.

view raw piechart.R hosted with ❤ by GitHub

Visits to top ten posts:

Visits by by pages:

Visits by top referring sites:


Visits by country:


Visits by browsers:


Visits by operating system:

Dynamic views don't support javascript-so reverting back to simple views

Sorry for the temporary loss of GitHub gists...Hopefully dynamic views will support javascript soon!!

I Work For The Internet !

Recology has moved, go to http://recology.info/2011/12/i-work-for-internet

UPDATE: code and figure updated at 647 AM CST on 19 Dec '11.  Also, see Jarrett Byrnes (improved) fork of my gist  here.


The site I WORK FOR THE INTERNET is collecting pictures and first names (last name initials only) to show collective support against SOPA (the Stop Online Piracy Act).  Please stop by their site and add your name/picture.

I used the #rstats package twitteR, created by Jeff Gentry, to search for tweets from people signing this site with their picture, then plotted using ggplot2, and also used Hadley's lubridate to round timestamps on tweets to be able to bin tweets in to time slots for plotting.

Tweets containing the phrase 'I work for the internet' by time:

Here's the code as a GitHub gist.   Sometimes the searchTwitter fxn doesn't returns an error, which I don't understand, but you can play with it:

	require(plyr); require(stringr); require(ggplot2); require(lubridate); require(twitteR)
	datout_1 <- searchTwitter("I work for the internet", n = 1500, since='2011-11-11', until='2011-12-12')
	datout_2 <- searchTwitter("I work for the internet", n = 1500, since='2011-11-13', until='2011-12-14')
	datoutdf <- ldply(c(datout_1, datout_2), function(x) x$toDataFrame(), .progress="text")
	actual <- grep("I work for the internet", datoutdf[,1], ignore.case=T)
	datoutdf2 <- datoutdf[actual,]
	datoutdf2$newtime <- round_date(datoutdf2[,4], "hour")
	datoutdf2_ <- ddply(datoutdf2, .(newtime), summarise,
	numtweets = length(newtime))
	ggplot(droplevels(datoutdf2_), aes(newtime, numtweets)) +
	theme_bw(base_size=18) +
	geom_line(color="red") +
	labs(x='\nTime', y='Number of tweets\n') +
	opts(title="Tweets containing 'I work for the internet'")
	ggsave("iworkfortheinternet.png")

view raw iworkfortheinternet.R hosted with ❤ by GitHub

rOpenSci won 3rd place in the PLoS-Mendeley Binary Battle!

Recology has moved go to http://recology.info/2011/11/ropensci-won-3rd-place-in-plos-mendeley

I am part of the rOpenSci development team (along with Carl Boettiger, Karthik Ram, and Nick Fabina).   Our website: http://ropensci.org/.  Code at Github: https://github.com/ropensci

We entered two of our R packages for integrating with PLoS Journals (rplos) and Mendeley (RMendeley) in the Mendeley-PLoS Binary Battle.  Get them at GitHub (rplos; RMendeley).

These two packages allow users (from R! of course) to search and retrieve data from PLoS journals (including their altmetrics data), and from Mendeley.  You could surely mash up data from both PLoS and Mendeley.  That's what's cool about rOpenSci - we provide the tools, and leave it up to users vast creativity to do awesome things.

3rd place gives us a $1,000 prize, plus a Parrot AR Drone helicopter.

Public vote open for Mendely-PLoS Binary Battle: vote rOpenSci!

Recology has moved, go to http://recology.info/2011/11/public-vote-open-for-mendely-plos
http://www.surveygizmo.com/s3/722753/Mendeley-PLoS-Binary-Battle-Public-Vote

My talk on doing phylogenetics in R

Recology has moved, go to http://recology.info/2011/11/my-talk-on-doing-phylogenetics-in-r

I gave a talk today on doing very basic phylogenetics in R, including getting sequence data, aligning sequence data, plotting trees, doing trait evolution stuff, etc.

Please comment if you have code for doing bayesian phylogenetic inference in R.  I know phyloch has function mrbayes, but can't get it to work...

Phylogenetics in R

View more presentations from schamber

Check out a video of my research at RocketHub

Recology has moved, go to http://recology.info/check-out-video-of-my-research-at

Okay, so this post isn't at all about R - but I can't resist begging my readers for some help.

I’m trying to get some crowdfunding for my research on the evolution of native plants in agricultural landscapes. My campaign is part of a larger project by about 50 other scientists and me to see how well it works to go straight to the public to get funding for science research. All these projects, including mine, are hosted at a site called RocketHub - a site that hosts crowdfunding projects of all sorts – and now they have science.

It is important to get a few bucks at the beginning so that the people that don’t know me with deep pockets will hopefully chip in once they see the money ball rolling.

The funding will go towards paying some students to collect data in the lab for me.

Here’s the link if you want to donate, or just to check out the video I made about my research!
http://www.rockethub.com/projects/3790-evolution-in-agriculture



And watch the video here too:

My little presentation on getting web data through R

Recology has moved, go to http://recology.info/my-little-presentation-on-getting-web
With examples from rOpenSci R packages. p.s. I am no expert at this...
Web data from R

View more presentations from schamber

Two new rOpenSci R packages are on CRAN

Recology has moved, go to http://recology.info/2011/10/two-new-ropensci-r-packages-are-on-cran

Carl Boettiger, a graduate student at UC Davis, just got two packages on CRAN.  One is treebase, which which handshakes with the Treebase API.  The other is rfishbase, which connects with the Fishbase, although I believe just scrapes XML content as there is no API.  See development on GitHub for treebase here, and for rfishbase here.  Carl has some tutorials on treebase and rfishbase at his website here, and we have an official rOpenSci tutorial for treebase here.

Basically, these two R packages let you search and pull down data from Treebase and Fishbase - pretty awesome.  This improves workflow, and puts your data search and acquisition component into your code, instead of being a bunch of mouse clicks in a browser.

These two packages are part of the rOpenSci project.

Two-sex demographic models in R

Recology has moved, go to http://recology.info/2011/10/two-sex-demographic-models-in-r

Tom Miller (a prof here at Rice) and Brian Inouye have a paper out in Ecology (paper, appendices) that confronts two-sex models of dispersal with empirical data.

They conducted the first confrontation of two-sex demographic models with empirical data on lab populations of bean beetles Callosobruchus.

Their R code for the modeling work is available at Ecological Archives (link here).


Here is a figure made from running the five blocks of code in 'Miller_and_Inouye_figures.txt' that reproduces Fig. 4 (A-E) in their Ecology paper (p = proportion female, Nt = density).  Nice!
A: Saturating density dependence
B: Over-compensatory density dependence
C: Sex-specific gamma's (but bM=bF=0.5)
D: Sex-specific b's (but gammaM=gammaF=1)
E: Sex-specific b's (but gammaM=gammaF=2)

New food web dataset

Recology has moved, go to http://recology.info/2011/10/new-food-web-dataset



So, there is a new food web dataset out that was put in Ecological Archives here, and I thought I would play with it. The food web is from Otago Harbour, an intertidal mudflat ecosystem in New Zealand. The web contains 180 nodes, with 1,924 links.

Fun stuff...

igraph, default layout plot

igraph, circle layout plot, nice

My funky little gggraph function plot

get the gggraph function, and make it better, here at Github

	# Blog post about a new food web
	> require(ggplot2); require(igraph); require(bipartite); require(NetIndices)
	# read in edge list
	> mygraph <- read.graph("/Mac/R_stuff/Blog_etc/newfoodweb/edgelist_new.txt",
	+ format = "edgelist")
	> str(mygraph)
	List of 9
	$ : num 181
	$ : logi TRUE
	$ : num [1:1924] 4 4 4 5 5 5 5 5 5 5 ...
	$ : num [1:1924] 1 4 6 1 2 3 4 5 6 12 ...
	$ : num [1:1924] 0 1 2 1812 1813 ...
	$ : num [1:1924] 0 3 61 142 184 186 404 407 413 415 ...
	$ : num [1:182] 0 0 0 0 0 13 61 61 109 126 ...
	$ : num [1:182] 0 0 21 91 114 120 158 220 220 220 ...
	$ :List of 4
	..$ : num [1:2] 1 0
	..$ : Named list()
	..$ : list()
	..$ : list()
	- attr(*, "class")= chr "igraph"
	# igraph graphs
	> plot.igraph(as.undirected(mygraph),
	+ vertex.label = NA,
	+ vertex.size = 1)
	> plot.igraph(as.undirected(mygraph),
	+ layout = layout.circle,
	+ vertex.label = NA,
	+ vertex.size = 1)
	# my own little function, gggraph
	> source("/Mac/R_stuff/Blog_etc/newfoodweb/gggraph.R")
	> gggraph(mygraph, vplace = rnorm, method = "igraph", trophic = "FALSE")

view raw newfoodwb.R hosted with ❤ by GitHub

Phylogenetic community structure: PGLMMs

Recology has moved, go to http://recology.info/2011/10/phylogenetic-community-structure-pglmms


So, I've blogged about this topic before, way back on 5 Jan this year.

Matt Helmus, a postdoc in the Wootton lab at the University of Chicago, published a paper with Anthony Ives in Ecological Monographs this year (abstract here).  The paper addressed a new statistical approach to phylogenetic community structure.

As I said in the original post, part of the power of the PGLMM (phylogenetic generalized linear mixed models) approach is that you don't have to conduct quite so many separate statistical tests as with the previous null model/randomization approach.

Their original code was written in Matlab.  Here I provide the R code that Matt has so graciously shared with me.  There are four functions and a fifth file has an example use case.  The example and output are shown below.

Look for the inclusion of Matt's PGLMM to the picante R package in the future.

Here are links to the files as GitHub gists: 
PGLMM.data.R:  https://gist.github.com/1278205
PGLMM.fit.R:  https://gist.github.com/1284284
PGLMM.reml.R:  https://gist.github.com/1284287
PGLMM.sim.R:  https://gist.github.com/1284288
PGLMM_example.R:  https://gist.github.com/1284442

Enjoy!


The example

	> # Simulate a dataset for the example
	> modelflag <- 1
	> sim.dat <- PGLMM.sim(stree(16, "balanced"), nsites = 30, modelflag = modelflag,
	+ second.env = TRUE, compscale = 1)
	> sim.dat$Vphylo # I think this is the variance-covariance matrix from the phylogeny
	t1 t2 t3 t4 t5 t6 t7 t8 t9 t10 t11 t12 t13 t14 t15 t16
	t1 1.00 0.75 0.50 0.50 0.25 0.25 0.25 0.25 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
	t2 0.75 1.00 0.50 0.50 0.25 0.25 0.25 0.25 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
	t3 0.50 0.50 1.00 0.75 0.25 0.25 0.25 0.25 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
	t4 0.50 0.50 0.75 1.00 0.25 0.25 0.25 0.25 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
	t5 0.25 0.25 0.25 0.25 1.00 0.75 0.50 0.50 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
	t6 0.25 0.25 0.25 0.25 0.75 1.00 0.50 0.50 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
	t7 0.25 0.25 0.25 0.25 0.50 0.50 1.00 0.75 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
	t8 0.25 0.25 0.25 0.25 0.50 0.50 0.75 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
	t9 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00 0.75 0.50 0.50 0.25 0.25 0.25 0.25
	t10 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.75 1.00 0.50 0.50 0.25 0.25 0.25 0.25
	t11 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.50 0.50 1.00 0.75 0.25 0.25 0.25 0.25
	t12 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.50 0.50 0.75 1.00 0.25 0.25 0.25 0.25
	t13 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.25 0.25 0.25 0.25 1.00 0.75 0.50 0.50
	t14 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.25 0.25 0.25 0.25 0.75 1.00 0.50 0.50
	t15 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.25 0.25 0.25 0.25 0.50 0.50 1.00 0.75
	t16 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.25 0.25 0.25 0.25 0.50 0.50 0.75 1.00
	> # PGLMM.sim also automatically outputs three figures, shown below
	> # Organizes the data so that PGLMM can be fit
	> dat <- PGLMM.data(modelflag=modelflag,sim.dat=sim.dat)
	> str(dat) # Structure of the data
	List of 3
	$ YY: num [1:304, 1] 0 0 0 0 1 1 1 0 0 0 ...
	$ VV:List of 2
	..$ Vfullspp : num [1:304, 1:304] 1 0.75 0.5 0.5 0.25 0.25 0.25 0.25 0 0 ...
	..$ Vfullsite: num [1:304, 1:304] 1 1 1 1 1 1 1 1 1 1 ...
	$ XX: num [1:304, 1:16] 1 0 0 0 0 0 0 0 0 0 ...
	> # Simulate a dataset for the example
	> # low number of iterations, maxit = 25 is probably good, but may need
	> # to increase exitcountermax
	> out <- PGLMM.fit(dat = dat, maxit = 25, exitcountermax = 30)
	Loading required package: corpcor
	exitcounter sigma 1 sigma 2 B 1 B 2 B 3 B 4
	0.000000 1.935156 0.315625 -1.321756 -2.140066 -1.029619 -1.029619
	[1] 1.000000000 2.531105003 0.003184433 -1.416240951 -2.293742824 -1.102306778 -1.101311848
	[1] 2.000000000 1.441218975 0.001114315 -1.497821485 -2.407445316 -1.174431967 -1.172422085
	[1] 3.0000000000 1.2060004973 0.0009031993 -1.4099134066 -2.2651017713 -1.1023243928 -1.1012119516
	[1] 4.0000000000 1.0707744528 0.0004762704 -1.4010924746 -2.2545781854 -1.0954740784 -1.0945435442
	[1] 5.0000000000 1.0624400709 0.0002050476 -1.3901393670 -2.2381045469 -1.0864806439 -1.0856860763
	[1] 6.000000e+00 1.141377e+00 3.644753e-05 -1.389817e+00 -2.237711e+00 -1.086224e+00 -1.085433e+00
	[1] 7.0000000000 1.2064541057 0.0000712796 -1.3962534327 -2.2474380739 -1.0915088447 -1.0906400661
	[1] 8.0000000000 1.2064541057 0.0000712796 -1.4017207377 -2.2557858161 -1.0959973521 -1.0950623973
	[1] 9.000000e+00 1.190569e+00 3.077974e-05 -1.401810e+00 -2.255942e+00 -1.096072e+00 -1.095136e+00
	[1] 10.0000000000 1.1694629518 0.0000782201 -1.4004829979 -2.2539137657 -1.0949830703 -1.0940630885
	[1] 11.0000000000 1.1694629518 0.0000782201 -1.3987246019 -2.2512324218 -1.0935389047 -1.0926402314
	> str(out) # Structure of the data
	List of 5
	$ B : num [1:16, 1] -1.399 -2.251 -1.094 -1.093 -0.822 ...
	$ B0 : num [1:16, 1] -1.322 -2.14 -1.03 -1.03 -0.773 ...
	$ s : num [1:2, 1:3] 1.17 7.82e-05 1.05 0.00 1.29 ...
	..- attr(*, "dimnames")=List of 2
	.. ..$ : NULL
	.. ..$ : chr [1:3] "" "0.05" "0.95"
	$ LL : num [1, 1] 454
	$ flag: chr "converged"
	> out$B # coefficients from PGLMM
	[,1]
	[1,] -1.3987338
	[2,] -2.2512486
	[3,] -1.0935466
	[4,] -1.0926478
	[5,] -0.8215762
	[6,] -1.0993153
	[7,] -2.2693828
	[8,] -2.2693828
	[9,] -2.2007681
	[10,] -1.7226929
	[11,] -2.9483086
	[12,] -2.1973124
	[13,] -2.2204310
	[14,] -1.7485642
	[15,] -2.2259044
	[16,] -2.2230874
	> out$B0 # coefficients from standard logistic regression
	[,1]
	[1,] -1.3217558
	[2,] -2.1400662
	[3,] -1.0296194
	[4,] -1.0296194
	[5,] -0.7731899
	[6,] -1.0296194
	[7,] -2.1400662
	[8,] -2.1400662
	[9,] -2.1400662
	[10,] -1.6739764
	[11,] -2.8903718
	[12,] -2.1400662
	[13,] -2.1400662
	[14,] -1.6739764
	[15,] -2.1400662
	[16,] -2.1400662
	> out$s # parameter(s) of the covariance matrix
	0.05 0.95
	[1,] 1.1694629518 1.05248 1.286471
	[2,] 0.0000782201 0.00000 1.153910
	> out$LL # log-likelihood
	[,1]
	[1,] 453.7648
	> out$flag # did model converge?
	[1] "converged"

view raw PGLMM_exampleoutput.R hosted with ❤ by GitHub

..and the figures...

R talk on regular expressions (regex)

Recology has moved, go to http://recology.info/2011/10/r-talk-on-regular-expressions-regex

Regular expressions are a powerful in any language to manipulate, search, etc. data.

For example:

> fruit <- c("apple", "banana", "pear", "pineapple")
> fruit
[1] "apple"     "banana"    "pear"      "pineapple"

> grep("a", fruit) # there is an "a" in each of the words
[1] 1 2 3 4
> 
> strsplit("a string", "s") # strsplit splits the string on the "s"
[[1]]
[1] "a "    "tring"

R base has many functions for regular expressions, see slide 9 of Ed's talk below.  The package stringr, created by Hadley Wickham, is a nice alternative that wraps the base regex functions for easier use. I highly recommend stringr.


Ed Goodwin, the coordinator of the Houston R Users group, gave a presentation to the group last night on regular expressions in R. It was a great talk, and he is allowing me to post his talk here.

Enjoy!  And thanks for sharing Ed!

regex-presentation_ed_goodwin

View more presentations from schamber.

R tutorial on visualizations/graphics

Recology has moved, go to http://recology.info/2011/09/r-tutorial-on-visualizationsgraphics

Rolf Lohaus, a Huxley postdoctoral fellow here in the EEB dept at Rice University, gave our R course a talk on basic visualizations in R this morning.

Enjoy!

	#######################################################
	#### Created by Rolf Lohaus, EEB Dept., Rice University
	####
	#### Part of the code modified and extended from
	#### Robert I. Kabacoff's Quick-R website
	#### (http://www.statmethods.net)
	#######################################################
	############### HIGH-LEVEL PLOTTING FUNCTIONS ###############
	# load some data (dataset of types of cars and their specs used earlier in our meetings)
	data(mtcars)
	# create a basic scatterplot (e.g weight vs. miles per gallon)
	plot(mtcars$wt, mtcars$mpg)
	# specify better axis labels and add a title (note plot always creates a new plot erasing the current one)
	plot(mtcars$wt, mtcars$mpg, xlab="Weight", ylab="Miles per gallon", main="Weight vs. Miles per Gallon")
	##### SAVING GRAPHS #####
	# save the graph either by using 'Export' in the 'Plots' tab in RStudio or using the following functions
	pdf("weight_mpg.pdf") # create PDF in current working directory
	plot(mtcars$wt, mtcars$mpg, xlab="Weight", ylab="Miles per gallon", main="Weight vs. Miles per Gallon")
	dev.off() # "save" the PDF
	# pdf("C:/Graphs/weight_mpg.pdf") # create PDF in some other directory
	# plot(mtcars$wt, mtcars$mpg, xlab="Weight", ylab="Miles per gallon", main="Weight vs. Miles per Gallon")
	# dev.off() # "save" the PDF
	# change the size of the PDF (in inches)
	pdf("weight_mpg.pdf", width=6, height=5) # create PDF in current working directory
	plot(mtcars$wt, mtcars$mpg, xlab="Weight", ylab="Miles per gallon", main="Weight vs. Miles per Gallon")
	dev.off() # "save" the PDF
	# change the font family of the text in the PDF (default is sans-serif "Helvetica")
	pdf("weight_mpg.pdf", width=6, height=5, family="Times") # create PDF in current working directory
	plot(mtcars$wt, mtcars$mpg, xlab="Weight", ylab="Miles per gallon", main="Weight vs. Miles per Gallon")
	dev.off() # "save" the PDF
	png("weight_mpg.png") # create PNG in current working directory
	plot(mtcars$wt, mtcars$mpg, xlab="Weight", ylab="Miles per gallon", main="Weight vs. Miles per Gallon")
	dev.off() # "save" the PNG
	# jpeg("weight_mpg.jpg") # create JPEG in current working directory
	# plot(mtcars$wt, mtcars$mpg, xlab="Weight", ylab="Miles per gallon", main="Weight vs. Miles per Gallon")
	# dev.off() # "save" the JPEG
	##### HISTOGRAMS #####
	# create a basic histogram
	hist(mtcars$mpg)
	hist(mtcars$mpg, xlab="Miles per gallon", main="Histogram of Miles per Gallon")
	# control number of bins with parameter 'breaks'
	hist(mtcars$mpg, breaks=12)
	hist(mtcars$mpg, breaks=seq(10,34))
	# get information about histogram
	hInfo <- hist(mtcars$mpg, breaks=seq(10,34))
	hInfo
	# plot probability density instead of frequencies
	hist(mtcars$mpg, breaks=seq(10,34), freq=F)
	# create a kernel density plot
	plot(density(mtcars$mpg))
	##### BAR PLOTS #####
	# create a basic barplot
	gearCounts <- table(mtcars$gear)
	gearCounts
	barplot(gearCounts)
	barplot(gearCounts, xlab="Number of gears")
	# label the bars individually with parameter 'names.arg'
	barplot(gearCounts, names.arg=c("3 Gears", "4 Gears", "5 Gears"))
	# horizontal bar plot
	barplot(gearCounts, names.arg=c("3 Gears", "4 Gears", "5 Gears"), horiz=TRUE)
	# stacked bar plot
	gearTransmissionCounts <- table(mtcars$am, mtcars$gear)
	gearTransmissionCounts
	barplot(gearTransmissionCounts, names.arg=c("3 Gears", "4 Gears", "5 Gears"))
	# add a legend with automatic labels with parameter 'legend.text'
	barplot(gearTransmissionCounts, names.arg=c("3 Gears", "4 Gears", "5 Gears"), legend.text=T)
	# add a legend with custom labels
	barplot(gearTransmissionCounts, names.arg=c("3 Gears", "4 Gears", "5 Gears"), legend.text=c("automatic", "manual"))
	# grouped bar plot with parameter 'beside=TRUE'
	barplot(gearTransmissionCounts, names.arg=c("3 Gears", "4 Gears", "5 Gears"), legend.text=c("automatic", "manual"), beside=TRUE)
	##### BOX PLOTS #####
	# create a boxplot for an individual variable (e.g. miles per gallon)
	boxplot(mtcars$mpg, ylab="Miles per gallon")
	# create a boxplot for a variable by group (e.g. miles per gallon by cylinders)
	boxplot(mpg ~ cyl, data=mtcars, xlab="Number of cylinders", ylab="Miles per gallon")
	# do not draw outliers
	boxplot(mpg ~ cyl, data=mtcars, xlab="Number of cylinders", ylab="Miles per gallon", outline=FALSE)
	##### LINE PLOTS #####
	# load some time data (dataset of growth of five orange trees over time)
	data(Orange)
	Orange
	# extract data for 1st tree
	tree1 <- subset(Orange, Orange$Tree==1)
	# plot age vs circumference
	plot(tree1$age, tree1$circumference, xlab="Age (days)", ylab="Circumference (mm)")
	# use lines with parameter 'type="l"'
	plot(tree1$age, tree1$circumference, xlab="Age (days)", ylab="Circumference (mm)", type="l")
	# use both points and lines with parameter 'type="b"'
	plot(tree1$age, tree1$circumference, xlab="Age (days)", ylab="Circumference (mm)", type="b")
	# use both points and lines overlayed with parameter 'type="o"'
	plot(tree1$age, tree1$circumference, xlab="Age (days)", ylab="Circumference (mm)", type="o")
	##### PLOTS OF THREE VARIABLES / 3D PLOTS #####
	# image(x, y, z, ...)
	# contour(x, y, z, ...)
	# persp(x, y, z, ...)
	############### LOW-LEVEL PLOTTING FUNCTIONS & GRAPHICAL PARAMETERS ###############
	##### plot lines and points separately #####
	# set up a plot without any points and lines
	plot(tree1$age, tree1$circumference, xlab="Age (days)", ylab="Circumference (mm)", type="n")
	# draw lines
	lines(tree1$age, tree1$circumference)
	# draw points
	points(tree1$age, tree1$circumference)
	# draw points with black border and white fill color
	plot(tree1$age, tree1$circumference, xlab="Age (days)", ylab="Circumference (mm)", type="n")
	lines(tree1$age, tree1$circumference)
	points(tree1$age, tree1$circumference, pch=21, bg="white") # specify point type and background fill color
	##### plot mutiple data in one graph (data for all trees) and add a legend #####
	# get the range for x and y axes
	xRange <- range(Orange$age)
	xRange
	yRange <- range(Orange$circumference)
	yRange
	# set up plot using data
	plot(xRange, yRange, xlab="Age (days)", ylab="Circumference (mm)", type="n")
	# set up plot using x and y axis limits
	plot(0, xlab="Age (days)", ylab="Circumference (mm)", type="n", xlim=xRange, ylim=yRange)
	# draw points and lines for all trees
	for (t in 1:5) {
	tree <- subset(Orange, Orange$Tree==t)
	lines(tree$age, tree$circumference)
	points(tree$age, tree$circumference, pch=21, bg="white")
	}
	# plot data for all trees using different line types
	plot(0, xlab="Age (days)", ylab="Circumference (mm)", type="n", xlim=xRange, ylim=yRange)
	# create a vector with 5 elements specifying the line type to be used for each tree
	lineTypes <- c(2, 1, 4, 5, 3)
	for (t in 1:5) {
	tree <- subset(Orange, Orange$Tree==t)
	lines(tree$age, tree$circumference, lty=lineTypes[t]) # specify line type
	points(tree$age, tree$circumference, pch=21, bg="white")
	}
	# add a legend
	legend("topleft", legend=1:5, lty=lineTypes, title="Tree")
	# plot data for all trees using different line types, point types and colors
	plot(0, xlab="Age (days)", ylab="Circumference (mm)", type="n", xlim=xRange, ylim=yRange)
	# create a vector with 5 elements specifying the point type to be used for each tree
	pointTypes <- seq(21,25)
	# create a vector with 5 elements specifying the color to be used for each tree
	colors <- rainbow(5)
	for (t in 1:5) {
	tree <- subset(Orange, Orange$Tree==t)
	lines(tree$age, tree$circumference, lty=lineTypes[t], col=colors[t])
	points(tree$age, tree$circumference, pch=pointTypes[t], bg="white", col=colors[t])
	}
	# add the corresponding legend but this time omit the legend box and title,
	# and place it in the bottom right corner
	legend("bottomright", legend=c("Tree 1", "Tree 2", "Tree 3", "Tree 4", "Tree 5"), lty=lineTypes, pch=pointTypes, col=colors, pt.bg="white", bty="n")
	##### COLORS #####
	# list all color names
	colors()
	# create color(s) by specifying red, green, blue intensities
	yellowColor = rgb(red=1, green=1, blue=0)
	yellowColor = rgb(red=255, green=255, blue=0, maxColorValue=255)
	transparentGrayColor = rgb(0.3, 0.3, 0.3, alpha=0.5)
	##### calculate and plot mean circumference of trees with error bars ######
	# calculate means at each age
	meanCircumferences <- tapply(Orange$circumference, Orange$age, mean)
	meanCircumferences
	# calculate standard deviation at each age
	sdCircumferences <- tapply(Orange$circumference, Orange$age, sd)
	sdCircumferences
	# get vector of ages
	ages <- as.numeric(names(sdCircumferences))
	ages
	# set up plot
	plot(0, xlab="Age (days)", ylab="Mean circumference (mm)", type="n", xlim=xRange, ylim=yRange)
	# draw lines for means
	lines(ages, meanCircumferences)
	# draw error bars as simple lines
	segments(ages, meanCircumferences-sdCircumferences, ages, meanCircumferences+sdCircumferences)
	# draw points for means
	points(ages, meanCircumferences, pch=21, bg="white")
	# draw error bars with cross bars instead
	plot(0, xlab="Age (days)", ylab="Mean circumference (mm)", type="n", xlim=xRange, ylim=yRange)
	lines(ages, meanCircumferences)
	# draw error bars using 'arrows' function
	arrows(ages, meanCircumferences-sdCircumferences, ages, meanCircumferences+sdCircumferences, angle=90, code=3, length=0.05)
	points(ages, meanCircumferences, pch=21, bg="white")
	##### AXES #####
	# set up plot without box around with parameter 'bty="n'
	plot(0, xlab="Age (days)", ylab="Mean circumference (mm)", type="n", xlim=xRange, ylim=yRange, bty="n")
	# use math and/or symbols in axis labels
	# (for more details see 'help(plotmath)', 'example(plotmath)' and 'demo(plotmath)')
	yLabel = expression(paste("Mean circumference ", bar(italic(C)), " (mm)"))
	plot(0, xlab="Age (days)", ylab=yLabel, type="n", xlim=xRange, ylim=yRange, bty="n")
	# set up plot without axes with parameter 'axes=FALSE'
	# alternatively:
	# - 'xaxt="n"' for no x axis
	# - 'yaxt="n"' for no y axis
	plot(0, xlab="Age (days)", ylab=yLabel, type="n", xlim=c(100, 1600), ylim=c(0, 250), axes=FALSE, bty="n")
	# add custom x axis ('side=1' for bottom side of the graph)
	axis(side=1, at=c(100, 600, 1100, 1600), labels=c("100", "600", "1,100", "1,600"))
	# add custom x axis ('side=2' for left side of the graph)
	axis(side=2, at=c(0, 50, 100, 150, 200, 250))
	lines(ages, meanCircumferences)
	segments(ages, meanCircumferences-sdCircumferences, ages, meanCircumferences+sdCircumferences)
	points(ages, meanCircumferences, pch=21, bg="white")
	##### MULTI-PANEL PLOTS #####
	# split the graph area horizontally in 2 "screens"" (1 row and 2 columns)
	split.screen(c(1,2))
	# select the 1st screen for plotting
	screen(1)
	plot(0, xlab="Age (days)", ylab=expression(paste("Circumference ", italic(C), " (mm)")), type="n", xlim=c(100, 1600), ylim=c(0, 250), axes=FALSE, bty="n")
	axis(side=1, at=c(100, 600, 1100, 1600), labels=c("100", "600", "1,100", "1,600"))
	axis(side=2, at=c(0, 50, 100, 150, 200, 250))
	for (t in 1:5) {
	tree <- subset(Orange, Orange$Tree==t)
	lines(tree$age, tree$circumference, lty=lineTypes[t])
	points(tree$age, tree$circumference, pch=pointTypes[t], bg="white")
	}
	legend("bottomright", legend=c("Tree 1", "Tree 2", "Tree 3", "Tree 4", "Tree 5"), lty=lineTypes, pch=pointTypes, pt.bg="white", bty="n")
	# select the 2nd screen for plotting
	screen(2)
	plot(0, xlab="Age (days)", ylab=yLabel, type="n", xlim=c(100, 1600), ylim=c(0, 250), axes=FALSE, bty="n")
	axis(side=1, at=c(100, 600, 1100, 1600), labels=c("100", "600", "1,100", "1,600"))
	axis(side=2, at=c(0, 50, 100, 150, 200, 250))
	lines(ages, meanCircumferences)
	segments(ages, meanCircumferences-sdCircumferences, ages, meanCircumferences+sdCircumferences)
	points(ages, meanCircumferences, pch=21, bg="white")
	# exit the split-screen mode (close all screens)
	close.screen(all = TRUE)
	# add panel labels
	split.screen(c(1,2))
	screen(1)
	plot(0, xlab="Age (days)", ylab=expression(paste("Circumference ", italic(C), " (mm)")), type="n", xlim=c(100, 1600), ylim=c(0, 250), axes=FALSE, bty="n")
	axis(side=1, at=c(100, 600, 1100, 1600), labels=c("100", "600", "1,100", "1,600"))
	axis(side=2, at=c(0, 50, 100, 150, 200, 250))
	for (t in 1:5) {
	tree <- subset(Orange, Orange$Tree==t)
	lines(tree$age, tree$circumference, lty=lineTypes[t])
	points(tree$age, tree$circumference, pch=pointTypes[t], bg="white")
	}
	legend("bottomright", legend=c("Tree 1", "Tree 2", "Tree 3", "Tree 4", "Tree 5"), lty=lineTypes, pch=pointTypes, pt.bg="white", bty="n")
	# add bold "A" on top left side of plot
	mtext(expression(bold(A)), side=3, adj=0, cex=2, line=1)
	screen(2)
	plot(0, xlab="Age (days)", ylab=yLabel, type="n", xlim=c(100, 1600), ylim=c(0, 250), axes=FALSE, bty="n")
	axis(side=1, at=c(100, 600, 1100, 1600), labels=c("100", "600", "1,100", "1,600"))
	axis(side=2, at=c(0, 50, 100, 150, 200, 250))
	lines(ages, meanCircumferences)
	segments(ages, meanCircumferences-sdCircumferences, ages, meanCircumferences+sdCircumferences)
	points(ages, meanCircumferences, pch=21, bg="white")
	# add bold "B" on top left side of plot
	mtext(expression(bold(B)), side=3, adj=0, cex=2, line=1)
	close.screen(all = TRUE)
	##### GENERAL GRAPHICAL PARAMETERS: MARGINS, TEXT & SYMBOL SIZE, LINE WIDTH, ... #####
	# (for more details see 'help(par)')
	pdf("orangetree_growth.pdf", width=11, height=5)
	# specify margin sizes ('mar', in lines), increase text size of axis labels ('cex.lab')
	# and text size of axes annotation/tick labels ('cex.axis')
	par(mar=c(4.5, 4.5, 2.5, 1.5), cex.lab=1.4, cex.axis=1.2)
	split.screen(c(1,2))
	screen(1)
	plot(0, xlab="Age (days)", ylab=expression(paste("Circumference ", italic(C), " (mm)")), type="n", xlim=c(100, 1600), ylim=c(0, 250), axes=FALSE, bty="n")
	# increase axes line width with 'lwd'
	axis(side=1, at=c(100, 600, 1100, 1600), labels=c("100", "600", "1,100", "1,600"), lwd=1.2)
	axis(side=2, at=c(0, 50, 100, 150, 200, 250), lwd=1.2)
	for (t in 1:5) {
	tree <- subset(Orange, Orange$Tree==t)
	# increase line width with 'lwd'
	lines(tree$age, tree$circumference, lty=lineTypes[t], lwd=1.2)
	# increase point border line width with 'lwd' and point size with 'cex'
	points(tree$age, tree$circumference, pch=pointTypes[t], bg="white", lwd=1.2, cex=1.25)
	}
	legend("bottomright", legend=c("Tree 1", "Tree 2", "Tree 3", "Tree 4", "Tree 5"), lty=lineTypes, pch=pointTypes, pt.bg="white", bty="n", inset=0.016, lwd=1.2, pt.cex=1.25)
	mtext(expression(bold(A)), side=3, adj=0, cex=2, line=1)
	screen(2)
	plot(0, xlab="Age (days)", ylab=yLabel, type="n", xlim=c(100, 1600), ylim=c(0, 250), axes=FALSE, bty="n")
	# increase axes line width with 'lwd'
	axis(side=1, at=c(100, 600, 1100, 1600), labels=c("100", "600", "1,100", "1,600"), lwd=1.2)
	axis(side=2, at=c(0, 50, 100, 150, 200, 250), lwd=1.2)
	# increase line width with 'lwd'
	lines(ages, meanCircumferences, lwd=1.2)
	# increase line width with 'lwd'
	segments(ages, meanCircumferences-sdCircumferences, ages, meanCircumferences+sdCircumferences, lwd=1.2)
	# increase point border line width with 'lwd' and point size with 'cex'
	points(ages, meanCircumferences, pch=21, bg="white", lwd=1.2, cex=1.25)
	mtext(expression(bold(B)), side=3, adj=0, cex=2, line=1)
	close.screen(all = TRUE)
	dev.off()
	############### ADDITIONAL RESOURCES ###############
	# Quick-R: http://www.statmethods.net/index.html
	# R Graph Gallery: http://addictedtor.free.fr/graphiques/
	# R Graphical Manual: http://www.oga-lab.net/RGM2
	# ggplot2 package: http://had.co.nz/ggplot2/
	# lattice package: http://cran.r-project.org/web/packages/lattice/index.html
	# grid package: http://www.stat.auckland.ac.nz/~paul/grid/grid.html
	# Edward Tufte: http://www.edwardtufte.com/tufte/
	# The Visual Display of Quantitative Information: http://www.edwardtufte.com/tufte/books_vdqi

view raw visualizations_tutorial.R hosted with ❤ by GitHub

Short on funding? Can't get a grant? Crowdfunding! #SciFund

Recology has moved, go to http://recology.info/2011/09/short-on-funding-cant-get-grant

Crowdsourced funding is becoming a sustainable way for various artists, entrepreneurs, etc. to get their idea funded from individuals. For example, think of Kickstarter and RocketHub.


Jai Ranganathan and Jarrett Byrnes have started an experiment to determine how well crowdfunding can work for scientists: The SciFund Challenge. Go here to signup and here for their website.


The deadline to sign up is Oct. 1

"Ten Simple Rules for OA Publishers" talk by Philip Bourne

Recology has moved, go to http://recology.info/2011/09/ten-simple-rules-for-oa-publishers-talk

Ten Simple Rules for Open Access Publishers

View more presentations from Philip Bourne

@drewconway interview on @DataNoBorders at the Strata conference

Recology has moved, go to http://recology.info/2011/09/drewconway-interview-on-datanoborders

The O'Reilly Media Strata Summit has many interviews on YouTube (just search YouTube for it)

Drew Conway is the author of a R packages, including infochimps, an R wrapper to the Infochimps API service.

The YouTube video:

Open science talk by Carl Boettiger

Recology has moved, go to http://recology.info/2011/09/open-science-talk-by-carl-boettiger

Carl Boettiger gave a talk on the topic of open science to incoming UC Davis graduate students.

Here is the audio: http://www.archive.org/details/ThingsIWishIKnewThreeYearsAgo-ByTheDavisOpenScienceGroup&reCache=1

Here are the slides: http://hazelnusse.github.com/DOS_WOW2011/#title-slide

My take on an R introduction talk

Recology has moved, go to http://recology.info/2011/09/my-take-on-r-introduction-talk

UPDATE: I put in an R tutorial as a Github gist below.


Here is a short intro R talk I gave today...for what it's worth...

R Introduction

View more presentations from schamber

	# Intro Tutorial to R
	# Scott Chamberlain <myrmecocystus@gmail.com>
	# w/ code edited from Nathan Kraft and Dan Bunker
	########## HELP ##########
	# Getting Help with R on listserves online, etc.
	# Go here for correct way to make reproducible examples to submit to help lists:
	# https://github.com/hadley/devtools/wiki/Reproducibility
	# Code writing style:
	# Google style guide: http://google-styleguide.googlecode.com/svn/trunk/google-r-style.html#generallayout
	# Henrik Bengtsson style guide: http://www1.maths.lth.se/help/R/RCC/
	# Hadley Wickham's style guide: https://github.com/hadley/devtools/wiki/Style
	# Generally use <- instead of = when assigning an object to a name
	########## Installing R ##########
	# go to www.r-project.org, and download a 'precompiled binary'
	# for your system from the 'cran' site nearest you
	# follow links to 'download'
	# follow the directions to install
	# starting R: double click the R icon, or double click a .R file (a script)
	########## R is case sensitive ##########
	myvector <- c(8, 9, 5)
	myvector # OK
	Myvector # wrong case, no result
	########## R as a calculator and simple operations ##########
	# try these simple commands by typing them into the 'console'
	2 + 2
	c(1, 2, 3)
	x = c(1, 2, 3)
	x
	x = c(1:10, 13)
	y = x^2
	y
	ls() # this lists all objects you have created that are still in memory
	rm(x) # Use the rm command to remove objects after you're finished.
	rm(list = ls()) # removes everything!!!
	ls()
	########## Assignments ##########
	# everything in R is an 'object'. That means that R knows something
	# about each object you create, and that there are 'methods'
	# available to work with objects of various classes.
	########## Scalars ##########
	a = 10; b <- 10
	# print to console
	b
	########## Vectors ##########
	x = c(3, 5, 2, 8, 5, 6, 8, 10)
	y = c(1, 3, 4, 5, 3, 9, 8, 3)
	x; y
	########### Lists ##########
	# lists are arbitrary collections of variables
	mylist = c(dd = "foo", ee = 66)
	mylist
	########## Character ##########
	char1 = "this is some text"
	char1
	########## Factors ##########
	# factors are like character vectors but only contain values in predefined 'levels'
	fac1 = factor(rep(c("a", "b"), 4))
	fac1
	########## Matrices ##########
	xx = matrix(10, 2, 3)
	xx
	########## Dataframes ##########
	# dataframes are special objects in R that we will use a lot.
	# They are special because each column of a dataframe can
	# itself be of a different class (UNLIKE MATRICES).
	df = data.frame(x, y, fac1)
	df # print data frame
	summary(df) # basic stats on each column
	str(df) # description of columns and format style
	head(df) # top 6 rows
	tail(df) # bottom 6 rows
	df[ , 3] # pick columns
	df[1 , ] # pick rows
	########## Immutable data frames ##########
	# Special case of a dataframe used in the plyr package created by Hadley
	# Wickham. When summarising data in plyr, can use immutable data frames to
	# do calculations faster, e.g., see the example below
	library(plyr)
	bigdf <- data.frame(a = rnorm(10000), b = rnorm(10000))
	bigdf_i <- idata.frame(bigdf) # make an immutable data frame
	system.time( replicate(100, (ddply(bigdf, .(), mean))) )
	system.time( replicate(100, (ddply(bigdf_i, .(), mean))) ) # faster!
	########## Transforming variables ##########
	# Many options: log10, log, log2, sqrt, asin, etc.
	df$lny=log10(df$y)
	df
	########## Indexing ##########
	# indexing dataframes
	df[1,2] # returns row 1, column 2
	df[1, 'y'] # same
	df[,2] # returns column y as a vector
	df[,'y'] # same
	df$y # same
	df[1:4,] # rows 1-4
	df[df$fac1=="a",] # rows where fac1 is 'a'
	df[df$x>mean(df$x),] # rows where x is above the mean for x
	df[order(df$x), ] # returns whole dataframe sorted by x
	df[df$y%in%c(3:7),] # rows where 'y' matches a list of values
	# indexing lists
	mylist
	mylist[1] # item 1 in list, with label
	mylist[[1]] # the value of item 1
	########## Sequences/regular patterns ##########
	####Sequence, or seq(), is useful for making vectors in regular patterns:
	seq(from=1, to=10, by=1)
	##can also be written as '1:10' in some contexts:
	print(1:10)
	#you can change the interval between elements:
	seq(from=-10, to=30, by=5)
	##same as:
	seq(-10, 30, 5)
	####repeat, or rep(), also useful for regular patterns
	#repeats something as many times as you want:
	rep(5, times=5)
	##same as:
	rep(5,5)
	##more complex:
	rep(c(1,4,5,2), 3)
	##also works with text (aka "strings"):
	rep(c("red", "blue", "gray"), 3)
	## vector function
	vector(mode="character", length=5)
	     ############ BOOLEANS (aka "TRUE" and "FALSE") and tests ##########
	########## Booleans ##########
	# Booleans are a special class in R, and include "TRUE" and "FALSE".
	# Note that they need to be typed in all caps. In an R editor they should
	# turn yellow (or a different color) when typed in. On the Mac,
	# the built-in editor works fine for this; on a PC, you can use an
	# editing program like TINN-R (http://www.sciviews.org/Tinn-R/).
	# You used to be able to abbreviate with T and F, but this is now
	# discouraged and will cause errors in some cases.
	##Many functions include boolean arguments:
	vector <- as.integer(c(seq(1:10), NA))
	sum(vector) ##in this case "sum" fails because of an NA in the data
	?sum ## see that there is an argument "na.rm" that we can set to TRUE
	sum(vector, na.rm=TRUE) ##this removes the NA and sums the remaining 10 integers
	##equivilently you can use the na.omit() command:
	sum(na.omit(vector))
	##Many tests in R return booleans:
	vector > 5 ##this returns true for each part of the vector greater than 5
	vector[which(vector<3)] ##which() can be used to get the indices of an object where the test is TRUE
	vector >= 5 ##etc
	vector == 7 ## this returns true for the element in the vector equal to 7
	##Note that two equals signs ("==") are used for testing equivelence, while one ("=") is used to assign something to a variable.
	##there are a bunch of commands that start with "is." that are very useful:
	is.na(vector) ##shows us that element 11 is an NA
	is.finite(vector) ## checks if values are finite numbers
	##Booleans can be added: TRUE=1,FALSE=0
	sum(is.na(vector)) ##means that there is one NA in the vector
	sum(vector<=3, na.rm=TRUE) ## shows that there are 3 elements of the vector <= to 3.
	## the "!" (known as a "bang" to people with a background in older
	## programming languages) is used as a "not" to reverse TRUE and FALSE:
	is.na(vector)
	!is.na(vector)
	vector != 10 ##which part of the vector is not equal to 10??
	########## "for" loops: ##########
	# for loops are a very powerful way of repeating commands in a simulation
	# or working through each part of an object.
	# the part after "for" sets up a counter variable 'i' (can be called anything),
	# which starts (in this case) at 1 and increments by 1 to 10. Each time through
	# the loop it does whatever is in the { }
	for(i in 1:10){
	print("hello")
	}
	##you can refer to the counter variable in the loop:
	for(j in 1:10){
	print(j)
	}
	##you can use loops to work through each element of a list or a vector:
	colors <- c("red", "blue", "green", "silver")
	colors
	length(colors)
	#this loops adds some text to each element of the 'colors' vector
	for(k in 1:length(colors)){
	print(colors[k])
	colors[k] <- paste(colors[k], "is a color")
	}
	colors
	# You can also use for loops to run through any arbitrary set of values:
	for (k in colors){
	print(k)
	}
	# compare this 'colors' loop with the previous one- here the colors
	# themselves (not the numbers used to index the colors vector) are what
	# the loop moves over.
	## your loop does not have to be sequential from 1:
	for (k in seq(2,10,2)) {
	print(k)
	}
	################ "If/ else" statments ################
	# "if" statements are similar in structure to for loops, and are useful for
	# doing something if a certain condition is met
	##lets set up a boolean (true/ false) variable:
	friendly <- TRUE
	##if statements just do whatever is in the {} if the statement in ()
	# following "if" is TRUE:
	if(friendly){
	print("hello")
	}
	##input the entire statment- you should get a greeting
	##lets change "friendly" to be false:
	friendly<-FALSE
	if(friendly){
	print("hello")
	}
	##now you should get no greeting at all
	## you can follow up the "if" with an "else" statement that gives an
	## alternative if the first test is false:
	if(friendly){
	print("hello")
	}else{
	print("grrrrr...")
	}
	##switch "friendly" and run it again:
	friendly <-!friendly ##remember the "!" means "not"
	if(friendly){
	print("hello")
	}else{
	print("grrrrr...")
	}
	#The statement in the parentheses can also be a comparison that generates a
	# true or false result. In this case, remember that the double equals is used
	# to make a comparison, single equals to assign a value:
	x <- 1
	if (x==1) print('x = 1')
	## notice that you can omit the {} if you only have one command to perform after
	## the if(test...)
	########## Functions ###############
	# Basic example
	greet=function(){
	print("hello there!")
	}
	# this is the format of function writing: a function name of your choosing-
	# "greet" in this case- (action verbs are recommended to keep them distinct
	# from variables), followed by the "=function" plus a () containing any
	# arguments that the function requires. Here the function doesn't require
	# an argument. The code inside the {} is performed whenever the function
	# is called.
	# You need to copy and paste the entire text of the function into the
	# console before you can use the function. There are other ways to do this
	# we'll cover in a little while. Once this is input, type:
	greet()
	## and you should see a message in the console.
	##Here is an example with an argument:
	square=function(x){
	x*x
	}
	##as with function names, you can call arguments anything you want. When you call functions, you can refer to the arguments explictly:
	square(x=2) ## should tell you "4"
	## this should also work without a hitch:
	square(2)
	##once a function is input, you can always see the code for it by typing it in without the ()
	square
	##you can write functions with multiple arguments:
	divide=function(x,y){
	x/y
	}
	##try it out:
	divide(x=8, y=2)
	##same as:
	divide(8,2)
	##not the same as:
	divide(2,8)
	## But if you explictly name the arguments, the order of the arguments does not matter:
	divide(y=2, x=8)
	##You need to specify something for all arguments in this kind of function. For example, this should fail:
	divide(4)
	##You can however specify defaults for certain arguments when you write a function. Here is an example:
	increase=function(x, factor=2){
	x*factor
	}
	##input the function then run these:
	increase(x=3, factor=2)
	##is the same as
	increase(x=3)
	##because we specified a default for 'factor'
	##same as:
	increase(3)
	##but you can also modify the argument "factor" when you call the function even though it has a default:
	increase(3, factor=10)
	##same as
	increase(3,10)
	##the "return" command is very useful in functions for "outputting" what you want from the function. Here we've modified the "increase" function to assign to output of the operation to a variable called "output"
	increase_V2=function(x, factor=2){
	output <- (x*factor)
	}
	##now run:
	increase_V2(50)
	##whoops! No result!
	##try calling the "output" object that we created:
	output ##(this should fail if you've run the script as is)
	##Just as with "for" loops, anything created in the function (such as the 'output' variable) vanishes when it is done running.
	## here is a third version that has added a line of code telling the function to "return" the output variable- that is- send the value of it out to the world once it is done:
	increase_V3=function(x, factor=2){
	output <- (x*factor)
	return(output)
	}
	##now run:
	increase_V3(50)
	##now we have the value of the "output" object to play with again
	##notice that we still do not have a variable named output:
	output ##should fail again
	## Assign the returned value to a variable for further use:
	result <- increase_V3(50)
	result
	##Here's one final example. It is very common to use booleans (true/ false) as arguments with defaults in a function. Lets say you want a function that returns the range of values in a vector, and sometimes you want the absolute minimum and maximum numbers, but sometime you want to know how far apart the extremes are. You could put both options into one function:
	get_range=function(vector, minMax=FALSE){
	vector_min <- min(vector)
	vector_max <- max(vector)
	if(minMax){
	temp <- data.frame(vector_min, vector_max)
	return(temp)
	}else{
	range <- vector_max-vector_min
	return(range)
	}
	}
	##you actually *don't* need the "else" here- but lets keep it in for clarity.
	##once it's input, try it out:
	##here is a practice vector:
	v <- seq(5,20)
	v
	get_range(vector=v, minMax=FALSE)
	##same as:
	get_range(v)
	##different from:
	get_range(v, minMax=TRUE)
	# notice that if minMax is FALSE the function returns a number, but if it is
	# TRUE the function returns a dataframe. This might cause errors in other
	# things you might be doing with the output of this function, so you might
	# want to avoid things like this when writing functions for your own use.
	# as a final note, you can place functions in .r text file and then
	# import them into your workspace with the source() command, such as:
	# source("my_favorite_functions.r"). This can be faster than copying and
	# pasting the functions you always use into the start of a script.
	############ Error catching ##########
	# try
	try()
	# tryCatch
	tryCatch()
	# traceback
	foo <- function(x) { print(1); bar(2) }
	bar <- function(x) { x + a.variable.which.does.not.exist }
	## Not run:
	foo(2) # gives a strange error
	traceback()
	## End(Not run)
	## 2: bar(2)
	## 1: foo(2)
	bar ## Ah, this is the culprit ...
	############ Data Manipulation ##########
	# The following are five packages that do data manipulation
	library(plyr); library(data.table); library(doBy); library(sqldf)
	library(reshape2)
	# 'Pivot tables'
	dataset <- data.frame(var1 = rep(c("a","b","c","d","e","f"), each = 4),
	var2 = rep(c("level1","level1","level2","level2"), 6),
	var3 = rep(c("h","m"), 12), meas = rep(1:12))
	dataset_d <- dcast(dataset, var1 ~ var2 + var3) # pivot table
	melt(dataset_d, id = 1) # pivot table back, not ideal, but works
	# reshape and summarise the dataframe
	data(iris)
	iris_m <- melt(iris, id="Species", na.rm=TRUE)
	head(iris_m)
	dcast(iris_m, Species ~ variable, length) # oops, just length of data vector
	dcast(iris_m, Species ~ variable, identity) # that's better
	# Summarise by grouping levels of variables
	ddply(dataset, .(var1, var2), summarise,
	mean = mean(meas),
	sd = sd(meas))
	# Good way to go from lists to dataframes
	mylist <- list(var1 = c(1,2,3,4),
	var2 = c(5,6,7,8))
	library(plyr) # load plyr package
	ldply(mylist) # makes a dataframe from the list
	####### Looking at your data - pairwise comparison plots #######
	require(lattice)
	require(ggplot2)
	require(car)
	# base graphics
	pairs(iris[1:4])
	# lattice
	splom(~iris[1:4])
	# ggplot2
	plotmatrix(iris[1:4])
	# other
	source("/Users/ScottMac/Dropbox/R Group/Week1_R-Intro/ggcorplot.R")
	ggcorplot(
	data = iris[1:4],
	var_text_size = 5,
	cor_text_limits = c(5,10))
	# other2
	panel.cor <- function(x, y, digits=2, prefix="", cex.cor)
	{
	usr <- par("usr"); on.exit(par(usr))
	par(usr = c(0, 1, 0, 1))
	r <- abs(cor(x, y))
	txt <- format(c(r, 0.123456789), digits=digits)[1]
	txt <- paste(prefix, txt, sep="")
	if(missing(cex.cor)) cex <- 0.8/strwidth(txt)
	test <- cor.test(x,y)
	# borrowed from printCoefmat
	Signif <- symnum(test$p.value, corr = FALSE, na = FALSE,
	cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
	symbols = c("***", "**", "*", ".", " "))
	text(0.5, 0.5, txt, cex = cex * r)
	text(.8, .8, Signif, cex=cex, col=2)
	}
	pairs(iris[1:4], lower.panel=panel.smooth, upper.panel=panel.cor)
	# car package, function scatterplotMatrix (spm abbreviation)
	spm(iris[1:4])
	# compare run times for different methods
	system.time(pairs(iris[1:4]))
	system.time(splom(~iris[1:4]))
	system.time(plotmatrix(iris[1:4]))
	system.time(ggcorplot(
	data = iris[1:4],
	var_text_size = 5,
	cor_text_limits = c(5,10)))
	system.time(pairs(iris[1:4], lower.panel=panel.smooth, upper.panel=panel.cor))
	system.time(spm(iris[1:4]))
	########## Comparison of speed of data summarisation packages #########
	# Modified from the Recipes, scripts, and genomics blog
	library(sqldf)
	library(doBy)
	library(plyr)
	library(data.table)
	n<-10000
	grp1<-sample(1:750, n, replace=T)
	grp2<-sample(1:750, n, replace=T)
	d<-data.frame(x=rnorm(n), y=rnorm(n), grp1=grp1, grp2=grp2, n,
	replace=T)
	# sqldf
	rsqldf <- system.time(sqldf("select grp1, grp2, sum(x), sum(y) from d
	group by grp1, grp2"))
	#doBy
	rdoby <- system.time(summaryBy(x+y~grp1+grp2, data=d, FUN=c(sum)))
	#aggregate
	raggregate <- system.time(aggregate(d, list(d$grp1, d$grp2),
	function(x)sum(x)))
	#plyr
	rplyr <- system.time(ddply(d, .(grp1, grp2), summarise, avx = sum(x),
	avy=sum(y)))
	#data.table
	DT = data.table(d)
	rdataT <- system.time(DT[,list(sum(x),sum(y)),by=list(grp1,grp2)])
	rsqldf
	rdoby
	raggregate
	rplyr
	rdataT
	install.packages("gplots")
	library(gplots)
	x <- c(rdataT[3],rsqldf[3],rdoby[3],raggregate[3],rplyr[3])
	balloonplot( rep("time.elapsed",5),c("data.table","sqldf","doBy","aggregate","plyr"),
	round(x,1), ylab ="Method", xlab="",sorted=F,dotcolor=rev(heat.colors(5)),
	main="time.elapsed for different methods of grouping")
	##### Comparison of lattice and ggplot2 graphics ########
	library(lattice)
	library(ggplot2)
	head(iris)
	# Base graphics
	# Not sure how to do the specific example here?
	# lattice graphics
	pl <- densityplot(~Sepal.Length, data = iris, groups = Species,
	plot.points = FALSE, ref = TRUE)
	print(pl)
	# ggplot2 graphics
	pg <- ggplot(iris, aes(Sepal.Length)) +
	stat_density(geom = "path", position = "identity",
	aes(colour = factor(Species)))
	print(pg)

view raw Rintrotutorial.R hosted with ❤ by GitHub

A Data Visualization Book

Recology has moved, go to http://recology.info/2011/09/data-visualization-book

Note: thanks to Scott for inviting me to contribute to the Recology blog despite being an ecology outsider; my work is primarily in atomic physics. -Pascal
A part of me has always liked thinking about how to effectively present information, but until the past year, I had not read much to support my (idle) interest in information visualization. That changed in the spring when I read Edward Tufte's The Visual Display of Quantitative Information, a book that stimulated me to think more deeply about presenting information. I originally started with a specific task in mind--a wonderful tool for focusing one's interests--but quickly found that Tufte's book was less a practical guide and more a list of general design principles. Then, a few months ago, I stumbled upon Nathan Yau's blog, FlowingData, and found out he was writing a practical guide to design and visualization. Conveniently enough for me, Yau's book, Visualize This, would be released within a month of my discovery of his blog; what follows are my impressions of Visualize This.
I have liked Visualize This a lot.  Yau writes with much the same informal tone as on his blog, and the layout is visually pleasing (good thing, too, for a book about visualizing information!).  The first few chapters are pretty basic if you have done much data manipulation before, but it is really nice to have something laid out so concisely.  The examples are good, too, in that he is very explicit about every step: there is no intuiting what that missing step should be.  The author even acknowledges in the introduction that the first part of the book is at an introductory level.
Early in the book, Yau discusses where to obtain data. This compilation of sources is potentially a useful reference for someone, like me, who almost always generates his own data in the lab. Unfortunately, Yau does not talk much about preparation of (or best practices for) your own data.  Additionally, from the perspective of a practicing scientist, it would have been nice to hear about how to archive data to make sure it is readable far into the future, but that is probably outside the scope of the book.
Yau seems really big into using open source software for getting and analyzing data (e.g. Python, R, etc…), but he is surprisingly attached to the proprietary Adobe Illustrator for turning figures into presentation quality graphics.  He says that he feels like the default options in most analysis programs do not make for very good quality graphics (and he is right), but he does not really acknowledge that you can generate nice output if you go beyond the default settings.  For me, the primary advantage of generating output programmatically is that it is easy to regenerate when you need to change the data or the formatting on the plot.  Using a graphical user interface, like in Adobe Illustrator, is nice if you are only doing something once (how often does that happen?), but when you have to regenerate the darn figure fifty times to satisfy your advisor, it gets tedious to move things around pixel by pixel.
By the time I reached the middle chapters, I started finding many of the details to be repetitive. Part of this repetition stems from the fact that Yau divides these chapters by the type of visualization. For example, "Visualizing Proportions" and "Visualizing Relationships" are two of the chapter titles. While I think these distinctions are important ones for telling the right story about one's data, creating figures for the different data types often boils down to choosing different functions in R or Python. People with less analysis and presentation experience should find the repetition helpful, but I increasingly skimmed these sections as I went along.  
Working through Yau's examples for steps you do not already know would probably be the most useful way of getting something out of the book.  So, for example, I started trying to use Python to scrape data from a webpage, something I had not previously done.  I followed the book's example of this data-scraping just fine, but as with most things in programming, you find all sorts of minor hurdles to clear when you try your own thing. In my case, I am re-learning the Python I briefly learned about 10 years ago--partly in anticipation of not having access to Matlab licenses once I vacate the academy--since I have forgotten a lot of the syntax.  A lot of this stuff would be faster if I were working in Matlab which I grew more familiar with in graduate school.
Overall, Visualize This is a really nice looking book and will continue to be useful to me as a reference. Yau concludes his book with a refreshing reminder to provide context for the data we present. This advice is particularly relevant when presenting to a wider or lay audience, but it is still important for us, as scientists, to clearly communicate our findings in the literature. Patterns in the data are not often self-evident, and therefore we should think carefully about which visualization tools will best convey the meaning of our results.
Edited to add a link to Visualize This here and in the introductory paragraph.

FigShare Talk

Recology has moved, go to http://recology.info/2011/09/figshare-talk

FigShare - I very much like this idea of a place to put your data online that is NOT published. Dryad is a nice place for datastes linked with published papers, but there isn't really a place for datasets that perhaps did not make the cut for a published paper, and if known to the scientific community, could potentially help resolve the "file-drawer" effect in meta-analyses. (wow, run on sentence)

 
"Figshare - Why don't you publish all your research?" Mark Hahnel Imperial College London from London Biogeeks on Vimeo.

Jonathan Eisen on the Fourth Domain and Open Science

Recology has moved, go to http://recology.info/2011/09/jonathan-eisen-on-fourth-domain-and




Stalking the Fourth Domain with Jonathan Eisen, Ph D from mendelspod on Vimeo.

rnpn: An R interface for the National Phenology Network

Recology has moved, go to http://recology.info/2011/08/rnpn-r-interface-for-national-phenology

The team at rOpenSci and I have been working on a wrapper for the USA National Phenology Network API. The following is a demo of some of the current possibilities. We will have more functions down the road. Get the publicly available code, and contribute, at Github here. If you try this out look at the Description file for the required R packages to run rnpn. Let us know at Github (here) or at our website  http://ropensci.org/, or in the comments below, or on twitter (@rOpenSci), what use cases you would like to see with the rnpn package.

Method and demo of each:
Get observations for species by day
From the documentation: "This function will return a list of species, containing all the dates which observations were made about the species, and a count of the number of such observations made on that date."

#### Note, the data below is truncated for blogging brevity...

> getobsspbyday(c(1, 2), '2008-01-01', '2011-12-31') # Searched for species 1 and 2 from Jan 1, 2008 to Dec 31, 2011
          date count   species
1   2009-03-08     2 species 1
2   2009-03-15     1 species 1
3   2009-03-22     1 species 1
4   2009-03-24     1 species 1
5   2009-03-26     1 species 1
6   2009-04-17     1 species 1
7   2009-04-24     1 species 1
8   2009-05-12     1 species 1
9   2009-05-20     1 species 1
10  2009-11-24     1 species 1
11  2009-12-07     1 species 1
12  2010-01-18     1 species 1
13  2010-01-23     1 species 1
62  2011-05-29     1 species 1
63  2011-06-27     1 species 1
64  2011-06-30     2 species 1
65  2009-03-17     1 species 2
66  2009-04-03     3 species 2
67  2009-04-05     3 species 2
68  2009-04-10     3 species 2
69  2009-04-17     3 species 2

Get individuals at specific stations
From the documentation: "This function returns all of the individuals at a series of stations."

> getindsatstations(c(507, 523)) # Searched for any individuals at stations 507 and 523
   individual_id individual_name species_id kingdom
1           1200         dogwood         12 Plantae
2           1197    purple lilac         36 Plantae
3           1193         white t         38 Plantae
4           3569     forsythia-1         73 Plantae
5           1206            jack        150 Plantae
6           1199      trout lily        161 Plantae
7           1198           dandy        189 Plantae
8           1192           red t        192 Plantae
9           1710    common lilac         36 Plantae
10          1711  common lilac 2         36 Plantae
11          1712       dandelion        189 Plantae

Get individuals of species at stations
From the documentation: "This function will return a list of all the individuals, which are members of a species, among  any number of stations."

> getindspatstations(35, c(60, 259), 2009)  # Search for individuals of species 35 at stations 60 and 259 in year 2009
  individual_id individual_name number_observations
1          1715            west                   5
2          1716            east                   5

Get observation associated with particular observation
From the documentation: "This function will return the comment associated with a particular observation."

> getobscomm(1938) # The observation for observation number 1938

$observation_comment
[1] "some lower branches are bare"

Tenure track position in systematics at the University of Vermont

Recology has moved, go to http://recology.info/2011/08/tenure-track-position-in-systematics-at


There is an awesome position opening up for an assistant professor in systematics at the University of Vermont. Below is the announcement, and see the original post at the Distributed Ecology blog. Why is this related to R? One can do a lot of systematics work in R, including retrieving scientific collections data through an upcoming package handshaking with VertNet (part of the rOpenSci project), managing large data sets, retrieval of GenBank data through the ape package (see fxn read.genbank), phylogenetic reconstruction and analysis, and more. So I am sure a systematist with R ninja skills will surely have a head up on the rest of the field. 


Assistant Professor in Systematics

Department of Biology
University of Vermont
Burlington, Vermont

The Department of Biology of the University of Vermont seeks applications for a tenure- track Assistant Professor position in Systematics and Evolutionary Biology of arthropods, especially insects. The position will be open in the fall of 2012. The successful candidate will have expertise in classical and molecular systematics, including analysis of complex data sets. Candidates pursuing phylogenomics and innovative methods in bioinformatics in combination with taxonomy are especially encouraged to apply. Department information at: http://www.uvm.edu/~biology/.


All applicants are expected to: 1) hold a Ph.D. degree in relevant disciplines and have two or more years of postdoctoral experience; 2) develop a competitively funded research program; 3) teach undergraduate courses (chosen from among general biology, evolution, systematic entomology, and others in the candidate's expertise); 4) teach, mentor and advise undergraduate and graduate students; and 5) oversee a natural history collection of historic significance.

Candidates must apply online: http://www.uvmjobs.com/. On left see "Search Postings" then find "Biology" under "HCM Department" then posting 0040090 (first posting). Sorry, but we cannot supply the direct link because it will time out.

Attach a cover letter with a statement of research focus and teaching interests (one document), a curriculum vitae, representative publications, and the contact information of three references.

Review of applications will begin on September 15, 2011, and will continue until the position is filled. Questions and up to three additional publications may be directed to Dr. Jos. J. Schall: jschall@uvm.edu.

The University of Vermont recently identified several "Spires of Excellence" in which it will strategically focus institutional investments and growth over the next several years. One spire associated with the position is Complex Systems. Candidates whose research interests align with this spire are especially encouraged to applyhttp://www.uvm.edu/~tri/.
The University seeks faculty who can contribute to the diversity and excellence of the academic community through their research, teaching, and/or service. Applicants are requested to include in their cover letter information about how they will further this goal. The University of Vermont is an Affirmative Action/Equal Opportunity employer. The Department is committed to increasing faculty diversity and welcomes applications from women, underrepresented ethnic, racial and cultural groups, and from people with disabilities.

Thursday at #ESA11

Recology has moved, go to http://recology.info/2011/08/thursday-at-esa11

Interesting talks/posters:

Richard Lankau presented research on trade-offs and competitive ability. He suggests that during range expansion selection for increased intraspecific competitive ability in older populations leads to loss of traits for interspecific competitive traits due to trade-offs between these traits.


Ellner emphatically states that rapid evolution DOES matter for ecological responses, and longer-term evolutionary patterns as well. [His paper on the talk he was giving came out prior to his talk, which he pointed out, good form sir]


Lauren Sullivan gave an interesting talk on bottom up and top down effects on plant reproduction in one site of a huge network of sites doing similar nutrient and herbivory manipulations around the globe - NutNet (go here: http://nutnet.science.oregonstate.edu/).


Laura Prugh shows in California that the engineering effects (i.e., the mounds that they make) of giant kangaroo rats are more important for the associated food web than the species interaction effects (the proxy used was just density of rats).


Kristy Deiner suggests that chironomids are more phylogenetic similar in lakes with stocked fish relative to fishless lakes, in high elevation lakes in the Sierra Nevada. She used barcode data to generate her phylogeny of chironomids. If you have barcode data and want to search BOLD Systems site, one option is doing it from R using rbold, a package under development at rOpenSci (code at Github).


Jessica Gurevitch presented a large working group's methods/approach to a set of reviews on invasion biology. We didn't get to see a lot of results from this work, but I personally was glad to see her explaining to a packed room the utility of meta-analysis, and comparing to the medical field in which meta-analysis is sort of the gold standard by which to draw conclusions.


Following Jessica, Jason Fridley told us about the Evolutionary Imbalance Hypothesis (EIH) (see my notes here). He posed the problem of, when two biotas come together, what determines which species are retained in this new community and which species are left out. He listed a litany of traits/responses to measure to get at this problem, but suggested that with a little bit of "desktop ecology", we could simply ask: Is the invasability of X region related to the phylogenetic diversity of that region? In three destination regions (Eastern Deciduous Forests, Mediterranean California, and the Czech Republic) out of four there was a positive relationship between proportion of invasive plant species in a source region and the phylogenetic diversity of the source regions.

Wednesday at #ESA11

Recology has moved, go to http://recology.info/2011/08/wednesday-at-esa11

Interesting talks/posters:

• Ethan White's poster describing EcologicalData.org was of course awesome given my interest in getting data into the hands of ecologists over at rOpenSci.org. Ethan also has software you can download on your machine to get the datasets you want easily - EcoData Retriever. [rOpenSci will try to take advantage of their work and allow you to call the retriever from R]
• Carl Boettiger's talk was awesome. He explained how we need better tools to be able to predict collapses using early warning signals. He developed a way to estimate the statistical distribution of probabilities of system collapse. 
• Jennifer Dunne: Explained how she put together an ancient network from Germany. Bravo. 
• Carlos Melian explained his model of network buildup that starts from individuals, allows speciation, and other evolutionary processes. 
• Rachel Winfree told us that in two sets of mutualistic plant-pollinator networks in New Jersey and California, that the least connected pollinator species were the most likely to be lost from the network with increasing agricultural intensity. 
• Dan Cariveau suggests that pollination crop services can be stabilized even with increasing agriculture intensity if in fact pollinator species respond in different ways. That is, some pollinators may decrease in abundance with increasing ag intensity, while other species may increase - retaining overall pollination services to crops.

Monday at #ESA11

Recology has moved, go to http://recology.info/2011/08/monday-at-esa11

Monday was a good day at ESA in Austin. There were a few topics I promised to report on in my blogging/tweeting.


...focused on open source data. Carly Strasser's presentation on guidelines for data management was awesome (including other talks in the symposium on Creating Effective Data Management Plans for Ecological Research). Although this was a good session, I can't help but wish that they had hammered home the need for open science more. Oh well. Also, they talked a lot about how, and not a lot of why we should properly curate data. Still, a good session. One issue Carly and I talked about was tracking code in versioning systems such as Github. There doesn't seem to be a culture of versioning code for analyses/simulations in ecology. But when we get there...it will be easier to share/track/collaborate on  code.

...used R software. David Jennings talked about a meta-analysis asking if phylogenetic distance influences competition strength in pairwise experiments. David used the metafor package in R to do his meta-analysis. Good form sir.

...did cool science. Matt Helmus presented a great talk on phylogenetic species area curves (likely using R, or Matlab maybe?).


p.s. We launched rOpenSci today.


Oddities:

• The Tilman effect - Tilman's talk was so packed it looked like there was a line waiting to get into a trendy bar. Here's a picture (credit: Jaime Ashander). Bigger room next time anyone? 
• Wiley came out with an open source journal called Ecology and Evolution. This brings them to 3 open source journals (the other two are in other fields). We (rOpenSci) will attempt to hand-shake with these journals. 
• The vegetarian lunch option was surprisingly good. Nice. 

(#ESA11) rOpenSci: a collaborative effort to develop R-based tools for facilitating Open Science

Recology has moved, go to http://recology.info/2011/08/esa11-ropensci-collaborative-effort-to_08

Our development team would like to announce the launch of rOpenSci. As the title states, this project aims to create R packages to make open science more available to researchers.

http://ropensci.org/

What this means is that we seek to connect researchers using R with as much open data as possible, mainly through APIs. There are a number of R packages that already do this (e.g., infochimps, twitteR), but we are making more packages, e.g., for Mendeley, PLoS Journals, and taxonomic sources (ITIS, EOL, TNRS, Phylomatic, UBio).

Importantly, we are creating a package called rOpenSci, which aims to integrate functions from packages for individual open data sources.

If you are somewhat interested, follow our progress on our website, on Twitter, or contact us. If you are really^2 interested you could go to Github and contribute. If  you are really^3 interested, join our development team.

Blogging/tweeting from #ESA11

Recology has moved, go to http://recology.info/2011/07/bloggingtweeting-from-esa11

I will be blogging about the upcoming Ecological Society of America meeting in Austin, TX. I will focus on discussing talks/posters that:

1. Have taken a cool approach to using data, or
2. Have focused on open science/data, or
3. Done something cool with R software, or
4. Are just exciting in general

I will also tweet throughout the meeting from @recology_ (yes the underscore is part of the name, recology was already taken). 

The hashtag for the meeting this year is #ESA11

Models in Evolutionary Ecology seminar, organized by Timothee Poisot

Recology has moved, go to http://recology.info/2011/07/models-in-evolutionary-ecology-seminar

Here is one of the talks by Thomas Brouquet, and see the rest here.



Thomas Broquet by mez_en_video

Archiving ecology/evolution data sets online

We now have many options for archiving data sets online:

Dryad, KNB, Ecological Archives, Ecology Data Papers, Ecological Data, etc.

However, these portals largely do not communicate with one another as far as I know, and there is no way to search over all data set sources, again, as far as I know. So, I wonder if it would ease finding of all these different data sets to get these different sites to get their data sets cloned on a site like Infochimps, or have links from Infochimps.  Infochimps already has APIs (and there's an R wrapper for the Infochimps API already set up here: http://cran.r-project.org/web/packages/infochimps/index.html by Drew Conway), and they have discussions set up there, etc.

Does it make sense to post data sets linked to published works on Infochimps? I think probably not know that I think about it. But perhaps it makes sense for other data sets, or subsets of data sets that are not linked with published works to be posted there as I know at least Dryad only accepts data sets linked with published papers.

One use case is there was a tweet from someone recently that his students were excited about getting their data sets on their resume/CV, but didn't think there was a way to put them any place where there wasn't a precondition that the data set was linked with a published work. Seems like this could be a good opportunity to place these datasets on Infcohimps, and at least they are available then where a lot of people are searching for data sets, etc.

What I think would be ideal is if Dryad, KNB, etc. could link their datasets to Infochimps, where they could be found, then users can either get them from Infochimps, or perhaps you would have to go to the Dryad site, e.g. But at least you could search over all ecological data sets then.

CRdata vs. Cloudnumbers

Cloudnumbers and CRdata are two new cloud computing services.


I tested the two services with a very simple script. The script simply creates a dataframe of 10000 numbers via rnorm, and assigns them to a factor of one of two levels (a or b). I then take the mean of the two factor levels with the aggregate function.


In CRdata you need to put in some extra code to format the output in a browser window. For example, the last line below needs to have '<crdata_object>' on both sides of the output object so it can be rendered in a browser. And etc. for other things that one would print to a console. Whereas you don't need this extra code for using Cloudnumbers.

 
dat <- data.frame(n = rnorm(10000), p = rep(c('a','b'), each=5000))
 
out <- aggregate(n ~ p, data = dat, mean)
 
#<crdata_object>out</crdata_object>

Here is a screenshot of the output from CRdata with the simple script above.

This simple script ran in about 20 seconds or so from starting the job to finishing. However, it seems like the only output option is html. Can this be right? This seems like a terrible only option.


In Cloudnumbers you have to start a workspace, upload your R code file.
Then, start a session...
choose your software platform...
choose packages (one at a time, very slow)...
then choose number of clusters, etc.
Then finally star the job.
Then it initializes, then finally you can open the console, and
Then from here it is like running R as you normally would, except on the web.


Who wins (at least for our very minimal example above)

1. Speed of entire process (not just running code): CRdata
2. Ease of use: CRdata
3. Cost: CRdata (free only)
4. Least annoying: Cloudnumbers (you don't have to add in extra code to run your own code)
5. Opensource: CRdata (you can use publicly available code on the site)
6. Long-term use: Cloudnumbers (more powerful, flexible, etc.)

I imagine Cloudnumbers could be faster for larger jobs, but you would have to pay for the speed of course. 

What I really want to see is a cloud computing service that accepts code directly run from R or RStudio. Hmmm...that would be so tasty indeed. I think Cloudnumbers may be able to do this, but haven't tested it yet.  

Perhaps using the server version of RStudio along with Amazon's EC2 is a better option than both of these. See Karthik Ram's post about using RStudio server along with Amazon's EC2. Even just running RStudio server on your Unbuntu machine or virtual machine is a pretty cool option, even without EC2 (works like a charm on my Parallels Ubuntu vm on my Mac). 

rbold: An R Interface for Bold Systems barcode repository

Have you ever wanted to search and fetch barcode data from Bold Systems?

I am developing functions to interface with Bold from R. I just started, but hopefully folks will find it useful.

The code is at Github here. The two functions are still very buggy, so please bring up issues below, or in the Issues area on Github. For example, some searches work and other similar searches don't. Apologies in advance for the bugs.

Below is a screenshot of an example query using function getsampleids to get barcode identifiers for specimens. You can then use getseqs function to grab barcode data for specific specimens or many specimens.

iEvoBio 2011 Synopsis

We just wrapped up the 2011 iEvoBio meeting. It was awesome! If you didn't go this year or last year, definitely think about going next year.

Here is a list of the cool projects that were discussed at the meeting (apologies if I left some out):

1. Vistrails: workflow tool, awesome project by Claudio Silva
2. Commplish: purpose is to use via API's, not with the web UI
3. Phylopic: a database of life-form silouhettes, including an API for remote access, sweet!
4. Gloome
5. MappingLife: awesome geographic/etc data visualization interace on the web
6. SuiteSMA: visualizating multiple alignments
7. treeBASE: R interface to treebase, by Carl Boettiger
8. VertNet: database for vertebrate natural history collections
9. RevBayes: revamp of MrBayes, with GUI, etc. 
10. Phenoscape Knowledge Base
• Peter Midford lightning talk: talked about matching taxonomic and genetic data
11. BiSciCol: biological science collections tracker
12. Ontogrator 
13. TNRS: taxonomic name resolution service
14. Barcode of Life data systems, and remote access
15. Moorea Biocode Project
16. Microbial LTER's data
17. BirdVis: interactive bird data visualization (Claudio Silva in collaboration with Cornell Lab of Ornithology)
18. Crowdlabs: I think the site is down right now, another project by Claudio Silva
19. Phycas: Bayesian phylogenetics, can you just call this from R?
20. RIP MrBayes!!!! replaced by RevBayes (see 9 above)
21. Slides of presentations will be at Slideshare (not all presentations up yet)          
22. A birds of a feather group I was involved in proposed an idea (TOL-o-matic) like Phylomatic, but of broader scope, for easy access and submission of trees, and perhaps even social (think just pushing a 'SHARE' button within PAUP, RevBayes, or other phylogenetics software)! 
23. Synopses of Birds of a Feather discussion groups: http://piratepad.net/iEvoBio11-BoF-reportouts

PLoS journals API from R: "rplos"

The Public Libraries of Science (PLOS) has an API so that developers can create cool tools to access their data (including full text papers!!).

Carl Boettiger at UC Davis and I are working on R functions that use the PLoS API. See our code on Github here. See the wiki at the Github page for examples of use. We hope to deploy rplos as a package someday soon. Please feel free to suggest changes/additions rplos in the comments below or on the Github/rplos site.

Get your own API key here.

ggplot2 talk by Hadley Whickam at Google

OpenStates from R via API: watch your elected representatives

I am writing some functions to acquire data from the OpenStates project, via their API. They have a great support community at Google Groups as well.

On its face this post is not obviously about ecology or evolution, but well, our elected representatives do, so to speak, hold our environment in a noose, ready to let the Earth hang any day.

Code I am developing is over at Github.

Here is an example of its use in R, in this case using the Bill Search option (billsearch.R on my Github site), and in this case you do not provide your API key in the function call, but instead put it in your .Rprofile file, which is called when you open R. We are searching here for the term 'agriculture' in Texas ('tx'), in the 'upper' chamber.

> temp <- billsearch('agriculture', state = 'tx', chamber = 'upper')
 
> length(temp)
[1] 21
 
> temp[[1]]
$title
[1] "Congratulating John C. Padalino of El Paso for being appointed to the United States Department of Agriculture."
 
$created_at
[1] "2010-08-11 07:59:46"
 
$updated_at
[1] "2010-09-02 03:34:39"
 
$chamber
[1] "upper"
 
$state
[1] "tx"
 
$session
[1] "81"
 
$type
$type[[1]]
[1] "resolution"
 
 
$subjects
$subjects[[1]]
[1] "Resolutions"
 
$subjects[[2]]
[1] "Other"
 
 
$bill_id
[1] "SR 1042"

Created by Pretty R at inside-R.org


Apparently, the first bill (SR 2042, see $bill_id at the bottom of the list output) that came up was to congratulate John Paladino for being appointed to the USDA.

The other function I have ready is getting basic metadata on a state, called statemetasearch.

I plan to develop more functions for all the possible API calls to the OpenStates project.

Just turned on 'mobile template' for this blog, let me know if it works

> 10,000 visits to Recology!!!!

Thanks so much everyone for stopping by!

How to fit power laws

A new paper out in Ecology by Xiao and colleagues (in press, here) compares the use of log-transformation to non-linear regression for analyzing power-laws.

They suggest that the error distribution should determine which method performs better. When your errors are additive, homoscedastic, and normally distributed, they propose using non-linear regression. When errors are multiplicative, heteroscedastic, and lognormally distributed, they suggest using linear regression on log-transformed data. The assumptions about these two methods are different, so cannot be correct for a single dataset.

They will provide their R code for their methods once they are up on Ecological Archives (they weren't up there by the time of this post).

searching ITIS and fetching Phylomatic trees

I am writing a set of functions to search ITIS for taxonomic information (more databases to come) and functions to fetch plant phylogenetic trees from Phylomatic. Code at github.

Also, see the examples in the demos folder on the Github site above.

	# A species list, which ecologies often have, but may not know the family names, or if they do,
	# they (meaning me) may spell the names wrong from time to time
	> splist <- c("annona_cherimola", 'annona_muricata', "quercus_robur",
	"shorea_robusta", "pandanus_patina", "oryza_sativa", "durio_zibethinus",
	"rubus_ulmifolius", "asclepias_curassavica", "pistacia_lentiscus")
	# Get the Taxonomic Serial Number for each species, which are unique (note that there could be many TSN's within a species
	# for each subspecies for example. So you may want to query for you species using get_itis_xml in the demos on the github
	# site and then make sure you get the TSN for the right subspecies or species
	> tsns <- laply(splist, get_tsn, searchtype = "sciname", by_ = "name",
	.progress = "text")
	|===========================================================================================| 100%
	# Get Phylomatic formatted strings for each species (e.g., family/genus/genus_epithet)
	> dat_ <- laply(tsns, get_phymat_format, format='rsubmit', .progress="text")
	|===========================================================================================| 100%
	# Send the strings to Phylomatic website via an API, and retrieve the newick tree (alternatively, get xml output)
	> tree <- get_phylomatic_tree(dat_, 'TRUE', 'GET', 'new', 'TRUE')
	Loading required package: ape
	# Plot the tree, or do anything else you want to do!
	> plot(tree)

view raw taxize_example.R hosted with ❤ by GitHub

phylogenetic signal simulations

I did a little simulation to examine how K and lambda vary in response to tree size (and how they compare to each other on the same simulated trees). I use Liam Revell's functions fastBM to generate traits, and phylosig to measure phylogenetic signal.

Two observations: 


First, it seems that lambda is more sensitive than K to tree size, but then lambda levels out at about 40 species, whereas K continues to vary around a mean of 1.

Second, K is more variable than lambda at all levels of tree size (compare standard error bars).

Does this make sense to those smart folks out there?

	#### Simulations
	install.packages(c("ape","reshape2","ggplot2"))
	require(ape); require(reshape2); require(ggplot2)
	source("http://anolis.oeb.harvard.edu/~liam/R-phylogenetics/phylosig/v0.3/phylosig.R")
	source("http://anolis.oeb.harvard.edu/~liam/R-phylogenetics/fastBM/v0.4/fastBM.R")
	# Simulation function physig_sim
	# input: x = number of species in tree
	# output: a vector length two with (K, lamba)
	physig_sim <- function(x) {
	tree <- rcoal(x)
	traits <- fastBM(tree)
	physig_k <- phylosig(tree, traits, method = "K")
	physig_l <- phylosig(tree, traits, method = "lambda")$lambda
	sigs <- c(physig_k, physig_l)
	return(sigs)
	}
	# Run simulation
	spnumbs <- c(5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 65, 80, 95)
	rands <- 1000
	siglist <- alply( spnumbs, 1, function(x) t(replicate(rands, physig_sim(x))), .progress="text")
	sigdf <- as.data.frame(do.call(rbind, siglist))
	names(sigdf) <- c("K", "L")
	sigdf$numsp_ <- rep(spnumbs, each=rands)
	sigdf_m <- melt(sigdf, id = 3)
	# Plot results
	plotdf <- ddply(sigdf_m, .(numsp_, variable), summarise,
	mean = mean(value),
	se = sd(value)/sqrt(length(value))
	)
	limits <- aes(ymax = mean + se, ymin = mean - se)
	dodge <- position_dodge(width=0.9)
	ggplot(plotdf, aes(x = numsp_, y = mean, shape = variable)) +
	geom_point(position=dodge) +
	geom_errorbar(limits, position=dodge, width=0.25) +
	geom_smooth()
	ggsave("physig_sim.jpeg")

view raw phylogeneticsignal_sim.R hosted with ❤ by GitHub