Noam Ross

Improved R Profiling Summaries

2013-05-02T16:48:01-07:00

In my last post I mentioned that I had improved on R’s summaryRprof() function with a custom function called proftable(). I’ve updated proftable() to take advantage of R 3.0.0’s ability to record line numbers while profiling. I’ve put it on github – you can get it there or below.

proftable reads in a file generated by Rprof() and creates an easy-to read table of the most time-consuming calls in your code, ordered from most time-consuming to least. Unlike summaryRprof(), it prints the entire call stack, so you can trace the origin of the time hogs. To make this easier to read, I lop off the “parent stack” common to all of the function calls, and display it just once, below the table. Here’s some example output:

> Rprof("profile1.out", line.profiling=TRUE)
> source("http://pastebin.com/download.php?i=KjdkSVZq")
> Rprof(NULL)
> proftable("profile1.out", lines=10)

 PctTime Call                                                      
 20.47   1#17 > [ > 1#17 > [.data.frame                            
  9.73   1#17 > [ > 1#17 > [.data.frame > [ > [.factor             
  8.72   1#17 > [ > 1#17 > [.data.frame > [ > [.factor > NextMethod
  8.39   == > Ops.factor                                           
  5.37   ==                                                        
  5.03   == > Ops.factor > noNA.levels > levels                    
  4.70   == > Ops.factor > NextMethod                              
  4.03   1#17 > [ > 1#17 > [.data.frame > [ > [.factor > levels    
  4.03   1#17 > [ > 1#17 > [.data.frame > dim                      
  3.36   1#17 > [ > 1#17 > [.data.frame > length                   

#File 1: http://pastebin.com/download.php?i=KjdkSVZq

Parent Call: source > withVisible > eval > eval >

Total Time: 5.96 seconds
Percent of run time represented: 73.8 %

Note that the “Parent Call” at the bottom shows some functions which RStudio wrapped my code in. Also, I chose only to display the top 10 time-consuming calls, but proftable told me that those 10 calls represent 73.8% of the run time. I find this display a lot more intuitive than summaryRprof()

Here’s the whole function. If you have improvements, fork it on github:

proftable <- function(file, lines=10) {
# require(plyr)
  interval <- as.numeric(strsplit(readLines(file, 1), "=")[[1L]][2L])/1e+06
  profdata <- read.table(file, header=FALSE, sep=" ", comment.char = "",
                         colClasses="character", skip=1, fill=TRUE,
                         na.strings="")
  filelines <- grep("#File", profdata[,1])
  files <- aaply(as.matrix(profdata[filelines,]), 1, function(x) {
                        paste(na.omit(x), collapse = " ") })
  profdata <- profdata[-filelines,]
  total.time <- interval*nrow(profdata)
  profdata <- as.matrix(profdata[,ncol(profdata):1])
  profdata <- aaply(profdata, 1, function(x) {
                      c(x[(sum(is.na(x))+1):length(x)],
                        x[seq(from=1,by=1,length=sum(is.na(x)))])
              })
  stringtable <- table(apply(profdata, 1, paste, collapse=" "))
  uniquerows <- strsplit(names(stringtable), " ")
  uniquerows <- llply(uniquerows, function(x) replace(x, which(x=="NA"), NA))
  dimnames(stringtable) <- NULL
  stacktable <- ldply(uniquerows, function(x) x)
  stringtable <- stringtable/sum(stringtable)*100
  stacktable <- data.frame(PctTime=stringtable[], stacktable)
  stacktable <- stacktable[order(stringtable, decreasing=TRUE),]
  rownames(stacktable) <- NULL
  stacktable <- head(stacktable, lines)
  na.cols <- which(sapply(stacktable, function(x) all(is.na(x))))
  stacktable <- stacktable[-na.cols]
  parent.cols <- which(sapply(stacktable, function(x) length(unique(x)))==1)
  parent.call <- paste0(paste(stacktable[1,parent.cols], collapse = " > ")," >")
  stacktable <- stacktable[,-parent.cols]
  calls <- aaply(as.matrix(stacktable[2:ncol(stacktable)]), 1, function(x) {
                   paste(na.omit(x), collapse= " > ")
                     })
  stacktable <- data.frame(PctTime=stacktable$PctTime, Call=calls)
  frac <- sum(stacktable$PctTime)
  attr(stacktable, "total.time") <- total.time
  attr(stacktable, "parent.call") <- parent.call
  attr(stacktable, "files") <- files
  attr(stacktable, "total.pct.time") <- frac
  cat("\n")
  print(stacktable, row.names=FALSE, right=FALSE, digits=3)
  cat("\n")
  cat(paste(files, collapse="\n"))
  cat("\n")
  cat(paste("\nParent Call:", parent.call))
  cat(paste("\n\nTotal Time:", total.time, "seconds\n"))
  cat(paste0("Percent of run time represented: ", format(frac, digits=3)), "%")

  invisible(stacktable)
}

FasteR! HigheR! StrongeR! - A Guide to Speeding Up R Code for Busy People

2013-04-25T14:49:34-07:00

This is an overview of tools for speeding up your R code that I wrote for the Davis R Users’ Group.

First, Ask “Why?”

It’s customary to quote Donald Knuth at this point, but instead I’ll quote my twitter buddy Ted Hart to illustrate a point:

I’m just going to say it.I like for loops in #Rstats, makes my code readable.All you [a-z]*ply snobs can shove it!
— Ted Hart (@DistribEcology) March 12, 2013

Code optimization is a matter is a matter of personal taste and priorities. There may be some ways of writing code that are better or worse, and there are definitely ways that will make it run faster, but before you dive into optimization, you should ask yourself these questions:

(1) Do you want your code to be readable?

If you need to explain your code to yourself or others, or you will need to return to it in a few months time and understand what you wrote, it’s important that you write it in a way that is easy to understand. Comments, indenting, and formatting go a long way, but your choices make a difference, as Ted notes above.

Some optimal code can be hard to read. For instance, and curly braces can be faster than parentheses. But which is more intuitive to read?¹

(2) Do you want your code to be sharable?²

Most of the considerations of (1) apply here, but they have to be balanced with the fact that, if your code is painfully slow, others are not going to want or have time to use it.

Also, some optimization strategies don’t transfer well. If you use parallel processing, it won’t be usable to others who don’t have multicore computers. Some tricks for speeding up code will work well as a quick-and-dirty fix, but are likely to break in new versions of R or are disallowed from CRAN packages altogether.

(3) Do you have other shit to do?

No? Please contact me and I’ll help you with that. Yes? You are among the 90% or R users who’s first priority is not computer programming. The time spent optimizing code is often longer than the computing time actually saved. Use the simple solutions and you can get on with your research/life.

With these in mind, I’ll go through some strategies for speeding up your R computations. Because I’m especially mindful of consideration #3, this guide is somewhat in the reverse of traditional tutorials on this topic. I’m going to talk about the blunt instruments first, then get into the idiosyncrasies of code.

Get a Bigger Computer

If you have a task that only needs to be accomplished a few times, you might be best off finding more resources rather than re-writing your code. You have an Amazon.com account, right? Amazon rents virtual machines on the web. It’s so easy to set one up that I accidentally started one today without meaning to while searching for the link. The cost of computing power roughly equivalent to my 2011 MacBook Pro is 13¢ per hour. You can have 100 MacBooks of computing power for less than minimum wage.³ These virtual machines can also have much more memory than your laptop or desktop.

Sometimes you don’t even need more computing power, you just need to run R on some other computer without tying up your machine for hours. In any case, consider this option. If you just need to run a task once, and you think you are going to have to spend more than an hour rewriting your code, you might as well use that hour getting a virtual machine set up and never worry about computing power again.

Cloud resources like Amazon aren’t the only way to get more computer power. If you’re at a university, there are probably high-power computing clusters available to you. Depending on your particular set-up, you may need to to run your program in batch mode to use such clusters. Personally, I think that a cloud-based service like Amazon avoids administrative hassle if you are on a short deadline or have a one-off task, but if you need lots of power regularly, getting in on a local cluster is worthwhile.

The most user-friendly way to use cloud computing is to set up a server to run RStudio, which can then be accessed by a web browser and works just like it would on your desktop. Bioconductor has a step by step guide to starting an Amazon Web Instance that they’ve pre-configured this way. There are a couple of other useful guides to doing the same thing by Louis Aslett (including a video) and Karthik Ram. Read all of these and you should be able to set yourself up in an hour or two.

A fancier alternative is the experimental Segue package, which provides an apply() function variant, emrapply(), that sends commands from your local R session to Amazon servers and returns the results. An introductory guide is here.

Of course, if you are writing software that you expect others to run on normal laptop and desktop computers, you’ll want to make sure it does so at a reasonable speed. In this case, running your software in the cloud may not make sense.

Parallel Processing

Parallel processing used to be a very fancy thing, but with recent R packages, it can be implemented with very little programming time. Parallel processing breaks up your task, splits it among multiple processors, and puts the components back together. This is very useful and easy if you have a task that can be split up, especially without the different parts needing to “talk to” each other.

On your typical laptop or desktop, implementing parallel processing in your R code might speed up your program by a factor of 2 to 4. However, if you have a big machine, or you followed my advice in the last section, you can take advantage of large computers with many CPUs or “cores”.

I’m not going to get into the nitty gritty of setting up parallel processes and the various packages that do so e.g. multicore, snow,snowfall, doParallel,foreach, plyr,…) because Jacob Teter will give a presentation on that at D-RUG next week. Instead I’ll discuss a bit about when this approach might be appropriate.

Essentially, parallel processing works best when you have a task that has to be completed many times, but each repeat is independent. For instance, if you are repreating a simulation, and in each case are drawing new parameters from a distribution, that is an easily parallelizable task. However, if you are simulating population dynamics, in which each time step you need the outputs of the previous time steps, than you can’t split up each time step into parallel components.

A good rule of thumb is that if you can wrap your task in an apply function or one of its variants, it’s a good option for parallelization. In fact most implementations of parallel processing in R are versions of apply.

Not every task runs better in parallel. When you run a task in parallel, your computer “dispatches” each task to a CPU core. This dispatching adds computational overhead. So, it’s usually best to try to minimize the number of dispatches. In most cases you are going to have a small number of computing cores relative to your tasks. A powerful laptop or desktop will have 2, 4, or 8 cores, and even the most powerful Amazon virtual machine has 88.

The following would be a poor use of parallel computing:

aaply(1:10000, 1, function(x) rnorm(1, mean=x), .parallel=TRUE)

This would dispatch the very small task of generating a random number 10,000 times. The dispatching would probably take more time than the actual calculations. Instead, organize the tasks into blocks, and dispatch each block of tasks, for instance:

aaply(seq(1,10000,100), function(x) rnorm(1, mean=x+(1:100), .parallel=TRUE)

This sends the task of generating 100 numbers to each core, reducing the overhead.

[Update 4/26/2013: As Vince Buffalo points out, random number generation requires that R keep track of state, so may not be a good example of a parallel-izable task]

In practice, this means that you want to “parallelize” your code at the highest, or chunkiest level. If you have a series of nested loops, for instance, you will want to run the highest-level loop in parallel.

Note that more cores is not always better. For reasons I don’t completely understand, the overhead of using more cores is not always worth it, so it may be better to use, say 16, rather than 64 codes. Some experimentation will help here.

One other thing: if you can limit the parallel processes to computational tasks, and avoid things like reading/writing to disk or downloading files in the parallel code, you’ll be less likely to run into problems.

A disadvantage of parallel computing is that it is not available to everyone. People run code on different machines. If you are writing software that you expect to share with others (and as a scientist, you should ALWAYS be prepared to share your work on request, and post your code as a supplement to your paper), you will want to provide an option to run it not in parallel. I like to do something like this:

#This script uses parallel processing if p.flag=TRUE.  Set up a parallel 
#cluster as appropriate for your machine as appropriate. (the commented code 
#below will use 2 cores on a multicore computer)
#library(doParallel)
#cl <- makeCluster(2)  # Use 2 cores
#registerDoParallel(cl) # register these 2 cores with the "foreach" package
library(plyr)
p.flag=FALSE  # Change to TRUE if using parallel processing
.
.
.
 aaply(seq(1,10000,100), function(x) rnorm(1, mean=x+(1:100), 
      .parallel=p.flag)

Compiling Functions

R code is interpreted when it is run, unlike some other languages, which are first compiled. This is one reason why functions written in C are often faster than functions written in R.

However, R has compiling ability, which can speed up functions by a factor or 3 or 4 in some cases. Compiling functions just requires calling cmpfun() in the base compiler package. Here’s an example:

library(compiler)
f <- function(n, x) for (i in 1:n) x = (1 + x)^(-1)
g <- cmpfun(f)

library(microbenchmark)
compare <- microbenchmark(f(1000, 1), g(1000, 1), times = 1000)

library(ggplot2)
autoplot(compare)


compare

## Unit: microseconds
##        expr   min    lq median     uq  max neval
##  f(1000, 1) 950.0 991.3 1024.7 1087.2 2799  1000
##  g(1000, 1) 218.2 228.0  238.2  257.3 1743  1000

(See the benchmarking section below for more on comparing functions with microbenchmark())

cmpfun() works well in cases like this - when you are defining a new function that is mostly numerical manipulation. It doesn’t help with functions that call a lot of other R functions, or involve manipulating and translating between data types. It also doesn’t work well on already-defined R functions, which have been pre-compiled since R 2.14 and are sometimes written in compiled languages like C. For instance, if I try to compile this paste()-based function, I get no improvement:

f <- function(n) paste(1:n, collapse = " ")
g <- cmpfun(f)
compare <- microbenchmark(f(1000), g(1000), times = 1000)
autoplot(compare)

compare

## Unit: microseconds
##     expr   min    lq median    uq  max neval
##  f(1000) 320.2 331.8  353.9 429.2 1889  1000
##  g(1000) 320.7 330.6  352.8 412.5 2530  1000

You can also enable just-in-time compiling in R, which will automatically compile EVERY function the first time that it is run. This will slow down R quite a bit at first, as each function must be compiled before it is run the first time, but but then speed it up later. Just add this to the start of your script

library(compiler)
enableJIT(1)

The numeric argument in enableJIT specifies the “level” of compilation. At level 2 and 3 more functions are compiled. In a few rare cases enableJIT will slow down code, particularly, if most of it is already pre-compiled or written in C, and/or if the code creates functions repeatedly, which then need to be compiled every time. This is more likely to happen with enableJIT(2) or enableJIT(3), though these have the potential to speed up code more, as well.

Diagnosis: Profiling and Benchmarking

OK, the tools I described above are all blunt instruments that could speed up many types of code. However, if they fail you or aren’t feasible, you’re going to need to know something about what is slowing your code down.

Benchmarking

Benchmarking is a way to test and compare the speed of functions.

Traditionally, system.time(replicate()) could be used to time functions, but the microbenchmark package is a better option for benchmarking individuals calls and functions. If you have identified a function that is called many times in your code and needs to be speeded up, you can write or try out several different versions, using microbenchmark to compare them. Here’s the example from above:

f <- function(n, x) for (i in 1:n) x = (1 + x)^(-1)
g <- cmpfun(f)

library(microbenchmark)
compare <- microbenchmark(f(1000, 1), g(1000, 1), times = 1000)

library(ggplot2)
autoplot(compare)

compare

## Unit: microseconds
##        expr   min     lq median     uq   max neval
##  f(1000, 1) 958.2 1064.2 1220.5 1795.7 11595  1000
##  g(1000, 1) 219.8  243.7  271.4  379.4  1370  1000

All microbenchmark needs as inputs are a series of expressions to compare. times is an optional parameter of the number of times to test each function (default 100). Microbenchmark randomizes the test order and discards some “warmup” tests. The full distribution is returned as a data frame, but only the summaries are printed.

The autoplot() command draws a comparison of distribution times using ggplot()

Profiling

Benchmarking is good for comparing small functions. Profiling examines your code to determine what parts of it are running slow. This is the essential tool for getting the best bang for your buck in speeding up code.

R’s code profiler samples the “call stack” at regular intervals. The call stack is the list of the current function running, the function that called it, and all the functions that those. By examining these samples, you can find which functions are taking up the largest amount of time, either because they are slow or are called many times.

You activate the profiler with Rprof("filename") which writes each sample of the call stack to “filename”. So do this:

Rprof("file.out")
Some_slow_code()
Rprof(NULL)
summaryRprof("file.out")

Rprof(NULL) stops sampling and summaryRprof("file.out") gives you a summary of the results. Here’s an example of profiling a random walk function (from a talk by Duncan Temple Lang):

rw2s1 <- function(n) {
    xpos = ypos = numeric(n)
    xdir = c(TRUE, FALSE)
    pm1 = c(1, -1)
    for (i in 2:n) if (sample(xdir, 1)) {
        xpos[i] = xpos[i - 1] + sample(pm1, 1)
        ypos[i] = ypos[i - 1]
    } else {
        xpos[i] = xpos[i - 1]
        ypos[i] = ypos[i - 1] + sample(pm1, 1)
    }
    list(x = xpos, y = ypos)
}

Rprof("out.out")
for (i in 1:1000) pos = rw2s1(1000)
Rprof(NULL)
summaryRprof("out.out")

## $by.self
##              self.time self.pct total.time total.pct
## "sample.int"      7.62    47.04       7.62     47.04
## "rw2s1"           5.28    32.59      16.20    100.00
## "sample"          2.90    17.90      10.52     64.94
## "-"               0.30     1.85       0.30      1.85
## "+"               0.04     0.25       0.04      0.25
## "list"            0.04     0.25       0.04      0.25
## "numeric"         0.02     0.12       0.02      0.12
## 
## $by.total
##              total.time total.pct self.time self.pct
## "rw2s1"           16.20    100.00      5.28    32.59
## "sample"          10.52     64.94      2.90    17.90
## "sample.int"       7.62     47.04      7.62    47.04
## "-"                0.30      1.85      0.30     1.85
## "+"                0.04      0.25      0.04     0.25
## "list"             0.04      0.25      0.04     0.25
## "numeric"          0.02      0.12      0.02     0.12
## 
## $sample.interval
## [1] 0.02
## 
## $sampling.time
## [1] 16.2

The two tables list the functions called. The $by.self component lists the time spent by functions alone, while the $by.total table lists the time spent by functions and all the functions they call. So you see that rw2s1, being the parent function, takes the most time in $by.total, but sample.int is the bottom-level function that is taking the most time.

I don’t find this summaryRprof() output very intuitive. There are some moderately better tools out there in the packages profr and proftool, but I don’t like them either. Instead I’ve written a custom function that you can get here called proftable. Here it parses the same output as above:

proftable("out.out")

##  PctTime Stack                      
##  47.037  rw2s1 > sample > sample.int
##  32.593  rw2s1                      
##  17.901  rw2s1 > sample             
##   1.852  rw2s1 > -                  
##   0.247  rw2s1 > +                  
##   0.247  rw2s1 > list               
##   0.123  rw2s1 > numeric            
## 
## Total Time: 16.2 seconds

Here you can see that >50% of time is being taken up by sample.int, which is itself usually called by sample.

Misc profiling things

The latest version of R (3.0.0) allows code profiling by line number using Rprof(file, line.profiling=TRUE) and summaryRprof(file, lines="show"). I don’t yet have a good example of it and haven’t incorporated this function into proftable() yet. There’s an example of use at Stack Overflow
Note that, depending on the R interface you are using, it may wrap your code in some other functions, which will appear as the topmost functions in the stack. For instance, if use the “Source” button in Rtudio, all your stacks will start with something like:
```
source > withVisible > eval > eval >
```
[Update 5/8/2003: The latest verion of proftable can deal with line numbers. It also deals with wrapper stacks like above.]
trace() can also be useful to track down where a particular slow function is being called. Try using trace(function, tracer=cat(as.character(sys.calls())), for instance, to print the stack whenever function is called.
R also has the ability to do memory profiling - tracking how your code uses memory. These include Rprofmem(), tracemem(), and Rprof(memory.profiling=TRUE), but the documentation for these isn’t great and I haven’t found good resources on this topic. Let me know if you can point me to some resources.
Jeromy Anglim has a blog post of an Rprof() case study that is worth reading

Find better packages

Before you go about changing your own code, take a look and see if you can use better code someone else has written. There are thousands of packages on CRAN, but many of them are redundant and implement the same thing. There’s tremendous variation in package quality and speed, so you can often speed things up by finding a faster package.

For instance, the xts and zoo packages both handle time series data, but much of xts is written in C and Fortran, and it is considerably faster. Here are some performance benchmarks (though they are a bit out of date):

Comparison of functions on xts time series objects against time series from other packages. From this presentation

Other packages are specifically designed to speed up certain operations. Some creative googling (e.g., r package CRAN fast “something”) will help you find them.

For instance. Using base functions like colMeans() is faster than using apply functions to apply mean. The matrixStats package extends these functions to include others like rowMaxs(), colRanks(), etc.

data.table is a package designed for speeding up data frame operations when working very large datasets. data.table’s biggest advantage is in tasks like sub-setting and indexing. Here I compare the two in looking up the value of a row:

library(data.table)
myDF <- as.data.frame(matrix(rnorm(1e+05), 10000, 10))
myDF$key <- as.character(1:nrow(myDF))
myDT <- data.table(myDF)
setkey(myDT, key)
microbenchmark(myDF[myDF$key == 2, ], myDT["2", ])

## Unit: microseconds
##                   expr    min     lq median   uq  max neval
##  myDF[myDF$key == 2, ] 2658.8 2885.3   3187 3755 6434   100
##            myDT["2", ]  724.3  899.7   1112 1393 4287   100

Here data.table is 3X times faster. It gets faster with larger data sets.

Improving your code

Once you know what parts of your code are slowing you down, and you have a method to test speed, you need tools and options for writing faster functions. Here are a few.

Vectorization

This is the first thing most guides to coding discuss, so I won’t go into in much. Basically, R is good at matrix/vector algebra, so if you can express something as

It will go faster than

\[\begin{aligned} x + a \\ y + b \\ z + c \\ \end{aligned}\]

Or, in R

x <- c(1,2,3) + c(1, 2, 3)

Instead of

for(i in 1:3) x[i] <- i+i

Most base R functions are designed to handle vectors as well as scalars.
However, using an apply function is not really vectorization. apply functions are actually wrappers for for loops
For more on vectorization, see Chapters 3 and 4 of The R Inferno

Pre-allocating memory

Again, a common topic I won’t dwell on too much. Essentially, R is bad at continually re-sizing objects, because it makes an extra copy of these objects each time. So if you have a loop that creates a vector or list, don’t append to the vector or list with each pass of the loop. Instead, make an empty object of the correct size first, then fill in its elements.

More here

The apply family of functions helps you avoid these problems by doing pre-allocation automatically. Plus, as I wrote above, if you use apply variants, it’s easy to parallelize your code.

Using simpler data structures

R has some very handy objects, lists and dataframes, that are flexible and can store multiple types of data at once. However, this flexibility comes at a cost. Simpler data structures that only store one type of data can be manipulated much faster. If you can represent data in a matrix instead of a data frame (e.g., if you can reduce it to all numbers or characters), you can speed things up considerably.

For instance, compare these operations for getting the standard deviation of 100,000 rows of 10 variables in both a matrix and a data frame:

myMatrix <- matrix(rnorm(1e+05), 10000, 100)
myDF <- as.data.frame(myMatrix)
microbenchmark(apply(myMatrix, 2, sd), apply(myDF, 2, sd))

## Unit: milliseconds
##                    expr   min    lq median    uq    max neval
##  apply(myMatrix, 2, sd) 29.70 31.60  33.04 34.65  94.21   100
##      apply(myDF, 2, sd) 39.69 70.14  96.61 98.81 250.12   100

The small stuff

Little things can add up, but the little things are most likely not your main problem. Here are a few odds and ends I’ve discovered. They’ll shave a few percent off your run time, but thinking about them too much will add hours to your coding time.

Extra parentheses

Sometimes when I write complex mathematical expressions I include extra parentheses to help make the code more readable. This has a small performance cost, as you can see here:

a <- function(x) x = 1/{
    1 + x
}
f <- function(x) x = 1/(1 + x)
g <- function(x) x = (1/(1 + x))
h <- function(x) x = ((1/(1 + x)))
k <- function(x) x = (((1/(1 + x))))
x <- sample(1:100, 100, replace = TRUE)
comp <- microbenchmark(a(x), f(x), g(x), h(x), k(x), times = 1e+05)
comp

## Unit: microseconds
##  expr   min    lq median    uq   max neval
##  a(x) 1.991 2.197  2.313 2.971 11068 1e+05
##  f(x) 1.960 2.174  2.292 2.990 10435 1e+05
##  g(x) 2.090 2.302  2.427 3.112 10755 1e+05
##  h(x) 2.188 2.406  2.536 3.178  9660 1e+05
##  k(x) 2.259 2.504  2.644 3.367 11228 1e+05

autoplot(comp)

The difference is modest, but each layer of parentheses adds a bit to performance time. Note that, amazingly, r speeds up when you use curly braces ({ and }) instead of parentheses! Why? Radford Neal explains:

The difference between parentheses and curly brackets comes about because R treats curly brackets as a “special” operator, whose arguments are not automatically evaluated, but it treats parentheses as a “built in” operator, whose arguments (just one for parentheses) are evaluated automatically, with the results of this evaluation stored in a LISP-style list. Creating this list requires allocation of memory and other operations which seem to be slow enough to cause the difference, since the curly bracket operator just evaluates the expressions inside and returns the last of them, without creating such a list.

Standard and internal functions

A few standard functions in R are pretty slow, often because they perform error checking and methods dispatch before actually doing arithmetic. mean is among the most notorious of these. mean(x)does a lot more than just sum(x)/length(x). It checks the form of the inputs, assigns the appropriate method, then dispatches .Internal(mean(x)). The latter is much faster, as it is written in C:

x <- sample(1:100, 100, replace = TRUE)
comp <- microbenchmark(mean(x), sum(x)/length(x), .Internal(mean(x)), times = 1e+05)
comp

## Unit: nanoseconds
##                expr  min   lq median   uq      max neval
##             mean(x) 5671 6136   6386 6792 61466725 1e+05
##    sum(x)/length(x)  720  878   1064 1211  1511792 1e+05
##  .Internal(mean(x))  513  556    614  677  1082112 1e+05

autoplot(comp)

You can sometimes get away with using the internal functions like .Internal(mean(x)). However, this is considered poor practice because (a) internals are subject to change and break your code with updates to R, and (b) this is disallowed in packages submitted to CRAN, so it makes code more difficult to share. That said, the internal version can be 9X faster. Your call.

[Update 4/26/2013: On the suggestion of Kevin Ushey, I’m adding the following section]

Other Languages: More reward for more effort

I’ve mentioned a few times that functions written in other languages, like C or C++, are much faster than other R functions. Obviously, learning another language is an undertaking, but thanks to some recent developments, writing R functions in C has become much easier. The Rcpp package lets you write such C++ functions without learning much about compiling or other peripheral issues, and Hadley Wyckham has written an excellent beginners’s guide to Rcpp that teaches you just enough C++ to get started. From the guide:

Typical bottlenecks that C++ can help with are:

Loops that can’t easily be vectorised because each iteration depends on the previous. C++ modifies objects in place, so there is little overhead when modifying a data structure many times.

Recursive functions, or problems which involve calling functions millions of times. The overhead of calling a function in C++ is much lower than the overhead of calling a function in R. To give you some idea of the magnitude, on my computer when writing this book the overhead in C++ was ~5ns compared to ~200ns for R.

Problems that require advanced data structures and algorithms that R doesn’t provide. Through the standard >template library (STL), C++ has efficient implementations of many important data structures, from ordered maps >to double ended queues.

I used to also have a reference here that x*x is faster than x^2 in R, but as Thomas Lumley points out, this hasn’t been true for a while, and some quick empirical testing confirmed that the two are pretty close.↩
If you are a scientist, you should ALWAYS be prepared to share your code, and should probably publish it as a supplement to your papers.↩
There are cheaper services than Amazon, but there aren’t as many resources on how to use them. If you really need to reduce costs, you can use spot instances, which allow you to set an allowable price, and run your program when something that cheap is available. Also, you could run many small Amazon virtual machines simultaneously on a free trial account, but if your goal is to do things quickly, it might be worth paying for a more powerful machine.↩

Setting up Automated Reprint Requests in Gmail

2013-04-24T17:37:22-07:00

I recently wrote about my first peer-reviewed publication and lamented that it’s behind a $1500 paywall. I am not allowed to post the file to my website, but I am able to distribute the paper by e-mail, so I’ve set up a system to automate requests for the paper.

Michael McCarthy wrote a similar post this fall for people using Outlook e-mail. My method is similar but is meant for Gmail. Here’s how to do it:

Generate a custom mailto link

I used this website to generate a custom mailto link that would open the users’ mail client and send a pre-formatted e-mail to me. The email template is:

To: noam.ross@gmail.com

Subject: Reprint Request: Economics and Ecology of Open-Access Fisheries

Body:

Dear Noam,

Please send me a reprint of:

Fuller, Kate, David Kling, Kaeilin Kroetz, Noam Ross, and James N. Sanchirico (2013) Economics and Ecology of Open-Access Fisheries. In Encyclopedia of Energy, Natural Resource, and Environmental Economics, Ed. Jason Shogren, p.39-49. Amsterdam: Elsevier. http://dx.doi.org/10.1016/B978-0-12-375067-9.00114-5

This generates the link:
```
mailto:noam.ross@gmail.com?Subject=Preprint%20Request%3A%20Economics%20and%20Ecology%20of%20Open-Access%20Fisheries&Body=Dear%20Noam%2C%0A%0APlease%20send%20me%20a%20reprint%20of%3A%0A%0AFuller%2C%20Kate%2C%20David%20Kling%2C%20Kaeilin%20Kroetz%2C%20Noam%20Ross%2C%20and%20James%20N.%20Sanchirico%20%282013%29%20Economics%20and%20Ecology%20of%20Open-Access%20Fisheries.%20In%20%20*Encyclopedia%20of%20Energy%2C%20Natural%20Resource%2C%20and%20Environmental%20Economics*%2C%20Ed.%20Jason%20Shogren%2C%20p.39-49.%20Amsterdam%3A%20Elsevier.%20http%3A//dx.doi.org/10.1016/B978-0-12-375067-9.00114-5
```
Put the link on your webpage in lieu of a link to a PDF

Here it is on my site:

Put the PDF in your Dropbox and get a link to the file

Unfortunately Gmail won’t created automated responses with file attachments¹, so you still need to create a link where the file can be downloaded, even if it’s not public. I use Dropbox:

Create a Gmail template

Gmail has an experimental feature in “labs” that creates templates. Turn it on:

Then write an email without a recipient and save it as a canned response:

Create an email filter

Create a filter specific to the custom mailto link you created (e.g., the same subject line). With canned responses enabled, you can choose “Send canned response” as an option on gmail filters. Select the canned response you saved. I also add a tag so I can track these requests easily.

That’s it! Thanks to Michael McCarthy for the idea.

I tried to get automated attachments working with this IFTTT recipe, but the reply always had mangled text. If someone gets it working, let me know.↩

My first peer-reviewed publication: Economics and Ecology of Open-Access Fisheries

2013-04-24T10:57:21-07:00

My first peer-reviewed publication is an encyclopedia article written with colleagues in economics in UC Davis:

Economics and Ecology of Open-Access Fisheries

K. Fuller, D. Kling, K. Kroetz, N. Ross and J.N. Sanchirico (2013)

In: Shogren, J.F., (ed.) Encyclopedia of Energy, Natural Resource, and Environmental Economics, Vol. 2 Encyclopedia of Energy, Natural Resource, and Environmental Economics p.39-49. Amsterdam: Elsevier. http://dx.doi.org/10.1016/B978-0-12-375067-9.00114-5

Abstract

In open-access fisheries, participation and harvest by fishermen are unregulated. Under these conditions, the amount of fishing is greater than the socially optimal level, the economic profits from fishing are dissipated, and fish stocks are depleted and can even be driven to extinction. A large body of theoretical and empirical research has demonstrated the robustness of these results.

Although the basic outcomes hold across a variety of systems, the behavior of an open-access fishery can be complex. Simple models of fish ecology coupled with harvester behavior can explain boom-and-bust cycles in a fishery and fishing communities and complex supply responses. More complex ecological models of open-access fisheries illustrate the interactions among fishing, fish populations, and other parts of the ecosystem, as well as explain why some fishing regulations produce unintended consequences. Studying open access provides a foundation for understanding fishery economics, a baseline from which to measure the gains of implementing fishing regulations, and lends insight to the design of better regulatory systems.

The process of writing this article was interesting. Two years ago, when I was a first-year graduate student, my economics advisor Jim Sanchirico began convening lab meetings. This isn’t common practice among economists as far as I know, and we were trying to figure out how to use the meetings constructively. Around the same time, Jim received a request to write a few articles for this encyclopedia, and we figured it would be a good exercise to read the original literature on this subject and write the articles together.

This was a great chance for me. I had zero exposure to the natural resource economics literature and wouldn’t get around to courses in it for another year. Also, as a first-year ecology grad student working with a group of more experienced economists, I was a pretty good proxy for the target audience of outsiders to the field. So I got play the role of naive reader, asking for clearer explanations and less jargon from my co-authors. I also wrote some of the introductory and ecological sections.

I’m grateful to Jim for putting me on board this project, and for the chance to work with great colleagues in a different field, but I probably won’t do something similar again. I don’t support the publication model. Ironically, while the article is about “Open-Access Fisheries,” it is about as far from an open-access publication as can be. It’s in a $1500 encyclopedia published by Elsevier. We received a very small honorarium¹, and Elsevier will reap profits by selling access to this encyclopedia to libraries as part of some overpriced bundle. I’d rather spend my taxpayer-funded fellowship creating resources available to all.

If you’d like a copy, email me by clicking here

If that’s what it is called. The arrangement was, I think, one or two hundred dollars per article, which the lab was supposed to split for dinner and beer. Still waiting on that part, Jim.↩

Debugging Tools in R with Michael Hannon

2013-04-18T13:44:59-07:00

Today at Davis R Users’ Group, Michael Hannon gave a great talk on how to use R’s native debugging functions. Here are his notes and code.

Introduction

This is a discussion of debugging techniques in R. It is based on a paper by Roger Peng, now at Johns Hopkins University (http://www.biostat.jhsph.edu/~rpeng/docs/R-debug-tools.pdf)

Focus on five functions:

traceback
debug
browser
trace
recover

Severity level

R mainly uses two ways of reporting a problem:

Warning: Does not stop execution
Error: More serious – does stop execution

Example of warning:

> zzz <- log(-1)

## Warning: NaNs produced

> zzz

## [1] NaN

Example of error:

#### Define a function to be used/abused later on.
message <- function(x) {
  if(x > 0)
    print("Hello")
  else
    print("Goodbye")
  }

Here’s the function in action:

x <- 42
message(x)

## [1] "Hello"


x <- log(-1)  #### Produces warning but does not stop

## Warning: NaNs produced

message(x)

## Error: missing value where TRUE/FALSE needed

Debugging tools

Print “call stack” with `traceback`

How did we get here? That is, what sequence of function calls led us to the point where an error occurred? Answer: the traceback function can tell us. Here’s an example:

>   message(log(-1))
Error in if (x > 0) print("Hello") else print("Goodbye") : 
  missing value where TRUE/FALSE needed
In addition: Warning message:
In log(-1) : NaNs produced
> 
>   traceback()
1: message(log(-1))

It’s hard to get “lost” when you’ve got only one function to begin with, but real life is more complicated. THe following isn’t exactly “real-life” either, but it better illustrates the point.

f <- function(x) {
    r <- x - g(x)
    r
}
g <- function(y) {
    r <- y * h(y)
    r
}
h <- function(z) {
    r <- log(z)
    if (r < 10) 
        r^2 else r^3
}

You can probably guess that the log function is going to be the trouble maker, but let’s see what happens.

> f(-1)
Error in if (r < 10) r^2 else r^3 : missing value where TRUE/FALSE needed
In addition: Warning message:
In log(z) : NaNs produced
> traceback()
3: h(y) at #2
2: g(x) at #2
1: f(-1)
>

This shows that we started out with a call to f, which generated a call to g, which generated a call to h. The error evidently occurred at statement #2 of h. Note, by the way, the following trick to view the line numbers in a function:

> as.list(body(h))
[[1]]
`{`

[[2]]
r <- log(z)

[[3]]
if (r < 10) r^2 else r^3

I.e., statement #2 of h is, as expected, the one that involves the log function.

The `debug` function

The traceback function helps us to locate the error, but then what do we do? Ideally, we could just stare at our code, then find and fix the error. But it often helps to have some detailed information as to what was going on in our code at the time of the error.

The debug function allows us to do just that. We can step through our code, a line at a time, and examine variables along the way.

The syntax is:

debug(function-name)

and to turn off debugging we use:

undebug(function-name)

(I.e., we might find and fix a problem on the third pass through a loop of 10000 iterations, after which we no longer need or want to debug for the other 9997 times through the loop.)

Here’s an example. The mathematical background is that we have N observations of a random variable that we assume to be normally distributed but with an unknown mean and variance. To estimate the parameters we need to compute the sum of the squared differences from the mean:

\[ SS = \sum_{n=1}^{N}(x_n - \mu)^2 \]

(Details at: http://en.wikipedia.org/wiki/Normal_distribution#Estimation_of_parameters)

Here we’ll focus on just one part of the calculation, computing the sum of the squared differences. Here’s an R function that will do that for us:

SS <- function(mu, x) {
    d <- x - mu
    d2 <- d^2
    ss <- sum(d2)
    ss
}

You can probably see that the function could be written much more compactly, but we’ll string it out for purposes of illustration.

Let’s generate some points from a Normal(0,1) distribution:

set.seed(314159)  # Set the RNG seed so results are reproducible
xPoints <- rnorm(100)

Let’s run SS to get a feel for its behavior:

SS(1, xPoints)

## [1] 202.5

muTest <- as.numeric(-5:5)
SSTest <- sapply(muTest, SS, xPoints)
plot(muTest, SSTest)

Now suppose we want to examine in detail the behavior of the SS function. First, we flag the function for debugging:

debug(SS)

> SS(1, xPoints)
debugging in: SS(1, xPoints)
debug at #1: {
    d <- x - mu
    d2 <- d^2
    ss <- sum(d2)
    ss
}
Browse[2]> n
debug at #2: d <- x - mu
Browse[2]> n
debug at #3: d2 <- d^2
Browse[2]> n
debug at #4: ss <- sum(d2)
Browse[2]> n
debug at #5: ss
Browse[2]> n
exiting from: SS(1, xPoints)
[1] 202.482

This isn’t a particularly interesting use of the debugger, as all we did was to step through the function a line at a time, using the ‘n’ command, but it does illustrate some of the basic ideas:

First, when we enter the SS function, the function itself is printed, and execution stops.
The R prompt, ‘>’, changes to ‘Browse[]>’

There are four basic commands you can use in the debugger:

n: execute the current line of code and step to the next
c: continue execution without stopping again
Q: quit debugging completely and to back to the top-level R prompt
where: similar to calling traceback

You can execute any R code at the “Browse” prompt, i.e., not just code related to the function. The most important example of this is probably the ls() command, which shows us what objects are defined in the current environment (i.e., within a given function, for instance):

    Browse[2]> ls()
    [1] "mu" "x" 

    Browse[2]> mu
    [1] 1

    Browse[2]> head(x)
    [1] -0.79928868 -0.73009689  1.43687692  0.30502316 -0.39728401
    [6]  0.08889039

But the facility is completely general. Here’s a trivial example of defining something unrelated to the current function:

    Browse[2]> zzz <- sqrt(64)
    Browse[2]> zzz
    [1] 8

Here’s another, somewhat more interesting session with SS:

    > SS(1, xPoints)
    debugging in: SS(1, xPoints)
    debug at #1: {
        d <- x - mu
        d2 <- d^2
        ss <- sum(d2)
        ss
    }
    Browse[2]> n
    debug at #2: d <- x - mu
    Browse[2]> n
    debug at #3: d2 <- d^2
    Browse[2]> d[1]  ## Print the value of the first element of d
    [1] -1.799289
    Browse[2]> n
    debug at #4: ss <- sum(d2)
    Browse[2]> hist(d2)
    Browse[2]> n
    debug at #5: ss
    Browse[2]> print(ss)
    [1] 202.482
    Browse[2]> ls()
    [1] "d"  "d2" "mu" "ss" "x" 
    Browse[2]> where
    where 1: SS(1, xPoints)
    
    Browse[2]> y <- x^2  ## Create new object
    Browse[2]> ls()
    [1] "d"  "d2" "mu" "ss" "x"  "y" 
    Browse[2]> head(y)
    [1] 0.638862388 0.533041463 2.064615286 0.093039131 0.157834589
    [6] 0.007901501
    Browse[2]> c  ## Execute the rest of the function w/o stopping
    exiting from: SS(1, xPoints)
    [1] 202.482

    > undebug(SS)  ## No need for further debugging

Note the explicit use of the print function in the above. At an interactive prompt, R will “autoprint” the value of a variable when it is used, so use of the print function is optional. But what if, for instance, you had n as a parameter in your function? Then to examine n you’d need to say explicitly:

print(n)

One final thing: if you want to debug a function, you have to say so before you call the function. This is more or less obvious.

Perhaps not so obvious is the fact that you can flag functions for debugging “on the fly”. I.e., you can turn on debugging for other functions while you’re in the process of debugging a given function. Example later.

Explicit Calls to Browser

Rather than flagging a function for debugging and then stepping through the function a line at a time, it is sometimes convenient to modify your code in such a way as to trigger debugging at the point where you suspect the trouble lies.

Consider the following modified version of SS, for instance. We might feel confident that the trouble occurs in our function only after d2 has been calculated, so we insert a call to the function browser at the point where we want to start, well, browsing.


SS <- function(mu, x) {
    d <- x - mu
    d2 <- d^2
    browser()
    ss <- sum(d2)
    ss
}

    > SS(2, xPoints)
    Called from: SS(2, xPoints)
    Browse[1]> ls()
    [1] "d"  "d2" "mu" "x" 
    Browse[1]> print(mu)
    [1] 2
    Browse[1]> mean(x)
    [1] -0.05553729
    Browse[1]> n
    debug at #5: ss <- sum(d2)
    Browse[2]> c
    [1] 513.5895

One drawback to this mode of operation is that you have to remember to remove the call to browser before you put the code into production.

One real advantage of this approach is that you can use it to debug “anonymous” functions, as:

    sapply(1:50, function_with_no_name {
        .
        .
        .
        browser()
        .
        .
        .
    })

Inserting Code with `trace`

As with debug, the function trace allows you to debug code, but trace is much more general and sophisticated. The simplest uses of trace don’t appear much different from the corresponding use of debug, although some differences are evident.

The following example is taken from the book Software for Data Analysis by John M. Chambers. Suppose we have defined three functions, called Fn1, Fn2, Fn3. (The details of the function don’t matter for the purposes of this example.) Then we might use trace as follows:

    trace(Fn1, recover)
    trace(Fn2, exit = browser)
    trace(Fn3, browser, exit = browser)

Here, browser is the same function we saw above. The function recover lets us examine the full sequence of calls that brought us to a given point in our code, and then browse inside any of the functions in the sequence.

After the above calls to trace, all future calls to Fn1 will begin with a call to recover(); calls to Fn2 will invoke browser() just before the function exits; calls to Fn3 will invoke browser() on both entry to the function and exit from the function. (In the case of Fn3 we might get tired of stepping through the function and hit c (continue), but stop before we leave the function, juuust to make sure everything is still OK before we lose all the local variables in the function.)

The trace function works by modifying the code of the function, then saving the modified function, and running the modified version when you invoke the function. The original code is not touched, and you can restore the original behavior by executing:

    untrace(Fn1)
    etc.

Here are some examples with our SS function:

trace(SS, recover)
SS(1, xPoints)
untrace(SS)

trace(SS, browser)
SS(2, xPoints)
untrace(SS)

trace(SS, exit = browser)
SS(3, xPoints)
untrace(SS)

The paper by Roger Peng, mentioned at the beginning of this note, has an example in which he minimizes the negative log-likelihood of a statistical model, using an R function that he calls nLL. The details of the discussion are somewhat complicated and will be omitted here.

The salient point is that the procedure in question involves many iterations, and is plagued by the same problem that we discussed above: the log of a negative number produces a so-called NaN (“not a number”) in some of the iterations. We would hope not to get NaN’s in our calculations, and we might, therefore, want to examine in detail the situations in which a NaN is produced. But there might be tens or hundreds or thousands of iterations, most of which do not generate NaN’s. Stepping through all the iterations would take a lifetime. What to do?

The solution is to trace the function using some of its more sophisticated options:

trace("nLL", quote(if(any(is.nan(lz))) { browser() }), at=4, print=F)

Here:

nLL is the name of the function to be traced
the stuff inside the quote function is code to be inserted (dynamically) into the nLL function; it could also be just the name of a function (as in our introductory examples)
at is the line number at which to insert the code, which the author determines from as.list(body(nLL)) (see above)
print=FALSE just suppresses the printing of a descriptive message

This is as if we had written:

    if(any(is.nan(lz))) {
        browser()
    }

at line 4 of the nLL function.

Here’s a simpler example using the functions we mentioned at the beginning of the discussion:


trace(h, quote(if (is.nan(r)) {
    browser()
}), at = 3)
f(7)
f(-3)
untrace(f)

Finally, if using the quoted expression for the second argument of trace seems too intimidating. let’s note that trace has an edit argument. If you call trace as:

    trace (someFn, ..., edit=TRUE)

you will dropped into an editor, where you can edit your code to your heart’s content before the (temporarily modified) function executes. BTW, you can also supply the name of an editor here, as:

    trace (someFn, ..., edit="/usr/bin/vim")

Editor’s Addendum: A quick way to enable bugging everywhere

One way to get in the habit of debugging is to drop this into your .Rprofile file or run it at the start of your R session:

options(error=utils::recover)

This will automatically invoke recover() every time your code gives an error, allowing you to view the stack of function calls and enter the browser mode in any of the functions. For instance:

> f("a")
Error in log(z) : non-numeric argument to mathematical function

Enter a frame number, or 0 to exit   

1: f("a")
2: #2: g(x)
3: #2: h(y)

Now selecting any of these will start browser mode within those functions. If you don’t need to use the browser (for instance, if the source of the error is obvious), just enter 0 and move on.

Demographic analysis using the `popbio` library and some other fun stuff

2013-04-04T11:16:55-07:00

This week at the Davis R Users’ Group we had a great presentation by Kara Moore O’Leary on using the popbio package to examine rare plant demographics. The following is her script run through knitr. You can download the original script and associated data here. Find out more about Kara and her work at her website here

Demographic analysis using the popbio library and some other fun stuff

A population viability type analysis for a rare herbacious perennial plant, Penstemon albomarginatus, for its only remaining California population

Objectives:

Estimating the population growth rate - deterministic and stochastic methods
Sensitivity and elasticity analysis - which transitions/rate are most variable, sensitive to change
Projecting the population stochastically for different scenarios, observing variation in growth rate
Quasi extinction probability
Estimating vital rates from real data, which will never be as good/much the data that you would like to have for a PVA - some therapy for me, and maybe you too

Load popbio and some other helpful libraries. I use the dropdown menu in RStudio to install libraries, it seems to work more consistently for my mac than scripting the installs, sorry.
Am I missing something basic about installing packages?

library(popbio)
library(plyr)
library(reshape)

I have already done a bit of magic on my data get here, most or all of which is suggested in the Morris and Doak 2005 book Quantitative population biology. That book was an excellent guide to this kind of analysis. If you want to go deep (and deeper) into this kind of analysis you should get borrow or steal Morris and Doak 2005 and Caswell 2002 from the start! They are both really excellent and straightforward and you’ll be able to answer a lot of the questions that you have along the way. Most of the functions in popbio are based on functions in these two books. Many of the popbio examples are on their data and generate the same graphs etc.

What you need for a class/stage structured demographic model:

a bunch of individuals (maybe in different populations)
annual survival rates
annual class or stage transition rates
annual fecundity, ie. probability of contribution to the juvenile class in the next year
these data for a lot of years, preferably 10 or more (but we do what we can with less)

The crux for plants is that its challenging to 1) count all seeds produced annually and 2) know how many seeds really yeilds a juvenile in any year. I am doing a bit of work on these issues from limited field data. If you have a well behaved penguin that you can collar and track and you know produces 2 live juveniles each season, your PVA might be a bit simpler.

Our field data had surivival for each year, mean plant diameter, and inflorescence count classes. From this I found the median inflorescence number for each class MEDIAN_INLF I assigned each plant to CLASS based on its xDIAM_cm. I convinced myself that this had biological meaning by looking at the relationship between size class and survival. Whatever classes or stages you use, you should be confident that they are meaningful for your study species. OR use an integral projection model (IPM) instead. These allow you to use continuous variables like size or age rather than classes. nice abbreviated dataset

andre <- read.csv("karadat.csv")
str(andre)

## 'data.frame':    1221 obs. of  4 variables:
##  $ PLANT_UNQ: Factor w/ 395 levels "1_703","1_704",..: 238 240 249 253 258 263 250 255 259 239 ...
##  $ YEAR     : int  1994 1994 1994 1994 1994 1994 1994 1994 1994 1994 ...
##  $ CLASS    : Factor w/ 6 levels "A1","A2","A3",..: 1 1 1 1 1 1 2 2 2 3 ...
##  $ SEEDS    : num  71.8 0 71.8 0 71.8 ...

head2(andre)

##      PLANT_UNQ YEAR CLASS SEEDS
## 1          A_1 1994    A1 71.75
## 2          A_3 1994    A1  0.00
## 3         B_13 1994    A1 71.75
## .            .    .     .     .
## 1221    6_772A 2012     J  0.00

Let’s look at the stages/classes

levels(andre$CLASS)

## [1] "A1"   "A2"   "A3"   "A4"   "dead" "J"

# dead. After that it can be omitted.

1. ESTIMATING FECUNDITY

The big challenge for plants = estimating seeds/indiv –> juveniles produced the next year based on a few fruit and seed counts and a lot of inflorescence class data

Part 1: How to estimate seed production per plant?

Multiply seeds/infl from field observations by inflorescence number. Ideally you would have seed counts for each plant, but in the absence of those nearly impossible data, I’m using a few seeds/fruit counts * fruits/infl counts from 2011 and 2012.

Part 2: How many seeds makes an emergent juvenile?

I’ll add some notes on how I made these estimates at the end of this script, but let’s start today with a dataset ready to go for popbio I’m sure that every PVA has a hurdle or five to leap before you have the vital rate data ready to go. I’ve added estimated seeds produced per individual per year as a column in this simple dataset more on how I generated that later if we are interested

str(andre)

## 'data.frame':    1221 obs. of  4 variables:
##  $ PLANT_UNQ: Factor w/ 395 levels "1_703","1_704",..: 238 240 249 253 258 263 250 255 259 239 ...
##  $ YEAR     : int  1994 1994 1994 1994 1994 1994 1994 1994 1994 1994 ...
##  $ CLASS    : Factor w/ 6 levels "A1","A2","A3",..: 1 1 1 1 1 1 2 2 2 3 ...
##  $ SEEDS    : num  71.8 0 71.8 0 71.8 ...

2. PICK A STARTING POPULATION VECTOR

this is the # of individuals in each class/stage at the start of your model. you should play around with this to see how it effects the outcome. My model is insensitive to realistic changes in this vector. Here are the “options”, the number in each class observed annually

n_options <- ddply(andre, c("YEAR"), function(df) return(table(df$CLASS)))
n_options

##    YEAR A1 A2 A3 A4 dead   J
## 1  1994  6  3  2  2    0  13
## 2  1995 31 17 13 11    5  81
## 3  1996 49 15 14  6   52  16
## 4  1997 36 19 10  4   30   2
## 5  1998 10 22 13 17    9   0
## 6  1999  8 16 24 11    3   0
## 7  2000  6 17 17  8   11   0
## 8  2001  1 11 16 14    6   0
## 9  2002 13 14  6  0    9   0
## 10 2003  4  8  6  0   14   1
## 11 2004  0  0  0  0   19   0
## 12 2011 22 23 12  6    0 173
## 13 2012  6  1  0  1  172  34

picked starting population vector from 1995, the first year with 9 observed populations

n95 <- c(81, 31, 17, 13, 11)
n = n95

you can also see that I only have J class individuals for a subset of the years. Check your data for missing bits like this. Perhaps your collaborator is forgetting to tell you that they didn’t survey for new juveniles from 1996-2003. It happens!

3. GENERATE A TRANSITION MATRIX

this matrix links each individual to its fate in the next year/cycle/season this is why you need each individual to be “dead” for a year, but no longer. If your raw data is like mine and only taken on live plants, “dead” might be something you have to add. make columns for year2, fate, and seeds2 for the whole census

trans <- subset(merge(andre, andre, by = "PLANT_UNQ", sort = FALSE), YEAR.x == 
    YEAR.y - 1)
head2(trans)

##      PLANT_UNQ YEAR.x CLASS.x SEEDS.x YEAR.y CLASS.y SEEDS.y
## 3          A_1   1994      A1   71.75   1995      A2   258.3
## 8          A_1   1995      A2  258.30   1996    dead     0.0
## 13         A_3   1996      A1   71.75   1997      A2   258.3
## .            .      .       .       .      .       .       .
## 6117     9_317   2011       J   71.75   2012       J     0.0

rename rows and columns to improve clarity (I use the names used by popbio, which are similar to Morris and Doak) I think I worked out somewhat painfully that you need these particular column names.

rownames(trans) <- 1:nrow(trans)
colnames(trans)[1:7] <- c("plant", "year", "stage", "seeds", "year2", "fate", 
    "seeds2")
head2(trans)

##     plant year stage  seeds year2 fate seeds2
## 1     A_1 1994    A1  71.75  1995   A2  258.3
## 2     A_1 1995    A2 258.30  1996 dead    0.0
## 3     A_3 1996    A1  71.75  1997   A2  258.3
## .       .    .     .      .     .    .      .
## 825 9_317 2011     J  71.75  2012    J    0.0

add individual fertility estimates from the calculations above

seedlingtrans <- 0.00305  # This is the rate at which a seed becomes a J individual (I estimated this elsewhere, see Appendix below)

adding in the number of J individuals produced by each individual

trans$J <- trans$seeds * seedlingtrans  # note that J is not an integer, which is totally fine, its a rate of J production
head2(trans)

##     plant year stage  seeds year2 fate seeds2      J
## 1     A_1 1994    A1  71.75  1995   A2  258.3 0.2188
## 2     A_1 1995    A2 258.30  1996 dead    0.0 0.7878
## 3     A_3 1996    A1  71.75  1997   A2  258.3 0.2188
## .       .    .     .      .     .    .      .      .
## 825 9_317 2011     J  71.75  2012    J    0.0 0.2188

4. GENERATE ANNUAL MATRICES — THE SIMPLE WAY for 3 easy years

NAME STAGES

stages <- c("J", "A1", "A2", "A3", "A4")

you must have a vector named stages in this way for your own clases or stages used in the “stage” vector of your data frame

SET ITERATIONS

it <- 100  # set the number of time steps for a deterministic model

Make a demographic projection matrix for each year like so: 1994

trans94 <- subset(trans, year == 1994, c(plant, stage, fate, J))
(proj94 <- projection.matrix(trans94, stage, fate, J, sort = stages))  #this gives you a projection matrix for 1994

##     
##            J     A1     A2     A3    A4
##   J  0.07692 0.1459 0.7878 0.7878 4.355
##   A1 0.53846 0.3333 0.3333 0.0000 0.000
##   A2 0.07692 0.3333 0.0000 0.0000 0.000
##   A3 0.00000 0.1667 0.3333 0.5000 0.000
##   A4 0.00000 0.0000 0.3333 0.5000 1.000

you can do a simple deterministic projection of the matrix for just this year

(p94 <- pop.projection(proj94, n, it))

## $lambda
## [1] 1.496
## 
## $stable.stage
##       J      A1      A2      A3      A4 
## 0.48307 0.24671 0.07983 0.06803 0.12235 
## 
## $stage.vectors
##     0     1      2      3      4      5     6      7      8      9     10
## J  81 82.29 142.62 223.46 338.78 502.17 745.4 1112.5 1665.3 2493.0 3729.6
## A1 31 59.62  69.70 108.76 167.98 256.23 383.2  570.4  850.9 1272.2 1903.7
## A2 17 16.56  26.20  34.21  53.44  82.05 124.0  185.1  275.7  411.7  615.8
## A3 13 17.33  24.12  32.41  45.74  68.68 104.4  157.4  235.5  351.4  525.0
## A4 11 23.17  37.35  58.15  85.76 126.44 188.1  281.7  422.1  631.7  944.7
##        11   12    13    14    15    16    17    18     19     20     21
## J  5577.3 8340 12472 18653 27896 41720 62394 93313 139553 208707 312129
## A1 2848.1 4260  6370  9526 14247 21307 31865 47656  71271 106588 159407
## A2  921.5 1378  2061  3083  4610  6895 10311 15421  23063  34492  51584
## A3  785.1 1174  1757  2627  3929  5876  8787 13142  19654  29393  43959
## A4 1412.4 2112  3159  4724  7065 10567 15803 23633  35345  52859  79053
##        22     23      24      25      26      27      28      29       30
## J  466801 698120 1044066 1561441 2335197 3492380 5222992 7811190 11681943
## A1 238400 356536  533214  797443 1192607 1783591 2667431 3989249  5966079
## A2  77146 115374  172547  258051  385925  577166  863175 1290912  1930610
## A3  65742  98319  147041  219905  328877  491848  735578 1100086  1645222
## A4 118227 176813  264431  395467  591436  884516 1322829 1978343  2958690
##          31       32       33       34       35        36        37
## J  17470807 26128282 39075876 58439515 87398601 130708056 195479056
## A1  8922507 13343961 19956421 29845615 44635295  66753845  99833010
## A2  2887304  4318077  6457855  9657977 14443886  21601401  32305748
## A3  2460494  3679766  5503236  8230306 12308748  18408218  27530217
## A4  4424838  6617520  9896762 14800999 22135477  33104480  49509056
##           38        39        40        41        42        43        44
## J  292346643 437215941 653873694 977893912 1.462e+09 2.187e+09 3.271e+09
## A1 149304206 223290332 333939503 499419705 7.469e+08 1.117e+09 1.671e+09
## A2  48314520  72256272 108062106 161611144 2.417e+08 3.615e+08 5.406e+08
## A3  41172526  61575137  92088048 137721310 2.060e+08 3.080e+08 4.607e+08
## A4  74042747 110733850 165606843 247671569 3.704e+08 5.540e+08 8.285e+08
##           45        46        47        48        49        50        51
## J  4.892e+09 7.316e+09 1.094e+10 1.636e+10 2.447e+10 3.660e+10 5.474e+10
## A1 2.498e+09 3.736e+09 5.588e+09 8.357e+09 1.250e+10 1.869e+10 2.795e+10
## A2 8.085e+08 1.209e+09 1.808e+09 2.704e+09 4.044e+09 6.049e+09 9.046e+09
## A3 6.890e+08 1.030e+09 1.541e+09 2.305e+09 3.447e+09 5.154e+09 7.709e+09
## A4 1.239e+09 1.853e+09 2.771e+09 4.144e+09 6.198e+09 9.270e+09 1.386e+10
##           52        53        54        55        56        57        58
## J  8.186e+10 1.224e+11 1.831e+11 2.738e+11 4.095e+11 6.124e+11 9.159e+11
## A1 4.181e+10 6.252e+10 9.351e+10 1.398e+11 2.091e+11 3.128e+11 4.678e+11
## A2 1.353e+10 2.023e+10 3.026e+10 4.525e+10 6.768e+10 1.012e+11 1.514e+11
## A3 1.153e+10 1.724e+10 2.579e+10 3.856e+10 5.767e+10 8.625e+10 1.290e+11
## A4 2.073e+10 3.101e+10 4.637e+10 6.935e+10 1.037e+11 1.551e+11 2.320e+11
##           59        60        61        62        63        64        65
## J  1.370e+12 2.049e+12 3.064e+12 4.582e+12 6.852e+12 1.025e+13 1.533e+13
## A1 6.996e+11 1.046e+12 1.565e+12 2.340e+12 3.500e+12 5.234e+12 7.827e+12
## A2 2.264e+11 3.386e+11 5.063e+11 7.572e+11 1.132e+12 1.694e+12 2.533e+12
## A3 1.929e+11 2.885e+11 4.315e+11 6.453e+11 9.651e+11 1.443e+12 2.158e+12
## A4 3.469e+11 5.188e+11 7.759e+11 1.160e+12 1.736e+12 2.596e+12 3.882e+12
##           66        67        68        69        70        71        72
## J  2.292e+13 3.428e+13 5.127e+13 7.667e+13 1.147e+14 1.715e+14 2.565e+14
## A1 1.171e+13 1.751e+13 2.618e+13 3.916e+13 5.856e+13 8.758e+13 1.310e+14
## A2 3.788e+12 5.665e+12 8.473e+12 1.267e+13 1.895e+13 2.834e+13 4.238e+13
## A3 3.228e+12 4.828e+12 7.220e+12 1.080e+13 1.615e+13 2.415e+13 3.612e+13
## A4 5.805e+12 8.682e+12 1.298e+13 1.942e+13 2.904e+13 4.343e+13 6.495e+13
##           73        74        75        76        77        78        79
## J  3.836e+14 5.736e+14 8.579e+14 1.283e+15 1.919e+15 2.870e+15 4.292e+15
## A1 1.959e+14 2.930e+14 4.381e+14 6.552e+14 9.799e+14 1.466e+15 2.192e+15
## A2 6.339e+13 9.480e+13 1.418e+14 2.120e+14 3.171e+14 4.742e+14 7.092e+14
## A3 5.402e+13 8.079e+13 1.208e+14 1.807e+14 2.702e+14 4.041e+14 6.044e+14
## A4 9.714e+13 1.453e+14 2.173e+14 3.249e+14 4.860e+14 7.268e+14 1.087e+15
##           80        81        82        83        84        85        86
## J  6.418e+15 9.599e+15 1.436e+16 2.147e+16 3.211e+16 4.802e+16 7.181e+16
## A1 3.278e+15 4.902e+15 7.331e+15 1.096e+16 1.640e+16 2.452e+16 3.668e+16
## A2 1.061e+15 1.586e+15 2.372e+15 3.548e+15 5.306e+15 7.936e+15 1.187e+16
## A3 9.039e+14 1.352e+15 2.022e+15 3.024e+15 4.522e+15 6.763e+15 1.011e+16
## A4 1.626e+15 2.431e+15 3.636e+15 5.437e+15 8.132e+15 1.216e+16 1.819e+16
##           87        88        89        90        91        92        93
## J  1.074e+17 1.606e+17 2.402e+17 3.592e+17 5.373e+17 8.035e+17 1.202e+18
## A1 5.485e+16 8.203e+16 1.227e+17 1.835e+17 2.744e+17 4.104e+17 6.137e+17
## A2 1.775e+16 2.654e+16 3.970e+16 5.937e+16 8.879e+16 1.328e+17 1.986e+17
## A3 1.513e+16 2.262e+16 3.383e+16 5.059e+16 7.566e+16 1.132e+17 1.692e+17
## A4 2.720e+16 4.068e+16 6.084e+16 9.099e+16 1.361e+17 2.035e+17 3.043e+17
##           94        95        96        97        98        99
## J  1.797e+18 2.688e+18 4.020e+18 6.011e+18 8.990e+18 1.345e+19
## A1 9.178e+17 1.373e+18 2.053e+18 3.070e+18 4.591e+18 6.867e+18
## A2 2.970e+17 4.442e+17 6.643e+17 9.935e+17 1.486e+18 2.222e+18
## A3 2.531e+17 3.785e+17 5.661e+17 8.466e+17 1.266e+18 1.894e+18
## A4 4.552e+17 6.807e+17 1.018e+18 1.522e+18 2.277e+18 3.405e+18
## 
## $pop.sizes
##   [1] 1.530e+02 1.990e+02 3.000e+02 4.570e+02 6.917e+02 1.036e+03 1.545e+03
##   [8] 2.307e+03 3.449e+03 5.160e+03 7.719e+03 1.154e+04 1.726e+04 2.582e+04
##  [15] 3.861e+04 5.775e+04 8.636e+04 1.292e+05 1.932e+05 2.889e+05 4.320e+05
##  [22] 6.461e+05 9.663e+05 1.445e+06 2.161e+06 3.232e+06 4.834e+06 7.230e+06
##  [29] 1.081e+07 1.617e+07 2.418e+07 3.617e+07 5.409e+07 8.089e+07 1.210e+08
##  [36] 1.809e+08 2.706e+08 4.047e+08 6.052e+08 9.051e+08 1.354e+09 2.024e+09
##  [43] 3.027e+09 4.528e+09 6.771e+09 1.013e+10 1.514e+10 2.265e+10 3.387e+10
##  [50] 5.066e+10 7.576e+10 1.133e+11 1.695e+11 2.534e+11 3.790e+11 5.668e+11
##  [57] 8.477e+11 1.268e+12 1.896e+12 2.836e+12 4.241e+12 6.342e+12 9.485e+12
##  [64] 1.419e+13 2.121e+13 3.173e+13 4.745e+13 7.096e+13 1.061e+14 1.587e+14
##  [71] 2.374e+14 3.550e+14 5.309e+14 7.940e+14 1.187e+15 1.776e+15 2.656e+15
##  [78] 3.972e+15 5.940e+15 8.884e+15 1.329e+16 1.987e+16 2.972e+16 4.444e+16
##  [85] 6.646e+16 9.940e+16 1.487e+17 2.223e+17 3.325e+17 4.973e+17 7.437e+17
##  [92] 1.112e+18 1.663e+18 2.488e+18 3.720e+18 5.564e+18 8.321e+18 1.244e+19
##  [99] 1.861e+19 2.783e+19
## 
## $pop.changes
##  [1] 1.300 1.508 1.523 1.514 1.497 1.492 1.493 1.495 1.496 1.496 1.496
## [12] 1.495 1.495 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496
## [23] 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496
## [34] 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496
## [45] 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496
## [56] 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496
## [67] 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496
## [78] 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496
## [89] 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496 1.496

(l94 <- p94$lambda)  # wow! if we looked only at 1994 based on these estimates the population would be booming!

## [1] 1.496

stable.stage shows the proportion of the population in each stage class at the mythical equilibrium, 48% of plants are juveniles in 100 years note that the population has gone from around 200 individuals in 1994 to 2.8e+19 in 100 years, nice!

Now make some matrices for other years 1995

trans95 <- subset(trans, year == 1995, c(plant, stage, fate, J))
(proj95 <- projection.matrix(trans95, stage, fate, J, sort = stages))

##     
##           J     A1     A2      A3      A4
##   J  0.2000 0.2626 1.0298 3.52160 6.21499
##   A1 0.2125 0.7419 0.4706 0.07692 0.00000
##   A2 0.0000 0.1290 0.4706 0.15385 0.09091
##   A3 0.0000 0.0000 0.0000 0.76923 0.36364
##   A4 0.0000 0.0000 0.0000 0.00000 0.54545

2011

trans11 <- subset(trans, year == 2011, c(plant, stage, fate, J))
(proj11 <- projection.matrix(trans11, stage, fate, J, sort = stages))

##     
##           J      A1     A2     A3    A4
##   J  0.1573 0.98329 4.5710 8.0648 7.922
##   A1 0.0000 0.05556 0.3077 0.1429 0.000
##   A2 0.0000 0.00000 0.0000 0.0000 0.500
##   A3 0.0000 0.00000 0.0000 0.0000 0.000
##   A4 0.0000 0.00000 0.0000 0.1429 0.000


p95 <- pop.projection(proj95, n, it)
(l95 <- p95$lambda)  # lambda is much lower in 1995

## [1] 0.9954


(p11 <- pop.projection(proj11, n, pi))

## $lambda
## [1] 0.3065
## 
## $stable.stage
##        J       A1       A2       A3       A4 
## 0.969161 0.021632 0.009207 0.000000 0.000000 
## 
## $stage.vectors
##     0       1       2
## J  81 312.916 97.7480
## A1 31   8.810  2.1818
## A2 17   5.500  0.9286
## A3 13   0.000  0.0000
## A4 11   1.857  0.0000
## 
## $pop.sizes
## [1] 153.0 329.1 100.9
## 
## $pop.changes
## [1] 2.1509 0.3065

(l11 <- p11$lambda)  # and based on 2011 alone extinction is eminent. The gist here is we need lots of years of data to make any decent estimation

## [1] 0.3065

# of what the population is really likely to do (ie more than the three
# here)

5. DETERMINISTIC ANALAYSIS OF THESE 3 WELL BEHAVED YEARS - basic non-stochastic PVA

thesearethemeanprojmats <- list(proj94, proj95, proj11)  # make a list of the three matrices
(meanxprojmat <- mean(thesearethemeanprojmats))  # make a mean of the three projection matrices for deterministic analysis

##     
##            J      A1     A2      A3     A4
##   J  0.14475 0.46393 2.1295 4.12473 6.1639
##   A1 0.25032 0.37694 0.3705 0.07326 0.0000
##   A2 0.02564 0.15412 0.1569 0.05128 0.1970
##   A3 0.00000 0.05556 0.1111 0.42308 0.1212
##   A4 0.00000 0.00000 0.1111 0.21429 0.5152


n  # recall what n is, our starting population vector, ie the # of individuals in each class at the start of the projection

## [1] 81 31 17 13 11

(pprojme <- pop.projection(meanxprojmat, n))  # do the deterministic projection, lambda is the dominant left eigenvalue

## $lambda
## [1] 1.089
## 
## $stable.stage
##       J      A1      A2      A3      A4 
## 0.61436 0.25469 0.06677 0.03730 0.02688 
## 
## $stage.vectors
##     0      1       2       3       4       5      6      7      8      9
## J  81 183.73 177.921 182.074 196.519 214.113 233.26 253.98 276.50 301.01
## A1 31  39.21  66.116  75.791  82.085  88.980  96.75 105.30 114.63 124.79
## A2 17  12.35  15.265  19.380  21.551  23.374  25.38  27.61  30.05  32.71
## A3 13  10.44   9.224  10.355  11.748  12.972  14.16  15.42  16.79  18.28
## A4 11  10.34   8.938   8.277   8.636   9.361  10.20  11.11  12.10  13.17
##        10     11     12     13     14     15     16     17     18     19
## J  327.68 356.72 388.33 422.74 460.21 500.99 545.38 593.72 646.33 703.61
## A1 135.85 147.88 160.99 175.26 190.79 207.69 226.10 246.14 267.95 291.69
## A2  35.61  38.77  42.20  45.94  50.02  54.45  59.27  64.53  70.24  76.47
## A3  19.90  21.66  23.58  25.67  27.94  30.42  33.11  36.05  39.24  42.72
## A4  14.33  15.60  16.99  18.49  20.13  21.92  23.86  25.97  28.27  30.78
## 
## $pop.sizes
##  [1]  153.0  256.1  277.5  295.9  320.5  348.8  379.7  413.4  450.1  490.0
## [11]  533.4  580.6  632.1  688.1  749.1  815.5  887.7  966.4 1052.0 1145.3
## 
## $pop.changes
##  [1] 1.674 1.083 1.066 1.083 1.088 1.089 1.089 1.089 1.089 1.089 1.089
## [12] 1.089 1.089 1.089 1.089 1.089 1.089 1.089 1.089

(DetLamb <- pprojme$lambda)

## [1] 1.089

6. RUN STOCHASTIC ANALYSES ON THESE 3 WELL BEHAVED YEARS

For a stochastic analysis, include all of the annual matrices, then make a random draw with replacement

For a series of time step, or until stable stage distribution is reached.

thesearethemeanprojmats  # our list of projection matrices

## [[1]]
##     
##            J     A1     A2     A3    A4
##   J  0.07692 0.1459 0.7878 0.7878 4.355
##   A1 0.53846 0.3333 0.3333 0.0000 0.000
##   A2 0.07692 0.3333 0.0000 0.0000 0.000
##   A3 0.00000 0.1667 0.3333 0.5000 0.000
##   A4 0.00000 0.0000 0.3333 0.5000 1.000
## 
## [[2]]
##     
##           J     A1     A2      A3      A4
##   J  0.2000 0.2626 1.0298 3.52160 6.21499
##   A1 0.2125 0.7419 0.4706 0.07692 0.00000
##   A2 0.0000 0.1290 0.4706 0.15385 0.09091
##   A3 0.0000 0.0000 0.0000 0.76923 0.36364
##   A4 0.0000 0.0000 0.0000 0.00000 0.54545
## 
## [[3]]
##     
##           J      A1     A2     A3    A4
##   J  0.1573 0.98329 4.5710 8.0648 7.922
##   A1 0.0000 0.05556 0.3077 0.1429 0.000
##   A2 0.0000 0.00000 0.0000 0.0000 0.500
##   A3 0.0000 0.00000 0.0000 0.0000 0.000
##   A4 0.0000 0.00000 0.0000 0.1429 0.000


stochme <- stoch.growth.rate(thesearethemeanprojmats, prob = NULL, maxt = 50000, 
    verbose = TRUE)

## [1] Calculating stochastic growth at time 1
## [1] Calculating stochastic growth at time 10000
## [1] Calculating stochastic growth at time 20000
## [1] Calculating stochastic growth at time 30000
## [1] Calculating stochastic growth at time 40000
## [1] Calculating stochastic growth at time 50000

note that these stochastic approximations of lambda are in log form (not immediately comparable to pop.project$lambda)

exp(stochme$approx)  # is the analytic approximation of lambda via Tuljapakar's method

## [1] 0.946

stochme$approx  # this is more accurate (perhaps) when there is a lot of covariation in matrix elements

## [1] -0.05556

exp(stochme$sim)  # gives stochastic growth rate by simulation, random draws of whole

## [1] 0.8884

stochme$sim

## [1] -0.1183

Tuljapurkar’s approximation takes into account how stochastic variation in the matrix elements affects the long-term stochastic growth rate (Caswell 2001). It can be more accurate in cases where there is covariation between matrix elements within the same year but may not be as accurate when there is a high level of temporal variation (Morris and Doak 2002, Stubben et al. 2012).

Fun with stochastic analyses!

Its easy to give years “weights” in the stochastic model. For example, you can increase the drought rate by weighting drought years (2011)

yearweight <- c(1, 1, 2)
moredrought <- stoch.projection(thesearethemeanprojmats, n, tmax = 50, prob = yearweight, 
    nreps = 500)

the output is population sizes, which are fun to graph when comparing models

yearweight <- c(1, 1, 0)
nodrought <- stoch.projection(thesearethemeanprojmats, n, tmax = 50, prob = yearweight, 
    nreps = 500)

par(mfrow = c(2, 1))
hist(log(apply(moredrought, 1, sum)), col = "blue", density = 50, ylim = c(0, 
    150), xlim = c(-1.3, 25), xlab = "", main = "More drought")
abline(v = log10(200), lty = 3)  # puts a line at the starting population size for reference
hist(log(apply(nodrought, 1, sum)), col = "green3", density = 50, ylim = c(0, 
    150), xlim = c(-1.3, 25), xlab = "", main = "No drought")
abline(v = log10(200), lty = 3)

you can get fancy and put these on the same graph too to compare outcomes. y axis is frequency of final population size at tmax.

7. QUASI EXTINCTION BASED ON THESE 3 WELL BEHAVED YEARS

Another useful way to think about populations. Since our ability to really estimate lambda is based on the assumption of equilibrium at stable stage, it might be more realistic to think about comparing extinction probabilities for different scenarios that comparing lambdas. these are based on stochastic runs


obsd <- stoch.quasi.ext(thesearethemeanprojmats, n, prob = c(1, 1, 1), Nx = 10, 
    tmax = 50, maxruns = 10, nreps = 500, sumweight = c(1, 1, 1), verbose = FALSE)

drt <- stoch.quasi.ext(thesearethemeanprojmats, n, prob = c(1, 1, 2), Nx = 10, 
    tmax = 50, maxruns = 10, nreps = 500, sumweight = c(1, 1, 1), verbose = FALSE)

par(mfrow = c(2, 1))
matplot(obsd, ylab = "Quasi-extinction probability", ylim = c(0, 1.1), type = "l", 
    lty = 1, col = rainbow(10), las = 1, main = "Observed climate", xlab = "Years")
matplot(drt, ylab = "Quasi-extinction probability", ylim = c(0, 1.1), type = "l", 
    lty = 1, col = rainbow(10), las = 1, main = "Double drought", xlab = "Years")

8. SENSITIVITY AND ELASTICITY

SENSITIVITY is a measure of the amount of change is lambda give a small change in a matrix element.

ELASTICITY is a measure of ``proportional’’ effect, i.e., the effect that a change in a given matrix element has as a proportional to the change in that element

meanxprojmat

##     
##            J      A1     A2      A3     A4
##   J  0.14475 0.46393 2.1295 4.12473 6.1639
##   A1 0.25032 0.37694 0.3705 0.07326 0.0000
##   A2 0.02564 0.15412 0.1569 0.05128 0.1970
##   A3 0.00000 0.05556 0.1111 0.42308 0.1212
##   A4 0.00000 0.00000 0.1111 0.21429 0.5152

for an overall look at sensitivity and elasticity use the mean projection matrix you could do separate analyses by year or type of year too to examine how sensitivity and elasticity vary among years

(eigout <- eigen.analysis(meanxprojmat))  # do the associated sensitivity analysis

## $lambda1
## [1] 1.089
## 
## $stable.stage
##       J      A1      A2      A3      A4 
## 0.61436 0.25469 0.06677 0.03730 0.02688 
## 
## $sensitivities
##     
##           J      A1      A2      A3       A4
##   J  0.2255 0.09351 0.02451 0.01369 0.009867
##   A1 0.6930 0.28729 0.07532 0.04208 0.000000
##   A2 1.5373 0.63734 0.16708 0.09334 0.067253
##   A3 0.0000 1.13131 0.29658 0.16569 0.119377
##   A4 0.0000 0.00000 0.38356 0.21428 0.154386
## 
## $elasticities
##     
##            J      A1      A2       A3      A4
##   J  0.02999 0.03985 0.04795 0.051889 0.05587
##   A1 0.15935 0.09948 0.02564 0.002832 0.00000
##   A2 0.03621 0.09023 0.02408 0.004397 0.01217
##   A3 0.00000 0.05773 0.03027 0.064394 0.01329
##   A4 0.00000 0.00000 0.03915 0.042180 0.07306
## 
## $repro.value
##      J     A1     A2     A3     A4 
##  1.000  3.072  6.816 12.099 15.647 
## 
## $damping.ratio
## [1] 4.215


colSums(eigout$sensitivities)  # this gives the cumulative sensitivity of each stage/class, fun to graph

##      J     A1     A2     A3     A4 
## 2.4559 2.1494 0.9471 0.5291 0.3509

calculate fertility and survival sums

(projsums <- colSums(meanxprojmat))

##      J     A1     A2     A3     A4 
## 0.4207 1.0505 2.8792 4.8866 6.9973

(fert_row <- meanxprojmat[1, ])

##      J     A1     A2     A3     A4 
## 0.1448 0.4639 2.1295 4.1247 6.1639

(surv_row <- projsums - fert_row)

##      J     A1     A2     A3     A4 
## 0.2760 0.5866 0.7496 0.7619 0.8333

can do this for sens and elas too to make a bar chart

make x where columes are name, sens and elas

for my data I’m summing elasticity and sensitivity for each class. You could also pick the vital rates that are most meaningful in your own analysis

(eigout <- eigen.analysis(meanxprojmat))

## $lambda1
## [1] 1.089
## 
## $stable.stage
##       J      A1      A2      A3      A4 
## 0.61436 0.25469 0.06677 0.03730 0.02688 
## 
## $sensitivities
##     
##           J      A1      A2      A3       A4
##   J  0.2255 0.09351 0.02451 0.01369 0.009867
##   A1 0.6930 0.28729 0.07532 0.04208 0.000000
##   A2 1.5373 0.63734 0.16708 0.09334 0.067253
##   A3 0.0000 1.13131 0.29658 0.16569 0.119377
##   A4 0.0000 0.00000 0.38356 0.21428 0.154386
## 
## $elasticities
##     
##            J      A1      A2       A3      A4
##   J  0.02999 0.03985 0.04795 0.051889 0.05587
##   A1 0.15935 0.09948 0.02564 0.002832 0.00000
##   A2 0.03621 0.09023 0.02408 0.004397 0.01217
##   A3 0.00000 0.05773 0.03027 0.064394 0.01329
##   A4 0.00000 0.00000 0.03915 0.042180 0.07306
## 
## $repro.value
##      J     A1     A2     A3     A4 
##  1.000  3.072  6.816 12.099 15.647 
## 
## $damping.ratio
## [1] 4.215

eigout$sensitivities

##     
##           J      A1      A2      A3       A4
##   J  0.2255 0.09351 0.02451 0.01369 0.009867
##   A1 0.6930 0.28729 0.07532 0.04208 0.000000
##   A2 1.5373 0.63734 0.16708 0.09334 0.067253
##   A3 0.0000 1.13131 0.29658 0.16569 0.119377
##   A4 0.0000 0.00000 0.38356 0.21428 0.154386


(fert_row_s <- eigout$sensitivities[1, ])

##        J       A1       A2       A3       A4 
## 0.225549 0.093506 0.024513 0.013695 0.009867

(surv_row_s <- eigout$sensitivities[2, ] + eigout$sensitivities[3, ] + eigout$sensitivities[4, 
    ] + eigout$sensitivities[5, ])

##      J     A1     A2     A3     A4 
## 2.2303 2.0559 0.9225 0.5154 0.3410

sensme <- t(cbind(fert_row_s, surv_row_s))


(fert_row_e <- eigout$elasticities[1, ])

##       J      A1      A2      A3      A4 
## 0.02999 0.03985 0.04795 0.05189 0.05587

(surv_row_e <- eigout$elasticities[2, ] + eigout$elasticities[3, ] + eigout$elasticities[4, 
    ] + eigout$elasticities[5, ])

##       J      A1      A2      A3      A4 
## 0.19556 0.24744 0.11913 0.11380 0.09852

(survtable <- t(rbind(surv_row_s, surv_row_e)))

##    surv_row_s surv_row_e
## J      2.2303    0.19556
## A1     2.0559    0.24744
## A2     0.9225    0.11913
## A3     0.5154    0.11380
## A4     0.3410    0.09852

(ferttable <- t(rbind(fert_row_s, fert_row_e)))

##    fert_row_s fert_row_e
## J    0.225549    0.02999
## A1   0.093506    0.03985
## A2   0.024513    0.04795
## A3   0.013695    0.05189
## A4   0.009867    0.05587


par(mfrow = c(2, 2))
mynames <- c("Sensitivity", "Elasticity")
barplot(t(survtable[, 1:2]), beside = TRUE, las = 1, ylim = c(0, 3.5), xlab = "Stage class", 
    main = "Growth and survival")
abline(h = 0)
barplot(t(ferttable[, 1:2]), beside = TRUE, las = 1, ylim = c(0, 0.25), xlab = "Stage class", 
    main = "Fertility")
abline(h = 0)
legend("topright", mynames, fill = grey.colors(2))

APPENDIX A. HOW I GENERATED A FECUNDITY ESTIMATE

I’m including this for plant folks who might like to see how I made fecundity estimates from real data I also have developed a script to simulate juvenile numbers and transition rates for the years in my dataset that are missing these data, and bootstraps of the whole model.

Load up the raw-ish data

andre <- read.csv("D__composite9_13_2012.csv")
str(andre)

## 'data.frame':    1693 obs. of  15 variables:
##  $ PLANT_UNQ   : Factor w/ 619 levels "1_701","1_702",..: 462 464 473 477 482 487 474 479 483 463 ...
##  $ YEAR        : int  1994 1994 1994 1994 1994 1994 1994 1994 1994 1994 ...
##  $ CLASS       : Factor w/ 6 levels "A1","A2","A3",..: 1 1 1 1 1 1 2 2 2 3 ...
##  $ OBS_ID      : int  1 3 20 24 11 16 21 26 12 2 ...
##  $ COHORT      : Factor w/ 18 levels "1","2","3","4",..: 10 10 11 11 11 11 11 11 11 10 ...
##  $ PLANT.NUM   : Factor w/ 460 levels "1","10","11",..: 1 23 5 9 41 168 6 11 42 12 ...
##  $ xDIAM_cm    : num  12 11 17 12 18 12 23 22 29 32 ...
##  $ SEEDS_DC    : num  26.7 26.7 26.7 26.7 26.7 ...
##  $ MEDIAN_INFL : int  5 0 5 0 5 5 18 18 18 18 ...
##  $ INFL_CLASS  : int  1 0 1 0 1 1 2 2 2 2 ...
##  $ MAX_INFL    : int  111 0 111 0 111 111 123 123 123 123 ...
##  $ OBS_INLF    : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ SEEDS_YR_INF: num  NA NA NA NA NA NA NA NA NA NA ...
##  $ CAGED       : Factor w/ 2 levels "N","Y": 1 1 1 1 1 1 1 1 1 1 ...
##  $ WATERED     : Factor w/ 2 levels "N","Y": 1 1 1 1 1 1 1 1 1 1 ...

BEST ESTIMATE OF SEED PRODUCTION: average of seeds/fruit, weighed average of fruits/infl where 0.14 is 1/8 drought years

andre$SEEDS <- andre$MEDIAN_INFL * 14.35  # = 14.35 seeds/infl

The big crux Part 2: How many seeds makes a J plant the next year? ie what is the transition rate or fecundity rate? observed juveniles in each year (sadly not all years have the same # of cohorts, so I adjust for that below)

seedling_yr <- ddply(andre, c("YEAR"), function(df) return(c(sdlgs = sum(df$CLASS == 
    "J"))))  # shows the OBSERVED # of J plants each year,

and some of these are not observation years and are omitted below

seedling_yr

##    YEAR sdlgs
## 1  1994    13
## 2  1995    81
## 3  1996    16
## 4  1997     2
## 5  1998     0
## 6  1999     0
## 7  2000     0
## 8  2001     0
## 9  2002     0
## 10 2003     1
## 11 2004     0
## 12 2011   362
## 13 2012    57

seedlings for 9 cohorts based on census from each year

sdl94 = seedling_yr$sdlgs[(seedling_yr$YEAR == "1994")] * 9/2  # this year had only 2 cohorts
sdl95 = seedling_yr$sdlgs[(seedling_yr$YEAR == "1995")]
sdl11 = seedling_yr$sdlgs[(seedling_yr$YEAR == "2011")]
sdl12 = 0  # in 2012 observed seedlings were 0

seedling_pick = c(sdl94, sdl95, sdl11, sdl12)
seedling_pick

## [1]  58.5  81.0 362.0   0.0

estimate seedlings as the mean of the 4 OBSERVED years: in the other years seedlings where not surveyed for

(seedlings = mean(seedling_pick))  # = 125.375 seedlings/year on average

## [1] 125.4

what’s the annual total seed production rate? add up all the seeds estimated to be produced in each year

seed_yr <- ddply(andre, c("YEAR"), function(df) return(c(sumseeds = sum(df$SEEDS))))
seed_yr

##    YEAR sumseeds
## 1  1994     4434
## 2  1995    45834
## 3  1996    13302
## 4  1997    15426
## 5  1998    51976
## 6  1999    45030
## 7  2000    18913
## 8  2001    45576
## 9  2002     2597
## 10 2003     6257
## 11 2004        0
## 12 2011   158065
## 13 2012     3071

adjust so that 1994 has an estimate of all cohorts based on the observed 2 cohorts

seed_yr$sumseeds[(seedling_yr$YEAR == "1994")] <- (seed_yr$sumseeds[(seedling_yr$YEAR == 
    "1994")] * 9/2)

get mean seeds/year

avg_seeds_p_yr <- mean(seed_yr$sumseeds)

for each of the 4 years in which seedlings were observed, calculate an estimate of the transition rate from seed –> J

str(sdl.trans94 <- sdl94/avg_seeds_p_yr)

##  num 0.00179

str(sdl.trans95 <- sdl95/avg_seeds_p_yr)

##  num 0.00247

str(sdl.trans11 <- sdl11/avg_seeds_p_yr)

##  num 0.011

sdl.trans12 <- 0  # no seedlings observed this year

looks at the options for transition rate

str(seedlingtrans_pick <- c(sdl.trans94, sdl.trans95, sdl.trans11, sdl.trans12))

##  num [1:4] 0.00179 0.00247 0.01105 0

pick the mean for this analysis

(seedlingtrans <- mean(seedlingtrans_pick))  # mean seedling transition rate  0.00305

## [1] 0.003826

Now remove watered and caged plants from main dataset, these have different survival and transition rates. I included them thus far because we needed to get seed production estimates for them. Since they might have contributed to the observed juveniles

andre <- andre[(andre$CAGED == "N"), ]
andre <- andre[(andre$WATERED == "N"), ]
head2(andre)

##      PLANT_UNQ YEAR CLASS OBS_ID COHORT PLANT.NUM xDIAM_cm SEEDS_DC
## 1          A_1 1994    A1      1      A         1     12.0     26.7
## 2          A_3 1994    A1      3      A         3     11.0     26.7
## 3         B_13 1994    A1     20      B        13     17.0     26.7
## .            .    .     .      .      .         .        .        .
## 1675    6_772A 2012     J    608      6      772A      4.7      0.0
##      MEDIAN_INFL INFL_CLASS MAX_INFL OBS_INLF SEEDS_YR_INF CAGED WATERED
## 1              5          1      111       NA           NA     N       N
## 2              0          0        0       NA           NA     N       N
## 3              5          1      111       NA           NA     N       N
## .              .          .        .        .            .     .       .
## 1675           0          0        0        0            0     N       N
##      SEEDS
## 1    71.75
## 2     0.00
## 3    71.75
## .        .
## 1675  0.00

reduce datafile to include only PLANT_UNQ, YEAR, CLASS and SEEDS

str(andre)

## 'data.frame':    1221 obs. of  16 variables:
##  $ PLANT_UNQ   : Factor w/ 619 levels "1_701","1_702",..: 462 464 473 477 482 487 474 479 483 463 ...
##  $ YEAR        : int  1994 1994 1994 1994 1994 1994 1994 1994 1994 1994 ...
##  $ CLASS       : Factor w/ 6 levels "A1","A2","A3",..: 1 1 1 1 1 1 2 2 2 3 ...
##  $ OBS_ID      : int  1 3 20 24 11 16 21 26 12 2 ...
##  $ COHORT      : Factor w/ 18 levels "1","2","3","4",..: 10 10 11 11 11 11 11 11 11 10 ...
##  $ PLANT.NUM   : Factor w/ 460 levels "1","10","11",..: 1 23 5 9 41 168 6 11 42 12 ...
##  $ xDIAM_cm    : num  12 11 17 12 18 12 23 22 29 32 ...
##  $ SEEDS_DC    : num  26.7 26.7 26.7 26.7 26.7 ...
##  $ MEDIAN_INFL : int  5 0 5 0 5 5 18 18 18 18 ...
##  $ INFL_CLASS  : int  1 0 1 0 1 1 2 2 2 2 ...
##  $ MAX_INFL    : int  111 0 111 0 111 111 123 123 123 123 ...
##  $ OBS_INLF    : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ SEEDS_YR_INF: num  NA NA NA NA NA NA NA NA NA NA ...
##  $ CAGED       : Factor w/ 2 levels "N","Y": 1 1 1 1 1 1 1 1 1 1 ...
##  $ WATERED     : Factor w/ 2 levels "N","Y": 1 1 1 1 1 1 1 1 1 1 ...
##  $ SEEDS       : num  71.8 0 71.8 0 71.8 ...

str(andre[, c(1:3, 16)])

## 'data.frame':    1221 obs. of  4 variables:
##  $ PLANT_UNQ: Factor w/ 619 levels "1_701","1_702",..: 462 464 473 477 482 487 474 479 483 463 ...
##  $ YEAR     : int  1994 1994 1994 1994 1994 1994 1994 1994 1994 1994 ...
##  $ CLASS    : Factor w/ 6 levels "A1","A2","A3",..: 1 1 1 1 1 1 2 2 2 3 ...
##  $ SEEDS    : num  71.8 0 71.8 0 71.8 ...

Go back to step 2.

Model Selection and Multi-Model Inference

2013-02-20T17:14:43-08:00

At D-RUG this week Rosemary Hartman presented a really useful case study in model selection, based on her work on frog habitat. Here is her code run through ‘knitr’. Original code and data are posted here.

(yes, I am just doing this for the flying monkey)

Editor’s note: we’re giving away flying monkey dolls from our sponsor, Revolution Analytics, to all our D-RUG presenters.

So, let’s say you want to find out where things are and why they are there. But there are a lot of reasons someone might be somewhere

Let’s say these things are fishermen, And we aren’t sure what the reasons are yet.

Start by making some good guesses. Maybe they all go to the big lakes. Maybe they all go to the lakes that CDFW stocks with fish. Maybe they go to the lakes with lots of frogs (not likely, but hey, why not try?) Once you have all your guesses, you need to figure out which one is right. There are two basic methods to use (or two that I have been exploring)

Method #1: The kitchen sink

load package “glmulti”

library(glmulti)

load the data

lakes.df2 <- read.csv("lakes.df2.csv")

create a model that has all the predictor variables you would like to test

global.model <-  glm(fish~ # number of fishermen
                       stocked+ # CDFW fish stocking y/n
                       cows+ # cattle grazing in area y/n
                       herp.rich+ # amphibian species richness
                       fish.rich+ # fish species richness
                       visitors+ # total number of visitors
                       logsize, # log of the area
                     data=lakes.df2) # data frame we got this from

Now we will use the “glmulti” function to find the best model; this goes through every possible model and finds the best one.

fish.model <- glmulti(global.model, # use the model with built as a starting point
                      level = 1,  #  just look at main effects
                      crit="aicc") # use AICc because it works better than AIC for small sample sizes

## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
## 
## After 50 models:
## Best model: fish~1+cows+herp.rich+visitors+logsize
## Crit= 201.820969266156
## Mean crit= 355.848076682524

## Completed.


summary(fish.model)

## $name
## [1] "glmulti.analysis"
## 
## $method
## [1] "h"
## 
## $fitting
## [1] "glm"
## 
## $crit
## [1] "aicc"
## 
## $level
## [1] 1
## 
## $marginality
## [1] FALSE
## 
## $confsetsize
## [1] 100
## 
## $bestic
## [1] 201.8
## 
## $icvalues
##  [1] 201.8 203.4 204.7 205.0 206.4 206.5 206.5 208.3 225.4 227.7 230.6
## [12] 231.2 232.7 233.2 233.5 235.1 236.0 236.4 238.2 238.3 238.4 238.8
## [23] 240.7 240.9 258.2 259.3 260.2 261.0 261.1 261.7 262.0 262.6 265.1
## [34] 267.6 268.1 270.5 276.0 277.4 278.3 279.8 310.2 311.5 313.6 314.0
## [45] 315.4 316.2 319.0 319.7 500.6 502.0 502.5 503.0 503.8 504.8 511.3
## [56] 513.2 599.7 601.4 610.5 614.9 615.3 624.1 629.2 638.6
## 
## $bestmodel
## [1] "fish ~ 1 + cows + herp.rich + visitors + logsize"
## 
## $modelweights
##  [1] 4.485e-01 2.024e-01 1.075e-01 9.285e-02 4.464e-02 4.343e-02 4.310e-02
##  [8] 1.767e-02 3.399e-06 1.065e-06 2.477e-07 1.895e-07 8.849e-08 6.754e-08
## [15] 6.011e-08 2.630e-08 1.704e-08 1.404e-08 5.716e-09 5.287e-09 4.990e-09
## [22] 4.134e-09 1.618e-09 1.436e-09 2.540e-13 1.465e-13 9.334e-14 6.306e-14
## [29] 5.969e-14 4.361e-14 3.790e-14 2.885e-14 8.119e-15 2.304e-15 1.837e-15
## [36] 5.540e-16 3.542e-17 1.739e-17 1.096e-17 5.151e-18 1.300e-24 7.003e-25
## [43] 2.397e-25 1.957e-25 9.964e-26 6.413e-26 1.594e-26 1.123e-26 5.852e-66
## [50] 2.976e-66 2.260e-66 1.801e-66 1.220e-66 7.203e-67 2.847e-68 1.069e-68
## [57] 1.836e-87 7.713e-88 8.222e-90 9.195e-91 7.394e-91 9.070e-93 7.028e-94
## [64] 6.396e-96
## 
## $includeobjects
## [1] TRUE

That showed us the best model, now lets look at some of the others

weightable(fish.model)

##                                                                     model
## 1                        fish ~ 1 + cows + herp.rich + visitors + logsize
## 2            fish ~ 1 + cows + herp.rich + fish.rich + visitors + logsize
## 3              fish ~ 1 + stocked + cows + herp.rich + visitors + logsize
## 4                                    fish ~ 1 + cows + visitors + logsize
## 5                        fish ~ 1 + cows + fish.rich + visitors + logsize
## 6                          fish ~ 1 + stocked + cows + visitors + logsize
## 7  fish ~ 1 + stocked + cows + herp.rich + fish.rich + visitors + logsize
## 8              fish ~ 1 + stocked + cows + fish.rich + visitors + logsize
## 9                         fish ~ 1 + stocked + cows + herp.rich + logsize
## 10            fish ~ 1 + stocked + cows + herp.rich + fish.rich + logsize
## 11                                  fish ~ 1 + cows + herp.rich + logsize
## 12                                    fish ~ 1 + stocked + cows + logsize
## 13                                              fish ~ 1 + cows + logsize
## 14                      fish ~ 1 + cows + herp.rich + fish.rich + logsize
## 15                        fish ~ 1 + stocked + cows + fish.rich + logsize
## 16                                  fish ~ 1 + cows + fish.rich + logsize
## 17                              fish ~ 1 + herp.rich + visitors + logsize
## 18                                          fish ~ 1 + visitors + logsize
## 19                  fish ~ 1 + herp.rich + fish.rich + visitors + logsize
## 20                              fish ~ 1 + fish.rich + visitors + logsize
## 21                    fish ~ 1 + stocked + herp.rich + visitors + logsize
## 22                                fish ~ 1 + stocked + visitors + logsize
## 23        fish ~ 1 + stocked + herp.rich + fish.rich + visitors + logsize
## 24                    fish ~ 1 + stocked + fish.rich + visitors + logsize
## 25                     fish ~ 1 + cows + herp.rich + fish.rich + visitors
## 26                                 fish ~ 1 + cows + herp.rich + visitors
## 27           fish ~ 1 + stocked + cows + herp.rich + fish.rich + visitors
## 28                                 fish ~ 1 + cows + fish.rich + visitors
## 29                       fish ~ 1 + stocked + cows + herp.rich + visitors
## 30                                             fish ~ 1 + cows + visitors
## 31                       fish ~ 1 + stocked + cows + fish.rich + visitors
## 32                                   fish ~ 1 + stocked + cows + visitors
## 33                               fish ~ 1 + stocked + herp.rich + logsize
## 34                   fish ~ 1 + stocked + herp.rich + fish.rich + logsize
## 35                                           fish ~ 1 + stocked + logsize
## 36                               fish ~ 1 + stocked + fish.rich + logsize
## 37                                                     fish ~ 1 + logsize
## 38                                         fish ~ 1 + herp.rich + logsize
## 39                                         fish ~ 1 + fish.rich + logsize
## 40                             fish ~ 1 + herp.rich + fish.rich + logsize
## 41                      fish ~ 1 + stocked + cows + herp.rich + fish.rich
## 42                                  fish ~ 1 + stocked + cows + fish.rich
## 43                                            fish ~ 1 + cows + fish.rich
## 44                                fish ~ 1 + cows + herp.rich + fish.rich
## 45                                  fish ~ 1 + stocked + cows + herp.rich
## 46                                              fish ~ 1 + stocked + cows
## 47                                                        fish ~ 1 + cows
## 48                                            fish ~ 1 + cows + herp.rich
## 49                            fish ~ 1 + herp.rich + fish.rich + visitors
## 50                  fish ~ 1 + stocked + herp.rich + fish.rich + visitors
## 51                                        fish ~ 1 + herp.rich + visitors
## 52                                        fish ~ 1 + fish.rich + visitors
## 53                              fish ~ 1 + stocked + herp.rich + visitors
## 54                              fish ~ 1 + stocked + fish.rich + visitors
## 55                                                    fish ~ 1 + visitors
## 56                                          fish ~ 1 + stocked + visitors
## 57                             fish ~ 1 + stocked + herp.rich + fish.rich
## 58                                         fish ~ 1 + stocked + fish.rich
## 59                                         fish ~ 1 + stocked + herp.rich
## 60                                                   fish ~ 1 + fish.rich
## 61                                       fish ~ 1 + herp.rich + fish.rich
## 62                                                     fish ~ 1 + stocked
## 63                                                   fish ~ 1 + herp.rich
## 64                                                               fish ~ 1
##     aicc   weights
## 1  201.8 4.485e-01
## 2  203.4 2.024e-01
## 3  204.7 1.075e-01
## 4  205.0 9.285e-02
## 5  206.4 4.464e-02
## 6  206.5 4.343e-02
## 7  206.5 4.310e-02
## 8  208.3 1.767e-02
## 9  225.4 3.399e-06
## 10 227.7 1.065e-06
## 11 230.6 2.477e-07
## 12 231.2 1.895e-07
## 13 232.7 8.849e-08
## 14 233.2 6.754e-08
## 15 233.5 6.011e-08
## 16 235.1 2.630e-08
## 17 236.0 1.704e-08
## 18 236.4 1.404e-08
## 19 238.2 5.716e-09
## 20 238.3 5.287e-09
## 21 238.4 4.990e-09
## 22 238.8 4.134e-09
## 23 240.7 1.618e-09
## 24 240.9 1.436e-09
## 25 258.2 2.540e-13
## 26 259.3 1.465e-13
## 27 260.2 9.334e-14
## 28 261.0 6.306e-14
## 29 261.1 5.969e-14
## 30 261.7 4.361e-14
## 31 262.0 3.790e-14
## 32 262.6 2.885e-14
## 33 265.1 8.119e-15
## 34 267.6 2.304e-15
## 35 268.1 1.837e-15
## 36 270.5 5.540e-16
## 37 276.0 3.542e-17
## 38 277.4 1.739e-17
## 39 278.3 1.096e-17
## 40 279.8 5.151e-18
## 41 310.2 1.300e-24
## 42 311.5 7.003e-25
## 43 313.6 2.397e-25
## 44 314.0 1.957e-25
## 45 315.4 9.964e-26
## 46 316.2 6.413e-26
## 47 319.0 1.594e-26
## 48 319.7 1.123e-26
## 49 500.6 5.852e-66
## 50 502.0 2.976e-66
## 51 502.5 2.260e-66
## 52 503.0 1.801e-66
## 53 503.8 1.220e-66
## 54 504.8 7.203e-67
## 55 511.3 2.847e-68
## 56 513.2 1.069e-68
## 57 599.7 1.836e-87
## 58 601.4 7.713e-88
## 59 610.5 8.222e-90
## 60 614.9 9.195e-91
## 61 615.3 7.394e-91
## 62 624.1 9.070e-93
## 63 629.2 7.028e-94
## 64 638.6 6.396e-96

So this is the best model

f <- glm(fish ~ 1 + cows + herp.rich + visitors + logsize, data = lakes.df2)
summary(f)

## 
## Call:
## glm(formula = fish ~ 1 + cows + herp.rich + visitors + logsize, 
##     data = lakes.df2)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -4.744  -1.350  -0.243   0.986  10.483  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.5262     0.8914   -0.59    0.559    
## cowsy        -1.5214     1.2088   -1.26    0.217    
## herp.rich     0.7614     0.3219    2.37    0.024 *  
## visitors      0.1799     0.0277    6.49  1.8e-07 ***
## logsize       1.0498     1.1590    0.91    0.371    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## (Dispersion parameter for gaussian family taken to be 7.225)
## 
##     Null deviance: 791.60  on 39  degrees of freedom
## Residual deviance: 252.87  on 35  degrees of freedom
##   (79 observations deleted due to missingness)
## AIC: 199.3
## 
## Number of Fisher Scoring iterations: 2

I hate trying to interpret models based on tables of coeficients. Let’s look at some graphs

library(visreg)
visreg(f)

But according to theory, models with AIC within two points of each other are basically equal. So what about the other models? Should we totally throw them out? Looking at the table of aicc weights, there is a pretty big jump between model 8 and model 9. So lets try averaging the top 8 models if we want to use this to make predictions, evaluate variable importance, etc.

library(MuMIn)

run the top 8 models

f2 <- glm(fish ~ 1 + cows + herp.rich + fish.rich + visitors + logsize, data = lakes.df2)
f3 <- glm(fish ~ 1 + stocked + cows + herp.rich + visitors + logsize, data = lakes.df2)
f4 <- glm(fish ~ 1 + cows + visitors + logsize, data = lakes.df2)
f5 <- glm(fish ~ 1 + cows + fish.rich + visitors + logsize, data = lakes.df2)
f6 <- glm(fish ~ 1 + stocked + cows + visitors + logsize, data = lakes.df2)
f7 <- glm(fish ~ 1 + stocked + cows + herp.rich + fish.rich + visitors + logsize, 
    data = lakes.df2)
f8 <- glm(fish ~ 1 + stocked + cows + fish.rich + visitors + logsize, data = lakes.df2)

average the models together

f.ave <- model.avg(f, f2, f3, f4, f5, f6, f7, f8)
summary(f.ave)

## 
## Call:
## model.avg.default(object = f, f2, f3, f4, f5, f6, f7, f8)
## 
## Component models:
##        df logLik  AICc Delta Weight
## 1346    6 -93.64 201.8  0.00   0.45
## 12346   7 -92.96 203.4  1.59   0.20
## 13456   7 -93.59 204.7  2.86   0.11
## 146     5 -96.60 205.0  3.15   0.09
## 1246    6 -95.95 206.4  4.61   0.04
## 1456    6 -95.97 206.5  4.67   0.04
## 123456  8 -92.93 206.5  4.68   0.04
## 12456   7 -95.39 208.3  6.47   0.02
## 
## Term codes:
##      cows fish.rich herp.rich   logsize   stocked  visitors 
##         1         2         3         4         5         6 
## 
## Model-averaged coefficients: 
##             Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)  -0.5422     1.2363      1.2673    0.43    0.669    
## cowsy        -1.3751     1.2596      1.3019    1.06    0.291    
## herp.rich     0.7523     0.3274      0.3393    2.22    0.027 *  
## visitors      0.1841     0.0305      0.0316    5.83   <2e-16 ***
## logsize       0.7763     1.2689      1.3119    0.59    0.554    
## fish.rich     0.7784     0.7245      0.7512    1.04    0.300    
## stockedy     -0.7316     1.6009      1.6539    0.44    0.658    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Full model-averaged coefficients (with shrinkage): 
##  (Intercept)  cowsy herp.rich visitors logsize fish.rich stockedy
##       -0.542 -1.375     0.603    0.184   0.776     0.240   -0.155
## 
## Relative variable importance:
##      cows   logsize  visitors herp.rich fish.rich   stocked 
##      1.00      1.00      1.00      0.80      0.31      0.21

Now we know the relative variable importance, all of the avereged coefficents, etc, and can use the “predict” funciton to predict how many fishermen will be at the next lake. Graphing these are more difficult, so we’ll skip this for now. Ask me later if you really want to know.

Model selection method #2: Use your brain We often can discard (or choose) some models a priori based on our knowlege of the system. It’s usually better to do it this way if you have several hundered possible combination of variables, or want to put in some interaction terms. I used this method for my frog data. load package bbmle

library(bbmle)

load the data

fishlakes <- read.csv("fishlakes.csv")

Decide on which set of models you want to use. This is hard. A statistician who knows a lot more than I do told me so. I spent a long time playing around with different transformations, predictor variables, and combinations before I came up with a set of hypotheses (models) that I was happy with.

m1 <- glm(treatment ~ 1 + logveg + bank.slope, family = binomial("logit"), data = fishlakes)
m2 <- glm(treatment ~ 1 + logveg + silt.total + bank.slope, family = binomial("logit"), 
    data = fishlakes)
m3 <- glm(treatment ~ 1 + raca.basin + logveg + bank.slope, family = binomial("logit"), 
    data = fishlakes)
m4 <- glm(treatment ~ 1 + BUBO.breeding + logveg + bank.slope, family = binomial("logit"), 
    data = fishlakes)
m5 <- glm(treatment ~ 1 + BUBO.breeding + logveg + bank.slope + raca.basin, 
    family = binomial("logit"), data = fishlakes)
m6 <- glm(treatment ~ 1 + raca.basin * BUBO.breeding, family = binomial("logit"), 
    data = fishlakes)
m7 <- glm(treatment ~ 1 + herbacious + raca.basin + lakes.basin, family = binomial("logit"), 
    data = fishlakes)
m8 <- glm(treatment ~ 1 + BUBO.breeding, family = binomial("logit"), data = fishlakes)
m9 <- glm(treatment ~ BUBO.breeding + herbacious + bank.slope + logveg, family = binomial("logit"), 
    data = fishlakes)
m10 <- glm(treatment ~ BUBO.breeding * raca.basin + herbacious + lakes.basin, 
    family = binomial("logit"), data = fishlakes)
m11 <- glm(treatment ~ 1 + bank.slope + logveg + BUBO.breeding:raca.basin, family = binomial("logit"), 
    data = fishlakes)

Now let’s rank them via AICc

AICctab(m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, base = T, weights = T, 
    nobs = length(fishlakes))

##     AICc df dAICc weight 
## m11 35.9 5   0.0  0.79445
## m4  40.5 4   4.6  0.08072
## m9  40.9 5   5.1  0.06336
## m1  42.9 3   7.0  0.02424
## m5  43.0 5   7.1  0.02275
## m2  45.3 4   9.4  0.00734
## m3  45.3 4   9.4  0.00708
## m8  56.6 2  20.7  < 0.001
## m10 57.1 6  21.2  < 0.001
## m6  57.3 4  21.4  < 0.001
## m7  62.1 4  26.2  < 0.001

Model m11 comes out ahead by quite a bit, but we’ll average the top two models, just to show you how its done.

m.ave <- model.avg(m4, m11)
summary(m.ave)

## 
## Call:
## model.avg.default(object = m4, m11)
## 
## Component models:
##     df logLik  AICc Delta Weight
## 134  5 -12.09 35.84  0.00   0.91
## 123  4 -15.68 40.43  4.59   0.09
## 
## Term codes:
##               bank.slope            BUBO.breeding                   logveg 
##                        1                        2                        3 
## BUBO.breeding:raca.basin 
##                        4 
## 
## Model-averaged coefficients: 
##                           Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept)                 -1.356      3.575       3.694    0.37    0.714
## bank.slope                 -19.935      8.233       8.503    2.34    0.019
## logveg                       1.783      1.366       1.409    1.26    0.206
## BUBO.breedingn:raca.basin   -0.965      0.973       1.006    0.96    0.337
## BUBO.breedingy:raca.basin    1.313      0.699       0.722    1.82    0.069
## BUBO.breedingy               2.077      1.055       1.090    1.91    0.057
##                            
## (Intercept)                
## bank.slope                *
## logveg                     
## BUBO.breedingn:raca.basin  
## BUBO.breedingy:raca.basin .
## BUBO.breedingy            .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Full model-averaged coefficients (with shrinkage): 
##  (Intercept) bank.slope  logveg BUBO.breedingn:raca.basin
##       -1.356    -19.935   1.783                    -0.877
##  BUBO.breedingy:raca.basin BUBO.breedingy
##                      1.193          0.190
## 
## Relative variable importance:
##               bank.slope                   logveg BUBO.breeding:raca.basin 
##                     1.00                     1.00                     0.91 
##            BUBO.breeding 
##                     0.09

# OK, there is a predictive model, but how good is it?

Let’s try cross-validation first. If this was a single model, we could try using the cv.glm function from the “boot” library like this:

library("boot")
cost <- function(r, pi = 0) mean(abs(r - pi) > 0.5)  ## cost function necessary for binomial data
m11.cv <- cv.glm(data = fishlakes, m11, cost, K = 42)  # use leave-one-out cross validation (can use K-fold cross validation for larger data sets)

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Now lets see what our error rate was:

m11.cv$delta

## [1] 0.2381 0.2438

That’s not too bad.

IF we want to check the error rate of an averaged model, you need to get more creative. I’ve written some code to do this for averaged models that only have two component models. It shouldn’t be too hard to adapt this for more models.

function for leave one out cross validation of averaged models

Cross <- function(model1, model2, data, cost) {
    library(MuMIn)
    nobs <- nrow(data)
    model.ave <- model.avg(model1, model2)
    values <- matrix(NA, nrow = nobs, ncol = 5)
    values[, 1] <- data$treatment
    values[, 2] <- predict(model.ave, type = "response")
    CV = 0
    for (i in 1:nobs) {
        data2 <- data[-i, ]  # leave out one observasion
        model12 <- glm(model1$formula, family = model1$family, data = data2)
        model22 <- glm(model2$formula, family = model2$family, data = data2)
        model.ave2 <- model.avg(model12, model22)
        values[i, 3] <- predict(model.ave2, newdata = data[i, ], type = "response")
        values[i, 4] <- round(values[i, 3])
        if (values[i, 4] == values[i, 1]) 
            (values[i, 5] = 1) else values[i, 5] = 0
        
    }
    cv = mean(abs(values[, 3] - values[, 1]) > 0.5)
    return(cv)
}

Use the function on the component models

Cross(m11, m4, fishlakes, cost)

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## [1] 0.2381

That looks reasonable

Another method to test model accurace is Area Under the Reciever Operater Curve (AUC) This is baisically a plot of true presences versus false presences in a presence-absense model.

Load the library “pROC”

library(pROC)

Make your reciever-operater curve

m.roc <- roc(fishlakes$treatment, predict(m.ave, backtransform = TRUE))
plot(m.roc)

## 
## Call:
## roc.default(response = fishlakes$treatment, predictor = predict(m.ave,     backtransform = TRUE))
## 
## Data: predict(m.ave, backtransform = TRUE) in 24 controls (fishlakes$treatment 0) < 18 cases (fishlakes$treatment 1).
## Area under the curve: 0.949

Looks like a pretty good fit. Not too bad for the small size of the data set.

And that’s all I got. Hopefully you will find it helpful. Let me know if there is anything else I have forgotton or done wrong.

~ Rosemary

rosehartman@gmail.com

Mason Earles on interfacing R with the Forest Vegetation Simulator

2013-02-13T15:06:42-08:00

Mason Earles gave a great presentation this week at Davis R Users’ Group about linking R with the Forest Vegetation Simulator (FVS). FVS is a model developed by the US Forest Service to simulate forest growth over time. It’s written in FORTRAN and has been around since the 1970s.

FVS has recently gone open-source (its repository is on google code), and now has an alpha-level API which can be called from R. Mason showed us how he used the API along with some text manipulation functions to make the programs work together. Then he showed us how he used this to do some Monte Carlo simulations to look at the effects of stochastic fire and drought events.

Here are Mason’s slides:

Ryan Peek on using xts and ggplot for time-series data

2013-02-06T15:44:25-08:00

At Davis R Users’ Group today, Ryan Peek gave a presentation on how he takes data from his field instruments and visualizes it in R. Here are his notes. The original *.Rmd file and data can be found here

SHORT HOW-TO ON USING XTS AND GGPLOT FOR TIME SERIES DATA

XTS is a very helpful package when working with time series data. I work with temperature and flow data frequently, so the ability to work with timeseries, and particularly to shift intervals (from 15 min to hourly or daily) can be very handy. Once data has been re-worked, ggplot2 is a great package to plot data of all types.

First read in some data and format to POSIXlt. There are 2 types of time, POSIXlt and POSIXct. POSIXct is generally what ggplot uses and is based on seconds, while POSIXlt can be re-classed into individual components quite easily, but year zero is calculated from 1900 (so 2011 would be year 111).

Read in some data

## READ IN FILE
inputfile <- "Air_Water_loggers_15min_example.csv"
file <- read.csv(inputfile, stringsAsFactors = FALSE)
## FORMAT TIME COLUMN TO POSIXlt
file$time15 <- strptime(file$time15, format = "%m/%d/%Y %H:%M")
## ROUND TO NEAREST HOUR
file$time15hr <- round(file$time15, "hour")

## TO VIEW TIME CLASSES USE 'str'
str(unclass(file$time15))

## List of 9
##  $ sec  : num [1:6096] 0 0 0 0 0 0 0 0 0 0 ...
##  $ min  : int [1:6096] 15 30 45 0 15 30 45 0 15 30 ...
##  $ hour : int [1:6096] 0 0 0 1 1 1 1 2 2 2 ...
##  $ mday : int [1:6096] 30 30 30 30 30 30 30 30 30 30 ...
##  $ mon  : int [1:6096] 3 3 3 3 3 3 3 3 3 3 ...
##  $ year : int [1:6096] 112 112 112 112 112 112 112 112 112 112 ...
##  $ wday : int [1:6096] 1 1 1 1 1 1 1 1 1 1 ...
##  $ yday : int [1:6096] 120 120 120 120 120 120 120 120 120 120 ...
##  $ isdst: int [1:6096] 1 1 1 1 1 1 1 1 1 1 ...

Notice that the $time15 column classes exist for all values even though data only exists to minutes. Also years are “112” and “111” instead of 2011 and 2012. Month and hour (24 hour clock) start at zero.

Load the XTS package and make into an xts object

library(xts)
## MAKE XTS

## make XTS format using time column for the data columns of interest
df.xts <- xts(x = file[, c(2:7)], order.by = file[, "time15"])
head(df.xts)  # notice the row names are now the 'time15' column

##                     RUB_sol MFA_sol NFA_sol NFY_sol SFY_baro_air
## 2012-04-30 00:15:00   11.93   12.69   11.26    8.19         9.43
## 2012-04-30 00:30:00   11.88   12.63   11.28    8.18         9.25
## 2012-04-30 00:45:00   11.83   12.57   11.26    8.21         9.03
## 2012-04-30 01:00:00   11.78   12.51   11.23    8.22         8.82
## 2012-04-30 01:15:00   11.73   12.46   11.17    8.23         8.60
## 2012-04-30 01:30:00   11.68   12.40   11.15    8.24         8.42
##                     NFA_baro_air
## 2012-04-30 00:15:00        13.13
## 2012-04-30 00:30:00        12.93
## 2012-04-30 00:45:00        12.74
## 2012-04-30 01:00:00        12.61
## 2012-04-30 01:15:00        12.46
## 2012-04-30 01:30:00        12.34

Notice that the row names exist as the time series data from the time15’ column.

Then run any function along a given time interval using XTS functions

apply.daily - makes daily data
apply.monthly - monthly
period.apply(data,endpoints(data,on=“interval”,1)) # Using interval of your choice

When using Rstudio, you can hit tab at any time to view potential options for a given command or function. Try apply+TAB.

# MAKE DAILY AVG, MEAN, AND MAX
dAvg <- apply.daily(df.xts, function(x) apply(x, 2, mean))  # daily mean
dMin <- apply.daily(df.xts, function(x) apply(x, 2, min))  # daily min
dMax <- apply.daily(df.xts, function(x) apply(x, 2, max))  # daily max

# OR MAKE HOURLY (OR PREFERRED TIME INTERVAL)
hrAvg <- period.apply(df.xts, endpoints(df.xts, on = "hours", 1), function(x) apply(x, 
    2, mean))  # avg over 1 hour period

# RENAME COLUMN NAMES TO MAKE RECOMBINING EASIER
names(dAvg) <- sub("_sol$", ".avg", names(dAvg))
names(dAvg) <- sub("_air$", ".avg", names(dAvg))
names(dMax) <- sub("_sol$", ".max", names(dMax))
names(dMax) <- sub("_air$", ".max", names(dMax))
names(dMin) <- sub("_sol$", ".min", names(dMin))
names(dMin) <- sub("_air$", ".min", names(dMin))

Day <- data.frame(datetime = index(dAvg), dAvg[, c(1:6)], dMax[, c(1:6)], dMin[, 
    c(1:6)], row.names = NULL)

# WRITE NEW FILE TO A CSV write.csv(Hour,
# file='Airwatertemps_hourly_example.csv', row.names=FALSE)

When making a dataframe from an xts object, important to note that you must use index() instead of row.names to pull the time stamp/date into a new dataframe.

Make a ggplot of daily temperatures showing a max and min using geom_ribbon

library(scales)
library(ggplot2)

This plot is complex, and code is provided to view as you like. Main components are: * ggplot: this sets the data and x acess * geom_line: this is to plot lines, custom colors are set by including color within the aes call * geom_ribbon: this is to plot the filled max and min ribbon in the background *scale_color/fill_manual: this allows you to define custom colors for each geom, must match exactly with the color=“" or fill=”"

RUB <- ggplot(Day, aes(x = datetime)) + geom_line(aes(y = NFA_baro.avg, color = "AIR"), 
    size = 1, linetype = 1, alpha = 1) + geom_ribbon(aes(ymax = NFA_baro.max, 
    ymin = NFA_baro.min, fill = "AIR Max/Min"), colour = "darkorange", linetype = 3, 
    alpha = 0.4) + geom_ribbon(aes(ymax = RUB.max, ymin = RUB.min, fill = "RUB Max/Min"), 
    colour = "darkgreen", linetype = 3, alpha = 0.7) + geom_line(aes(y = RUB.avg, 
    color = "RUB-Water"), size = 0.9, linetype = 1, alpha = 1) + geom_ribbon(aes(datetime, 
    ymax = MFA.max, ymin = MFA.min, fill = "MFA Max/Min"), colour = "gray40", 
    linetype = 3, alpha = 0.5) + geom_line(aes(y = MFA.avg, color = "MFA-Water"), 
    size = 0.9, linetype = 1, alpha = 1) + ylab(expression(paste("Daily Temperature ( ", 
    degree, "C)", sep = ""))) + xlab("") + scale_x_datetime(limits = c(as.POSIXct("2012-05-10 12:00"), 
    as.POSIXct("2012-05-30 12:00")), labels = date_format("%b-%d"), breaks = date_breaks("7 days")) + 
    scale_y_continuous(breaks = seq(0, 42, by = 4), limits = c(0, 42), labels = seq(0, 
        42, by = 4)) + theme_bw() + labs(title = "Daily Air and Water Temperature: Rubicon") + 
    scale_color_manual("Mean", breaks = c("AIR", "RUB-Water", "MFA-Water"), 
        values = c(`RUB-Water` = "darkgreen", `MFA-Water` = "darkmagenta", AIR = "orange")) + 
    scale_fill_manual(" Max/Min", breaks = c("AIR Max/Min", "RUB Max/Min", "MFA Max/Min"), 
        values = c(`AIR Max/Min` = "gold", `MFA Max/Min` = "darkmagenta", `RUB Max/Min` = "darkolivegreen3")) + 
    theme(legend.direction = "horizontal", legend.position = "bottom", legend.key = element_blank(), 
        legend.background = element_rect(fill = "white", colour = "gray30")) + 
    guides(fill = guide_legend(keywidth = 0.9, keyheight = 1))

RUB

A Draft Dissertation Proposal: Forest Disease Dynamics and Management

2013-01-27T19:37:16-08:00

A first draft, just out to my committee. Feedback welcome! A modifiable souce version can be found on github here

Abstract

Forest disease spreads through plant communities structured by species composition, age distribution, and spatial arrangement. I propose to examine the consequences of the interaction of these components of population structure Phytophthora ramorum invasion of California redwood forests. First, I will compare the dynamic behavior of a series of epidemiological models that include different combinations of population structure. Then I will fit these models to time-series data from a network of disease monitoring plots to determine what components of forest population structure are most important for prediction of disease spread. Using the most parsimoniuous models, I will determine optimal schedules of treatment to minimize the probability of precipitous decline of host populations.

Introduction: Modeling Disease in Complex Populations for Management

Forest diseases can radically transform ecosystems. In North America, chestnut blight, Dutch elm disease and beech bark disease have caused precipitous declines in their hosts, leading to changes in the structure and function of forests. Changes from such diseases may be the dominant force changing the face of some forests in coming decades, overwhelming other forms of rapid environmental change such as climate change (Lovett et al. 2006).

Even simple host-disease systems exhibit complex dynamics (Kermack and Mckendrick 1927). Forest disease dynamics include additional complexity driven by pathogen life cycles, variation in pathogen and host populations, spatial strucure, and demographic and environmental stochasticity (Hansen and Goheen 2000, Holdenrieder et al. 2004). These factors yield greater complexity in disease dynamics. This complexity, combined with limited monitoring data and the uncertainty of the parameters of emergent diseases, makes prediction of disease dynamics difficult. Yet such prediction is needed to allocate limited resources for treatment and prevention.

One source of complexity in forest disease dynamics is structured variation in host populations. Host tree populations may consist of multiple species or genotypes, and consist of multiple sizes, each of which may interact with disease or each other differently (Gilbert and Hubbell 1996). All forest tree populations are spatially structured; non-motile individuals can not be “well-mixed” (Filipe and Maule 2003). Each of these forms of population structure have consequences for population dynamics that have been examined in theoretical models (Park et al. 2001, 2002, Dobson 2004, Klepac and Caswell 2010). However, dynamics resulting from these components of population structure interacting are not well-characterized.

Understanding modeled dynamics of disease in structured populations may provide insights useful for management. However, the ability to use such models for management planning depends requires confidence in their predictive power, particularly under novel management treatments. Recent developments in particle filter techniques allow estimation of the likelihoods of complex dynamic models (Arulampalam et al. 2002, Ionides et al. 2006, Knape and de Valpine 2012), and subsequent comparison of models in out-of-sample predictive performance (Vehtari and Ojanen 2012). Fitting dynamic models direclty provides more confidence in the dynamic results than deriving madel parameters individually from data, because allows tests confidence in models’ emergent dynamics rather than separate components.

Sudden oak death (SOD) is an emerging forest disease in California and Oregon that threatens populations of tanoak (Notholithocarpus densiflorus) (Rizzo and Garbelotto 2003), has the potential to modify community structure (Metz et al. 2012), and cause significant economic damage (Kovacs et al. 2011). Silvicultural, chemical, and other control techniques can modify the progression of disease (Swiecki and Bernhardt 2013), but eradication of SOD is unlikely [@Cobb2013] even at local scales. Long and short- and long-term solutions to limit the damage of disease will require continuous treatment and monitoring regimes. Under budget constraints, planning for such treatment can be informed by dynamics models via optimal control techniques.

I aim to answer the following overall questions:

How do interactions between community, size, and spatial structure affect disease dynamics in forests?
What aspects of population structure are most important in determining the probability and size of forest disease outbreaks in theoretical models?
What aspects of population structure are most important in predicting the probability and size of SOD outbreaks?
What control strategy minimizes the probability of SOD outbreaks at least cost?

Background

Sudden Oak Death

Sudden Oak Death (SOD) is an emerging forest disease in California and Oregon that poses risks to forests across North America and Europe. First observed in California in the mid-1990s, SOD is caused by the water mold Phytophthora ramorum (Rizzo et al. 2002). P. ramorum often kills tanoak (Notholithocarpus densiflorus), which provide habitat and food to many vertebrate species, and are the primary host of ectomycorrhizal fungi in redwood forests (Rizzo and Garbelotto 2003). Loss of this species may have cascading effects on other species. The disease has also caused significant economic damage through removal costs and property value reduction (Kovacs et al. 2011). It has the potential to spread to species in other regions, such as northern red oak (Quercus rubra), one of the most important eastern timber species (Rizzo et al. 2002).

Population Structure: Multiple Host Species

Phytophthora ramorum has over 100 host species in 40 genera, which fall into several functional types (Swiecki and Bernhardt 2013). In canker hosts, which are all members of Fagaceae in California, pathogen produces cankers on the trunk. Few if any spores are produced from the cankers, and these hosts are generally dead ends. SOD is fatal to many canker costs. In foliar hosts, P. ramorum resides and reproduces in leaves and twigs. Foliar hosts vary greatly in susceptibility to disease and spore production from diseased individuals. Some hosts are dead ends; others are major drivers of disease spread and survival.

In Redwood forests in Northern California, the hosts of importance are tanoak and bay laurel (Umbellularia californica). Tanoak has the dubious distinction of being the only species that is both a canker and foliar host, and the disease is generally lethal in these trees. Bay laurel is a foliar host in which P. ramorum produces prolifically, and bay suffers no harm from the disease (DiLeo et al. 2009). Redwood and other species are unsusceptible or dead-end hosts where P. ramorum has little detrimental effect. In forests of the redwood-tanoak-bay complex, bay Laurel acts as the primary disease reservoir [@Davidso2008]. Tanoak infection and mortality is dependent on bay density (Cobb et al. 2012), and removal of bay is a commonly suggested treatment to reduce risk of SOD outbreaks (@Filipe2013; Swiecki and Bernhardt 2013).

Dobson (2004) examined a model of disease in multiple host species based on an the SIR (Kermack and Mckendrick 1927) framework. In it, species had density dependence but interacted only though disease transmission. Importantly, transmission between species was always equal or less than transmission between species. Dobson found that in when transmission was density-dependent, species had a complementary effect on the disease’s rate of reproduction ; that is, a mix of both species led to faster outbreaks. However, in the frequency-dependent case, was greatest when one species dominated. Transmission between species dampens and synchronizes disease oscillations in each species, and at very high contact rates between species the most susceptible species are more like to go extinct. Craft and Hawthorne (2008) implemented a stochastic version of this model in a system of mamallian carnivores, finding that species with low intra-species contact rates acted as sinks for disease, and their contact with more social species increased spread disease within these low-contact species.

SOD dynamics differ from these systems in a few important ways. First, the contact structure among species may not be symmetrical. Contact rates among trees are driven by differences in species transmissivity (amount of spore produced) and susceptibility (probability of infection per spore) than social grouping. Where one species with high transmissivity coexists with a species of high susceptibility, the probability of transmission may be higher between species than within. Also, the spatial structure of species, and the dispersal pattern of spores, may result in a hybrid between density and frequency-dependent disease transmissions.

Population Structure: Size/Stage Classes

Within the tanoak population, epidemiological characteristics vary with tree size. Trees of larger size classes are more likely to be infected and die more quickly than smaller tree (Cobb et al. 2012), Possibly due to the vulnerability of cracked bark and the amount of bark tissue available for invasion (Swiecki and Bernhardt 2005). There is no evidence currently for variation in spore production or other physiological effects across tree sizes [@Davidson2008].

When disease dynamics occur much faster than demographic processes, we can treat size classes similarly to different species, and examine the progression of disease in a constant population structure. However, when disease progression and growth occur at similar rates, these processes can interact to produce complex dynamics (Klepac and Caswell 2010). This is the case with SOD; infectious periods can last many years, during which trees may continue to grow.

The effects of age-based population structure and contact rates have been studied extensively in the context of human disease and vaccination programs (Anderson and May 1985, Metcalf et al. 2012). Reducing in-class transmission rates at young ages increases the average age of infection, independent of the exact contact structure. When only the susceptibility, and not transmissivity, of each age class varies. Klepac and Caswell (2010) found that increases in within-class transmission increases both the infected and recovered population of those classes, and classes with high transissivity have higher equilibrium populations relative to low-transmissivity classes.

Population Structure: Space

P. Ramorum dispersal interacts with host population structure at many scales. The pathogen spreads between trees via wind-blown rain and splash, limiting most spores to spread within 15m of host plants (Davidson et al. 2005). However, occaisional weather events, such as fog, can transport spores up to 3 km (Rizzo et al. 2005). Over longer distances, P. ramorum can be transported in streams or spread via human-mediated vectors such as nursery plant trade [@Osterbauer2004]. Meentemeyer et al. (2011) found that the best fit kernel for a statewide model of P. ramorum spread was the sum of two Cauchy kernels, one on the scale of tens of meters and another on the scale of tens of kilometers.

Despite the ability to spread long distance, strong meter-scale gradients mean that individual trees in stands can not be considered “well mixed” in terms of contact rates between trees. Rather, contact patterns across space arise from the interaction of spatial clustering of trees and the dispersal kernel of the disease. Neither frequency- nor density-dependent transmission characterizes this arrangement.

The effect of such spatial structure on development of epidemics has been studied with continuous populations in space [@Bolker1999], metapopulation models (Park et al. (2001); Park et al. (2002)), models of individuals on a lattice (Filipe and Maule 2003, 2004) and discrete individuals in continuous space (Brown and Bolker (2004)). These approaches have reached similar conclusions of the effect of spatial structure in non-mobile populations. In all cases, the threshold of disease growth for a global outbreak is greater than the threshold for a local epidemic . For instance, in metapopulations on a lattice, , where is the average number of connections between patches and is the inter-patch contact rates (Park et al. 2001). In models of discrete individuals in space, , where would be the pathogen growth rate in a well-mixed population, and is the dispersal kernel-weighted spatial correlation between susceptible and infected individuals. evolves over time but reaches a minimum that represents the threshold required for a global epidemic (Brown and Bolker 2004). Also, both clustering and anti-clustering (oversdispersal) reduce the rate of increase of disease, except when clustering occurs at the scale of dispersal, in which case it can accellerate spread. Fat-tailed dispersal kernels accellerate spread, as well.

Study Approach

I will examine the comparative effects of the above components of population structure on SOD epidemics by (1) characterizing the dynamic behavior of models that include each type of structure and their combinations and (2) identifying the model structure that best predicts data of disease spread over time in redwood-tanoak-bay forests. Using the best-fit model, I determine optimal control plans for stand protection efforts to minimize the risk of loss of host plants.

Determining the Consequences of Population Structure via Comparative Dynamics

1. How do interactions between community, size, and spatial structure affect disease dynamics in forests?

2. What aspects of population structure are most important in determining the probability and size of forest disease outbreaks in theoretical models?

To characterize the effects of the population structure on disease dynamics, I will compare four models. All are extensions of the epidemiological model of Cobb et al. (2012).

Model A is the simplest of the four models, only representing structure in the form of differences between species of trees. It is described by this system of stochastic difference equations:

\[\begin{aligned} \boldsymbol{S}_{t+1} &\sim \boldsymbol{S_t} + \overbrace{ \text{Pois}\left[\boldsymbol{b(S_t+I_t)}\left(1-\sum_{i=1}^n w_i (S_{it} + I_{it}) \right) + \boldsymbol{rm_I I_t}\right]}^{\text{new recruits}} - \overbrace{\text{Binom}_1\left((\underbrace{1 - e^{-\boldsymbol{\beta I_t} - \boldsymbol{\lambda_{ex}}}}_{\text{force of infection}})\boldsymbol{S_t} \right)}^{\text{infections}} -\overbrace{\text{Binom}_2(\boldsymbol{m_S S_t})}^{\text{mortality}} \\ \boldsymbol{I}_{t+1} &\sim \boldsymbol{I_t} + \text{Binom}_1\left((1 - e^{-\boldsymbol{\beta I_t} - \boldsymbol{\lambda_{ex}}})\boldsymbol{S_t} \right) - \text{Binom}_3(\boldsymbol{m_I} \boldsymbol{S_t}) \end{aligned}\]

Here and are vectors of the populations of susceptible and infected individuals of each species, and is the matrix of contact rates. Other parameters are listed below in Table 1.

Model A Parameters
Parameter Vector Symbol	Description
	Probability of death per year, for both susceptible and infectious trees
	Fecundity per individual per year.
	Probability of resprouting after death by disease.
	Competitive coefficient (relative contribution to density-dependent recruitment)
	Rate of contact of each species with pathogen spores from outside the site

New susceptible trees enter the system via density-dependent seedling recruitment and resprouting from recently killed trees. Susceptible trees become infected at rates proportional to contact with infected trees (density-dependent), and pathogen migrating into the system. Both susceptible and infected trees die at constant, but different, rates.

The model is modified from Cobb et al. (2012) in several ways: (1) inclusion of demographic stochasticity, (2) conversion from continuous to discrete time, (3) the inclusion of , the force of infection from areas outside the system, and (4) the exclusion of recovery of infected trees, which has not been observed in the field. Finally, the only population structure in Model A is the difference in species parameters. It excludes age structure and spatial structure, assuming a well-mixed population and frequency-dependent transmission (Following Diekmann et al. (1990))

Model B adds stage structure to Model A. In this case, the vectors and represent the population divided by species and size class as follows

\[\begin{aligned} S'_t &\sim S_t - \text{Binom}_1\left((1 - e^{-\boldsymbol{\beta I_t} - \boldsymbol{\lambda_{ex}}})\boldsymbol{S_t} \right) \\ I'_t &\sim I_t + \text{Binom}_1\left((1 - e^{-\boldsymbol{\beta I_t} - \boldsymbol{\lambda_{ex}}})\boldsymbol{S_t} \right) \\ S_{t+1} &\sim \text{Multinom}(\boldsymbol{A_S(S'_t,I'_t)S'_t}) + \text{Pois}\left[\boldsymbol{b(S'_t+I'_t)}\left(1-\sum_{i=1}^n w_i (S'_{it} + I'_{it}) \right) + \boldsymbol{rm_I I'_t}\right] \\ I_{t+1} &\sim \text{Multinom}(\boldsymbol{A_I(S'_t,I'_t)I'_t}) \end{aligned}\]

Here is the matrix of demographic rates specifying transitions between classes and mortality. is a block-diagonal matrix of size transition matrices, with no transitions between species classes.

Disease transitions are separated in time from demmographic transitions (Klepac and Caswell 2010). In California, P. Ramorum reproduces and spreads during the winter rainy season while most tree growth occurs in the spring and fall.

Models 3 and 4 modify Models 1 and 2 by adding spatial structure in the form of a lattice metapopulation. and become and , matrices of each species and class at each location. Furthermore, the force of infection is replaced with an overall force of infection at each location

Here is the location, is the dispersal kernel of P. ramorum from all locations to , and is a vector of parameters of the dispersal kernel. remains the contact matrix, representing physiological and vertical dispersal components of spore dispersal. I will test the models using normal, exponential, and power-law disperal kernels.

Using these four models, I will determine:

How global epidemic growth rate , probability of global epidemic, the epidemic size, and time to extinction of both hosts and pathogens vary between the models
How the presence of each form of population structure affect these values when global species densities and mean parameter values remain identical
Whether the effects of population structure are additive or if they interact in a more complex way.

The model will be parameterized with three species: tanoak, bay, and redwood, (the latter representing all non-host species). Parameters and initial conditions will be drawn from Cobb et al. (2012) and @Filipe2013.

Selecting the Model that Best Repesents Observed Behavior

3. What aspects of population structure are most important in predicting the probability and size of SOD outbreaks?

To apply the understanding developed in exploring model structure to management, I must determine which model best represents observed dynamics of SOD in forests. I will fit the four models above to data of SOD spread in redwood-tanoak-bay forests in California and compare their ability to predict disease development

Data

Data to fit the models comes from a collaboration with the Rizzo lab’s disease monitoring plot network [@Ref; @Ref; 2012]. The data consist of observations tree size, alive/dead status, and disease status for trees greater than 1 cm diameter from the years 2002-2007. Trees were observed in 14 sites along the California coast from X (36.16°N) to Y (38.35°N), each of which contains an average of 8 500 m², seperated by minimium distances of 100m and arranged in a random linear sequence. All plots are in forests dominated by redwood, tanoak, and bay laurel.

Size, health, and disease status of all trees were measured in 2002 and 2007, as well as for a random sample of 5 previously infected and 5 previously uninfected trees of each species per plot in the years 2003-2006.

Since the time scale of the are not suffiently long to data to determine transition between size classes, the parameters of transition matrices will be drawn from other data from other sites and included as strong priors on these parameters.

The Observation Model

The process generating the data above may be considered hidden Markov process (Gimenez et al. 2012), that is, processes where the state of the system is dependent on its previous state, but where states are partially or imperfectly observed. Imperfect observation in this case is due to incomplete counting of trees and error in determination of tree disease status. Tree disease is determined by subjective examination of symptoms in the field, and then confirmed in the lab for those individuals with sufficiently suggestive symptoms. There are both false positives and false negatives:

\[\begin{aligned} p_{OI} &= p(\text{Observing that tree is diseased}|\text{tree is diseased}) \\ p_{OS} &= p(\text{Observing that tree is healty}|\text{tree is healthy}) \end{aligned}\]

Then, for each sub-population (by species, size class, or location) observed in each time period, then, the observed number of susceptible and infected trees are:

\[\begin{aligned} S_{obs} \sim \text{Binom}\left(p_{OS} S + (1-p_{OI}) I \right) \\ N_{obs} \sim \text{Binom}\left(p_{OI} I + (1-p_{OS}) S \right) \end{aligned}\]

OMITTED HERE: OBSERVATION MODEL FOR PARTIAL CENSUS YEARS THAT I REALIZED WAS WRONG

Finally, for models with size-class structure, the observations omit the smallest size class (seeds and seedlings < 1 cm), as these were not measured.

Model Fitting

Particle filtering determines the likelihood of a hidden Markov process by taking advantage of the fact that, at each time step, the likelihood of a system state is dependent on the previous state , and the observation .

While is hidden, may be approximated by by simulation of the model process. Repeated simulations (particles) provide a distribution of outcomes from which average likelihood are determined from the observation model . For each time step, the particles are re-sampled according to their individual likelihoods, in order to prevent particle degeneracy - most particles approaching zero likelihood. The product of all steps’ averaged particle likelihoods is an unbiased estimate of the model likeihood.(Arulampalam et al. 2002)

With an estimate of likelihood available, one needs an approach to determine the maximum-likelihood estimate of the model parameters . Several methods are available. Iterated filtering (IF) estimates by replacing constant parameter values with a random walk of decreasing variance over time, and determining the “states” of this walk as the variance approaches zero. IF is computationally efficient but requires long time series.

Alternatively, Particle Filter Markov Chain Monte Carlo (PFMCMC) uses a Markov Chain sampler (e.g., Metropolis-Hastings), calculating the likelihood at each iteration using a particle filter. While computationally expensive, it does not require long time series and may be used with multiple time series. PFMCMC also integrates particle filtering into a Bayesian context, allowing the use of informative priors. I will used PFMCMC to determine maximum-likelihood parameterizations for all four models.

When fitting models 1 and 2 to the data, populations from each plot in each site will be aggregated. For models 3 and 4, each plot will represent a single sub-population on a lattice spanning the site. Initial density, species, age and disease distributions for other sub-populations will be randomly generated from other points on the lattice using the proportions from the measured plots. Epidemiological parameters will remain constant across sites except for the external force of disease , which will be represented as a noramlly distributed random variable to account for site effects. varies.

Model selection

Model selection criteria have multiple purposes: to determine which model best represents “true” processes, and to estimate the utility of models for the purposes of prediction. In this case, I am interested in both - determining which components population structure matter, and also providing an estimate of predictive stength of the best model. For these purposes, I will use the deviance information criterion (Spiegelhalter et al. 2002, DIC). DIC has several advantages. Since it estimates the model complexity, or effective number of parameters, directly from the likelohood distribution, it makes it easy to incorporate random parameters such as that have non-integer contributions to model complexity. Secondly, it estimates the how priors on parameter values effect their contribution to model complexity.

DIC is not suitable for model-averaging. However, model-averaged results will not likely be useful for the optimization analysis below. Using the best-DIC model for prediction is a compromise between predictive ability and tractable analysis of the model used for prediction. Use of DIC assumes that there is a “best” or “real” model within the comparison. While no models are “true” (Box 1976), this excercise is selecting a “best” or “most useful” model for scenario exploration.

Using Dynamic Optimization to Determine the Best Control Strategy

4. What control strategy minimizes the probability of SOD outbreaks at least cost?

The best-fit model that captures essential aspects of the disease dynamics will be used to develop an optimal control solution to minimization of disease risk over time.

While forest disease management can occur under multiple, conflicting and uncertain goals and priorities, here I consider the case where conservation of tanoak and its ecosystem services are the goal. This is approximately the case in Humboldt County, California, where large landowners have agreed to work in tandem to control the disease. In this situation, managers goals can be simplified to three components: (1) The minimum acceptable density of the species of interest , (2) the cumulative acceptable probability that the species will fall below that density during the management period, and (3) the cost of disease control . Using acceptable probabilities avoids the challenge of defining and valuing the flows of ecosystem services as a function of the forest composition (Hartman 1976).

Tools for SOD control can be categorized into two types. Silvicultural (“host-centric”) treatments change forest composition. “Pathogen-centric” treatments, such as quarantine, equipment cleaning and spraying, reduce the probability of spore arrival (Swiecki and Bernhardt 2013).

I represent these two types of treatment with two control functions. Cutting is a periodic control treatment that can be imposed in any year. I define the cut as a reduction in species density to the levels that would minimize the probability of species densities falling below the acceptable level over a management period without further controls. This will be a value that balances stochastic risk of extinction with risk of diseae outbreak, bounded by the minimum acceptable species density. As tanoak and bay laurel trees have little timber value, each cut has a cost . For the purpose of this model is independent of the extent of density reduction (that is, the variable costs of cutting are negligible. Thus, the annual cost of cutting is , where is a zero-one control variable

Spraying is an ongoing control treatment that represents the combined annual effort to reduce the arrival of new spores at the site. This modifies the external force of infection

where is the effectiveness of treatment, is the effort spent at spraying and is the annual cost per unit effort.

The dynamic optimization problem is to determine the combination of treatment effort over time that minimizes the cost of treatment while maintaining the probability of species loss at an acceptable level. That is,

An alternative formulation, also useful to managers, is to minimize the probability of outbreak given budget constraints . Budget constraints could be forumulated as total budgets for the entire mangement period, or annual limits on spending. That is

References

Anderson, R. M., and R. M. May. 1985. Age-related changes in the rate of disease transmission: implications for the design of vaccination programmes. J. Hyg. Camb.

Arulampalam, M. S., S. Maskell, N. Gordon, and T. Clapp. 2002. A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking. IEEE Transactions on Signal Processing 50:174–188.

Box, G. E. P. 1976. Science and Statistics. Journal of the American Statistical Association 71:791–799.

Brown, D. H., and B. M. Bolker. 2004. The effects of disease dispersal and host clustering on the epidemic threshold in plants. Bulletin of mathematical biology 66:341–71.

Cobb, R. C., J. A. N. Filipe, R. K. Meentemeyer, C. A. Gilligan, and D. M. Rizzo. 2012. Ecosystem transformation by emerging infectious disease: loss of large tanoak from California forests. Journal of Ecology:no–no.

Craft, M. E., and P. L. Hawthorne. 2008. Dynamics of a multihost pathogen in a carnivore community. Journal of Animal Ecology 77:1257–1264.

Davidson, J. M., A. C. Wickland, H. A. Patterson, K. R. Falk, and D. M. Rizzo. 2005. Transmission of Phytophthora ramorum in Mixed-Evergreen Forest in California. Phytopathology 95:587–596.

DiLeo, M. V., R. M. Bostock, and D. M. Rizzo. 2009. Phytophthora ramorum does not cause physiologically significant systemic injury to California bay laurel, its primary reservoir host.. Phytopathology 99:1307–11.

Diekmann, O., J. A. P. Heesterbeek, and J. A. J. Metz. 1990. On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations. Journal of Mathematical Biology 28:365–382.

Dobson, A. 2004. Population dynamics of pathogens with multiple host species.. The American Naturalist 164 Suppl:S64–78.

Filipe, J. A. N., and M. M. Maule. 2003. Analytical methods for predicting the behaviour of population models with general spatial interactions. Mathematical Biosciences 183:15–35.

Filipe, J. A. N., and M. M. Maule. 2004. Effects of dispersal mechanisms on spatio-temporal development of epidemics. Journal of Theoretical Biology 226:125–141.

Gilbert, G. S., and S. P. Hubbell. 1996. Diseases and the Tropical of Conservation Forests. BioScience 46:98–106.

Gimenez, O., J. D. Lebreton, J. M. Gaillard, R. Choquet, and R. Pradel. 2012. Estimating demographic parameters using hidden process dynamic models. Theoretical Population Biology 82:307–316.

Hansen, E. M., and E. M. Goheen. 2000. Phellinus Weirii and other Native Root Pathogens as Determinants of Forest Structure and Process in Western North America. Annual review of phytopathology 38:515–539.

Hartman, R. 1976. The Harvesting Decision when a Standing Forest has Value. Economic Inquiry 14:52–58.

Holdenrieder, O., M. Pautasso, P. J. Weisberg, and D. Lonsdale. 2004. Tree diseases and landscape processes: the challenge of landscape pathology.. Trends in ecology & evolution (Personal edition) 19:446–52.

Ionides, E. L., C. Bretó, and a a King. 2006. Inference for nonlinear dynamical systems. Proceedings of the National Academy of Sciences of the United States of America 103:18438–43.

Kermack, W. O., and A. G. Mckendrick. 1927. A Contribution to the Mathematical Theory of Epidemics. Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 115:700–721.

Klepac, P., and H. Caswell. 2010. The stage-structured epidemic: linking disease and demography with a multi-state matrix approach model. Theoretical Ecology 4:301–319.

Knape, J., and P. de Valpine. 2012. Fitting complex population models by combining particle filters with Markov chain Monte Carlo.. Ecology 93:256–63.

Kovacs, K., T. Václavík, R. G. Haight, A. Pang, N. J. Cunniffe, C. a Gilligan, and R. K. Meentemeyer. 2011. Predicting the economic costs and property value losses attributed to sudden oak death damage in California (2010-2020).. Journal of environmental management 92:1292–302.

Lovett, G. M., C. D. Canham, M. A. Arthur, C. Kathleen, R. D. Fitzhugh, and K. C. Weathers. 2006. Forest Ecosystem Responses to Exotic Pests and Pathogens in Eastern North America. BioScience 56:395–405.

Meentemeyer, R. K., N. J. Cunniffe, A. R. Cook, J. A. N. Filipe, R. D. Hunter, D. M. Rizzo, and C. a Gilligan. 2011. Epidemiological modeling of invasion in heterogeneous landscapes: spread of sudden oak death in California (1990–2030). Ecosphere 2:1–24.

Metcalf, C. J. E., J. Lessler, P. Klepac, A. Morice, B. T. Grenfell, and O. N. Bjø rnstad. 2012. Structured models of infectious disease: Inference with discrete data. Theoretical population biology 82:275–282.

Metz, M. R., K. M. Frangioso, A. C. Wickland, R. K. Meentemeyer, and D. M. Rizzo. 2012. An emergent disease causes directional changes in forest species composition in coastal California. Ecosphere 3:art86.

Park, A. W., S. Gubbins, and C. A. Gilligan. 2002. Extinction times for closed epidemics: the effects of host spatial structure. Ecology Letters 5:747–755.

Park, A. W., S. Gubbins, and C. a Gilligan. 2001. Invasion and persistence of plant parasites in a spatially structured host population. Oikos 94:162–174.

Rizzo, D. M., M. Garbelotto, J. M. Davidson, G. W. Slaughter, and S. T. Koike. 2002. Phytophthora ramorum and sudden oak death in California: I. Host relationships. Pages 733–740 in Fifth Symposium on Oak Woodlands: Oaks in California’s Changing Landscapes. . USDA Forest Service, San Diego, CA.

Rizzo, D. M., M. Garbelotto, and E. M. Hansen. 2005. Phytophthora ramorum: integrative research and management of an emerging pathogen in California and Oregon forests. Annual review of phytopathology 43:309–35.

Rizzo, D. M., and M. Garbelotto. 2003. Sudden oak death: endangering California and Oregon forest ecosystems. Frontiers in Ecology and the Environment 1:197–204.

Spiegelhalter, D. J., N. G. Best, B. P. Carlin, and A. van der Linde. 2002. Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 64:583–639.

Swiecki, T. J., and E. Bernhardt. 2005. Disease risk factors and disease progress in coast live oak and tanoak affected by Phytophthora ramorum canker (sudden oak death). Pages 383–411 in Proceedings of Sudden Oak Death Science Symposium II. PJ Shea and M. Haverty, eds. Pacific Southwest Research Station, Forest Service, US Department of Agriculture, Albany, CA.

Swiecki, T. J., and E. A. Bernhardt. 2013. Managing sudden oak death in California: before, during, and after Phytophthora ramorum invasion. . U.S. Department of Agriculture; Forest Service; Pacific Southwest Research Station.

Vehtari, A., and J. Ojanen. 2012. A survey of Bayesian predictive methods for model assessment, selection and comparison. Statistics Surveys 6:142–228.

Noam Ross

Improved R Profiling Summaries

FasteR! HigheR! StrongeR! - A Guide to Speeding Up R Code for Busy People

First, Ask “Why?”

(1) Do you want your code to be readable?

(2) Do you want your code to be sharable?2

(3) Do you have other shit to do?

Get a Bigger Computer

Parallel Processing

Compiling Functions

Diagnosis: Profiling and Benchmarking

Benchmarking

Profiling

Misc profiling things

Find better packages

Improving your code

Vectorization

Pre-allocating memory

Using simpler data structures

The small stuff

Extra parentheses

Standard and internal functions

Other Languages: More reward for more effort

Setting up Automated Reprint Requests in Gmail

My first peer-reviewed publication: Economics and Ecology of Open-Access Fisheries

Economics and Ecology of Open-Access Fisheries

Abstract

Debugging Tools in R with Michael Hannon

Introduction

Focus on five functions:

Severity level

Debugging tools

Print “call stack” with traceback

The debug function

Explicit Calls to Browser

Inserting Code with trace

Editor’s Addendum: A quick way to enable bugging everywhere

Demographic analysis using the `popbio` library and some other fun stuff

Demographic analysis using the popbio library and some other fun stuff

1. ESTIMATING FECUNDITY

Part 1: How to estimate seed production per plant?

Part 2: How many seeds makes an emergent juvenile?

2. PICK A STARTING POPULATION VECTOR

3. GENERATE A TRANSITION MATRIX

4. GENERATE ANNUAL MATRICES — THE SIMPLE WAY for 3 easy years

NAME STAGES

SET ITERATIONS

5. DETERMINISTIC ANALAYSIS OF THESE 3 WELL BEHAVED YEARS - basic non-stochastic PVA

6. RUN STOCHASTIC ANALYSES ON THESE 3 WELL BEHAVED YEARS

Fun with stochastic analyses!

7. QUASI EXTINCTION BASED ON THESE 3 WELL BEHAVED YEARS

8. SENSITIVITY AND ELASTICITY

APPENDIX A. HOW I GENERATED A FECUNDITY ESTIMATE

Model Selection and Multi-Model Inference

Mason Earles on interfacing R with the Forest Vegetation Simulator

Ryan Peek on using xts and ggplot for time-series data

SHORT HOW-TO ON USING XTS AND GGPLOT FOR TIME SERIES DATA

Read in some data

Load the XTS package and make into an xts object

Then run any function along a given time interval using XTS functions

Make a ggplot of daily temperatures showing a max and min using geom_ribbon

A Draft Dissertation Proposal: Forest Disease Dynamics and Management

Abstract

Introduction: Modeling Disease in Complex Populations for Management

Background

Sudden Oak Death

Population Structure: Multiple Host Species

Population Structure: Size/Stage Classes

Population Structure: Space

Study Approach

Determining the Consequences of Population Structure via Comparative Dynamics

Selecting the Model that Best Repesents Observed Behavior

Data

The Observation Model

Model Fitting

Model selection

Using Dynamic Optimization to Determine the Best Control Strategy

References

(2) Do you want your code to be sharable?²

Print “call stack” with `traceback`

The `debug` function

Inserting Code with `trace`