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.