Blog Archives

Lags and Moving Means in dplyr

I’ve been a big fan of dplyr from the start. While learning any new package can be frustrating, overall dplyr is fast and consistent. The Nonstandard Evaluation (see vignette(“nse”) for details) is both a blessing and a curse, overall dply is making life with big(ish) data much easier. The rapid development makes it a challenge to keep up with new functions and capabilities, but it also makes dplyr more useful every day.

My favorite use of the week has been to create lagged terms and moving means. For example, I often need new variables like yesterday’s air temperature or the average temperature over the past 5 days or even the total amount of precipitation in the preceding 30 days. I had been using the slide function within the DataCombine package, but it was somewhat slow when dealing ~10 million rows of data. Also, the slide function was fine for creating a lag but it didn’t create a moving mean or moving sum. Here’s an example using the slide function


library(dplyr)
library(zoo)
library(DataCombine)

df = expand.grid(site = factor(seq(10)),
                    year = 2000:2004,
                    day = 1:50)
# use Poisson to make math easy to check moving means of temperature
df$temp = rpois(dim(df)[1], 5) 
# Assume rains 33% of the days and averages 5 mm each time but highly variable
df$precip = rbinom(dim(df)[1], 1, 1/3) * rlnorm(dim(df)[1], log(5), 1)

# Order by group and date
df = df[order(df$site, df$year, df$day), ]
df.slide = slide(df, Var = "temp", GroupVar = c("site", "year"), slideBy = -1, NewVar='temp.lag1')
head(df.slide, 75)
 

The slide function does give the nice reminder to

Remember to order test by site and the time variable before running.

However, I do like the consistency and speed of dplyr since I’m using it for so much. It’s also easy to do multiple things in one pipe series, such as creating a lag and a moving mean. Here’s how I’m doing it in dplyr:

# moving mean for that day and previous days (e.g. 5 represents the mean of that day and the for previous days)
df2 = df %>%
  group_by(site, year) %>%
  arrange(site, year, day) %>%
  mutate(temp.5 = rollmean(x = temp, 5, align = "right", fill = NA))
head(df2, 75)

# moving mean for the previous days not including the current day (e.g. 5 represents the mean of the 5 previous days)
df2 = df2 %>%
  mutate(temp.lag1 = lag(temp, n = 1)) %>%
  mutate(temp.5.previous = rollapply(data = temp.lag1, 
                            width = 5, 
                            FUN = mean, 
                            align = "right", 
                            fill = NA, 
                            na.rm = T))
head(df2, 75)

This easily gives the option to create a moving mean based on different starting points (that day or the previous day) using the lag term that’s created in the same pipeline. It is also easy to use dplyr and zoo (for the rollapply function) to create a 30-day rolling sum of precipitation as an indication of drought status.

df2 = df2 %>%
mutate(precip.30 = rollsum(x = precip, 30, align = "right", fill = NA))
head(df2, 75)
/

The nice thing about dplyr and magittr (where it gets the %>% function) is that all of this can be easily combined into one pipeline, as such

df2 = df %>%
  group_by(site, year) %>%
  arrange(site, year, day) %>%
  mutate(temp.5 = rollsum(x = temp, 5, align = "right", fill = NA),
         temp.5.previous = rollapply(data = temp.lag1, 
                                     width = 2, 
                                     FUN = mean, 
                                     align = "right", 
                                     fill = NA, 
                                     na.rm = T)
         precip.30 = rollsum(x = precip, 30, align = "right", fill = NA))
head(df2, 75)

Output

Source: local data frame [75 x 8]
Groups: site, year

   site year day temp     precip temp.lag1 temp.2 precip.30
1     1 2000   1    7   39.41016        NA     NA        NA
2     1 2000   2    9    0.00000         7    7.0        NA
3     1 2000   3    6  132.17521         9    8.0        NA
4     1 2000   4    4    0.00000         6    7.5        NA
5     1 2000   5    5  619.73564         4    5.0        NA
6     1 2000   6    7    0.00000         5    4.5        NA
7     1 2000   7    2  122.49375         7    6.0        NA
8     1 2000   8    3    0.00000         2    4.5        NA
9     1 2000   9    6    0.00000         3    2.5        NA
10    1 2000  10    5  365.61625         6    4.5        NA
11    1 2000  11    6    0.00000         5    5.5        NA
12    1 2000  12    5    0.00000         6    5.5        NA
13    1 2000  13    7    0.00000         5    5.5        NA
14    1 2000  14    4    0.00000         7    6.0        NA
15    1 2000  15    8    0.00000         4    5.5        NA
16    1 2000  16    7    0.00000         8    6.0        NA
17    1 2000  17    5  130.09170         7    7.5        NA
18    1 2000  18    6    0.00000         5    6.0        NA
19    1 2000  19    7  305.02369         6    5.5        NA
20    1 2000  20    4    0.00000         7    6.5        NA
21    1 2000  21    7    0.00000         4    5.5        NA
22    1 2000  22    4  459.30043         7    5.5        NA
23    1 2000  23    3  117.71606         4    5.5        NA
24    1 2000  24    5  328.37299         3    3.5        NA
25    1 2000  25    3    0.00000         5    4.0        NA
26    1 2000  26    6    0.00000         3    4.0        NA
27    1 2000  27    3 1279.23160         6    4.5        NA
28    1 2000  28    4    0.00000         3    4.5        NA
29    1 2000  29    1   88.35674         4    3.5        NA
30    1 2000  30    5    0.00000         1    2.5  3987.524
31    1 2000  31    4    0.00000         5    3.0  3948.114
32    1 2000  32    4    0.00000         4    4.5  3948.114
33    1 2000  33    3    0.00000         4    4.0  3815.939
34    1 2000  34    3    0.00000         3    3.5  3815.939
35    1 2000  35    6   63.27625         3    3.0  3259.479
36    1 2000  36    6  119.32251         6    4.5  3378.802
37    1 2000  37    9    0.00000         6    6.0  3256.308
38    1 2000  38    4 3047.85230         9    7.5  6304.161
39    1 2000  39    5    0.00000         4    6.5  6304.161
40    1 2000  40    5    0.00000         5    4.5  5938.544
41    1 2000  41    4    0.00000         5    5.0  5938.544
42    1 2000  42    7  511.26682         4    4.5  6449.811
43    1 2000  43    7    0.00000         7    5.5  6449.811
44    1 2000  44    5   48.02796         7    7.0  6497.839
45    1 2000  45    4    0.00000         5    6.0  6497.839
46    1 2000  46    1    0.00000         4    4.5  6497.839
47    1 2000  47    6    0.00000         1    2.5  6367.747
48    1 2000  48    4   56.90025         6    3.5  6424.648
49    1 2000  49    6    0.00000         4    5.0  6119.624
50    1 2000  50    7    0.00000         6    5.0  6119.624
51    1 2001   1    6    0.00000        NA     NA        NA
52    1 2001   2    6    0.00000         6    6.0        NA
53    1 2001   3    7  302.76758         6    6.0        NA
54    1 2001   4    3   52.95309         7    6.5        NA
55    1 2001   5    0  261.93515         3    5.0        NA
56    1 2001   6    6    0.00000         0    1.5        NA
57    1 2001   7    8    0.00000         6    3.0        NA
58    1 2001   8    5    0.00000         8    7.0        NA
59    1 2001   9    5   58.45233         5    6.5        NA
60    1 2001  10    4  267.24615         5    5.0        NA
61    1 2001  11    6  128.48269         4    4.5        NA
62    1 2001  12    4  330.67765         6    5.0        NA
63    1 2001  13    7  966.20692         4    5.0        NA
64    1 2001  14    8    0.00000         7    5.5        NA
65    1 2001  15    4    0.00000         8    7.5        NA
66    1 2001  16    4    0.00000         4    6.0        NA
67    1 2001  17    5    0.00000         4    4.0        NA
68    1 2001  18    3   86.79884         5    4.5        NA
69    1 2001  19    3  797.27860         3    4.0        NA
70    1 2001  20    6    0.00000         3    3.0        NA
71    1 2001  21    2    0.00000         6    4.5        NA
72    1 2001  22    3    0.00000         2    4.0        NA
73    1 2001  23    7    0.00000         3    2.5        NA
74    1 2001  24    4    0.00000         7    5.0        NA
75    1 2001  25    6    0.00000         4    5.5        NA

dplyr is great…but

I have been loving Hadley Wickham’s new dplyr package for R. It creates a relatively small number of verbs to quickly and easily manipulate data. It is generally as fast as data.table but unless you’re already very familiar with data.table the syntax is much easier. There are a number of great introductions and tutorials, which I won’t attempt to recreate but here are some links:

http://rpubs.com/justmarkham/dplyr-tutorial

http://cran.rstudio.com/web/packages/dplyr/vignettes/introduction.html

http://blog.rstudio.org/2014/01/17/introducing-dplyr/

http://datascience.la/hadley-wickhams-dplyr-tutorial-at-user-2014-part-1/

The only challenge with going to dplyr has been dealing with the non-standard evaluation (see vignette(“nse”) and the new tbl_df class that gets created. In theory the resulting dataframes created from dplyr should act both as a data.frame object and as a new tbl_df object. There are some nice new features associated with tbl_df objects but also some idiosyncrasies that can break existing code.

For example, I had a function:

stdCovs <- function (x, y, var.names) {
    for (i in 1:length(var.names)) {
        x[, var.names[i]] <- (x[, var.names[i]] - mean(y[, var.names[i]], na.rm = T)) / 
            sd(y[, var.names[i]], na.rm = T)
    }
    return(x)
}

that took the standardized covariates from one dataframe and applied them to standardize another dataframe for a particular set of covariates. This was particularly useful for doing predictions for new data using coefficients from a different dataset that had been standardized (y.standard). The code above worked when I was using data.frame objects. However, I have now started using dplyr for a variety of reasons but the code no longer works. I believe this is due to how tbl_df objects handle bracketing (e.g. y[ , var]). There is an issue on GitHub related to this.

In the end it turned out that double bracketing to extract the column works: Rather than the y[ , var] formulation it should be y[[var]]. Not too hard but the type of thing that wastes an hour when going to a new package. Despite this challenge and similar things, I highly recommend dplyr. Now hopefully the simplicity of the verbs and speed of computation can make up for the few hours spent learning and adjusting code.

Here’s my final code:

stdCovs <- function(x, y, var.names){
  for(i in 1:length(var.names)){
    x[ , var.names[i]] <- (x[ , var.names[i]] - mean(y[[var.names[i]]], na.rm = T)) / 
        sd(y[[var.names[i]]], na.rm = T)
  }
  return(x)
}

One of the best features of dplyr is it’s use of magrittr for piping functions together with the %>% symbol. It’s not shown above but it’s worth learning. Previously, you would either have a bunch of useless intermediate objects created or a long string of nested functions that were hard (impossible at times) to decipher. The “then” command (%>%) does away with this as shown here by Revolution Analytics:

hourly_delay <- filter( 
  summarise(
    group_by( 
      filter(
        flights, 
        !is.na(dep_delay)
      ), 
      date, hour
    ), 
    delay = mean(dep_delay), 
    n = n()
  ), 
  n > 10 
)

Here’s the same code, but rather than nesting one function call inside the next, data is passed from one function to the next using the %>% operator:

hourly_delay <- flights %>% 
 filter(!is.na(dep_delay)) %>% 
 group_by(date, hour) %>% 
 summarise( 
   delay = mean(dep_delay), 
   n = n() ) %>% 
 filter(n > 10)

Beautiful, right. Pipes are fantastic and combined with dplyr’s simple verbs make coding a lot easier and more readable once you get the hang of it.

Blue Jay and Scrub Jay : Using rvertnet to check the distributions in R

Daniel Hocking:

This is a great post about a great package. Unlike some other broad vertebrate databases, VertNet has good QAQC. It’s great being able to access it directly in R where the data manipulation, analysis, and even plotting will occur. Thank you to the ROpenSci team for development of the rvertnet package and Vijay Barve for the clear, simple tutorial.

Originally posted on Vijay Barve:

As part of my Google Summer of Code, I am also working on another package for R called rvertnet. This package is a wrapper in R for VertNet websites. Vertnet is a vertebrate distributed database network consisting of FishNet2MaNISHerpNET, and ORNIS. Out of that currently Fishnet, HerpNET and ORNIS have their v2 portals serving data. rvertnet has functions now to access this data and import them into R data frames.

Some of my lab mates faced a difficulty in downloading data for Scrub Jay (Aphelocoma spp. ) due to large number of records (220k+) on ORNIS, so decided to try using rvertnet package which was still in development that time. The package really helpe and got the results quickly. So while ecploring that data I came up with this case study.

So here to get data for Blue Jay (Cyanocitta cristata) which is distributed in…

View original 208 more words

RStudio Presentations

After finally getting around to to posting about knitr and Markdown just yesterday, RStudio comes out with an update that makes it even easier: http://www.rstudio.com/ide/docs/presentations/overview

P.S. I Love RStudio

Knitting beautiful documents in RStudio

Title: Knitting Beautiful Documents in RStudio

Header: Markdown

This document was written using Markdown in RStudio. RStudio is a wonderful IDE for writing R scripts and keeping track of variables, dataframes, history, packages and much more. Markdown is a simple formatting syntax for authoring web pages. It allows for easy formatting, for example

# Header 1
## Header 2
### Header 3
#### Header 4
#### Header 5

Header 1

Header 2

Header 3

Header 4

Header 5

Code can be embedded in documents with nice formatting using backticks (symbol below the tilda in the upper left of an American qwerty keyboard). For example, if you were to type the following (without the >)

>```{r}
>install.packages("lme4", repos = 'http://cran.us.r-project.org')
>```

you could run that block in R and install the lme4 package. In the markdown document it would show up as:

install.packages("lme4", repos = "http://cran.us.r-project.org")

## 
## The downloaded binary packages are in
##  /var/folders/r5/l8t2qxpj3m1_gnlf6ys7fwdm0000gq/T//RtmpiqoAcq/downloaded_packages

You should try this with the knitr package [UPDATE: You will likely have to restart RStudio after installing and loading the knitr library because it changes options in RStudio]. knitr is a fantastic package for converting the R markdown file (.Rmd) into an HTML document. The design of knitr allows any input languages (e.g. R, Python and awk) and any output markup languages (e.g. LaTeX, HTML, Markdown, AsciiDoc, and reStructuredText). I find Markdown to be the easies language to use to create documents for coursework/assignments. knitr is dynamic, so it will run your code and print the results below it. The best part is if you get new data or find an error in your code, you just re-knit the file and it’s updated a few seconds later.

In RStudio, once you install knitr a button shows up that you can click to knit the markdown file (produce your html document). When you click the Knit HTML button a web page will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)

##      speed           dist    
##  Min.   : 4.0   Min.   :  2  
##  1st Qu.:12.0   1st Qu.: 26  
##  Median :15.0   Median : 36  
##  Mean   :15.4   Mean   : 43  
##  3rd Qu.:19.0   3rd Qu.: 56  
##  Max.   :25.0   Max.   :120

These will be embedded and run in the document, but you can also run single lines or blocks (chunks) of code in RStudio to make sure they work before knitting the document. You can also embed plots, for example:

plot(cars)

unnamed-chunk-3

This is great but after kniting you are left with an HTML document. What if you want a PDF or Word Document? There may be other ways to do this, but I use pandoc to covert the HTML file into other forms such as PDF or MS Word Documents. After installing pandoc, you use it via the command line. I use a Mac, so I just open the Terminal to get to my Bash Shell and cd to the folder containing my document, then just run

pandoc Knitting_Beautiful_Documents_in_RStudio.md -o 
Knitting_Beautiful_Documents_in_RStudio.pdf

and I have a PDF. Unfortunately, the pandoc default is to have huge margins that look ridiculous. To correct for that you can save a style file in the folder with the document you’re converting. To do that you just create a text file with the line:

\usepackage[vmargin=1in,hmargin=1in]{geometry}

and save it as something like margins.sty. Then in the command line you run

pandoc Knitting_Beautiful_Documents_in_RStudio.md -o 
Knitting_Beautiful_Documents_in_RStudio.pdf -H margins.sty

Now when you view the PDF it should have 1 inch margins and look much nicer. That is about it for making attractive, dynamic documents using R, RStudio, knitr, and pandoc. It seems like a lot but after doing it once or twice it becomes very quick and easy. I have yet to use it but there is also a package pander, which allows you to skip the command line step and do the file conversion right in R (RStudio).

Using Markdown in RStudio with knitr is wonderful for creating reports, class assignments, and homework solutions (as a student, you’d really impress your instructor), but other people take it even further. Karthik Ram of rOpenSci fame has a great blog post of how to exend this system to include citations and bibliographies that allow you to ditch word.

For a nice Markdown reference list, check out this useful cheat sheet

Finally, if you want to see my Markdown code and everything that went into making this, you can find in my GitHub repository. Specifically, you can find the .Rmd, .md, margins.sty, and built PDF files.

[UPDATE: On a Mac R/RStudio/pandoc may have trouble finding your LaTeX distribution to use to build the PDF (even though using Markdown not LaTeX directly). A potential solution I found when building an R package was to specify the path to LaTeX: https://danieljhocking.wordpress.com/2012/12/15/missing-pdflatex/]

 

[UPDATE 2: From RProgramming.net I found a nice way to create a PDF without leaving R (calling pandoc from R)]

# Create .md, .html, and .pdf files
knit("File.Rmd")
markdownToHTML('File.md', 'File.html', options=c("use_xhml"))
system("pandoc -s File.html -o File.pdf")
Follow

Get every new post delivered to your Inbox.

Join 67 other followers