# R bloggers

R news and tutorials contributed by (573) R bloggers
Updated: 25 min 33 sec ago

### Comment on Overnight SPY Anomaly

Mon, 2015-11-16 12:25

This post is in response to Michael Harris’ Price Action Lab post, where he uses some simple R code to evaluate the asymmetry of returns from the day’s close to the following day’s open.  I’d like to respond to his 3 notes, which I’ve included below.

1. The R backtest assumes fractional shares. This means that equity is fully invested at each new position. This is important because it affects drawdown calculations.
2. When calculating the Sharpe ratio, the “geometric = FALSE” option must be used otherwise the result may not be correct. It took some time to figure that out.
3. The profit factor result in R does not reconcile with results from other platforms or even from excel. PF in R is shown as 1.23 but the correct value is 1.17. Actually, the profit factor is calculated on a per share basis in R, although returns are geometric.

I completely agree with the first point.  I’m not sure Mike considers the output of  SharpeRatio.annualized with geometric=TRUE to be suspect (he doesn’t elaborate).  The overnightRets are calculated as arithmetic returns, so it’s proper to aggregate them using geometric chaining (i.e. multiplication).

I also agree with the third point, because the R code used to calculate profit factor is wrong.  My main impetus to write this post was to provide a corrected profit factor calculation.  The calculation (with slightly modified syntax) in Mike’s post is:

require(quantmod)
getSymbols(‘SPY’, from = ‘1900-01-01′)
overnightRets <- na.omit(Op(SPY)/lag(Cl(SPY)) – 1)
posRet <- overnightRets > 0
profitFactor <- -sum(overnightRets[posRet])/sum(overnightRets[!posRet]) Note that profit factor in the code above is calculated by summing positive and negative returns, when it should be calculated using positive and negative P&L.  In order to do that, we need to calculate the equity curve and then take its first difference to get P&L.  The corrected calculation is below, and it provides the correct result Mike expected.

grossEquity <- cumprod(overnightRets+1)
grossPnL <- diff(grossEquity)
grossProfit <- sum(grossPnL[grossPnL > 0])
grossLoss <- sum(grossPnL[grossPnL < 0])
profitFactor <- grossProfit / abs(grossLoss) I’d also like to respond to Mike’s comment:

Since in the past I have identified serious flaws in commercially available backtesting platforms, I would not be surprised if some of the R libraries have some flaws.

I’m certain all of the backtesting R packages have flaws/bugs.  All software has bugs because all software is written by fallible humans.  One nice thing about (most) R packages is that they’re open source, which means anyone/everyone can check the code for bugs, and fix any bugs that are found.  With closed-source software, commercial or not, you depend on the vendor to deliver a patched version at their discretion and in their timing. Now, I’m not making an argument that open source software is inherently better. I simply wanted to point out this one difference.  As much as I love open source software, there are times where commercial vendor-supported software presents a more appealing set of tradeoffs than using open source software.  Each situation is different.

Categories: Methodology Blogs

### How to Search for Census Data from R

Mon, 2015-11-16 08:55

(This article was first published on AriLamstein.com » R, and kindly contributed to R-bloggers)

In my course Learn to Map Census Data in R I provide people with a handful of interesting demographics to analyze. This is convenient for teaching, but people often want to search for other demographic statistics. To address that, today I will work through an example of starting with a simple demographic question and using R to answer it.

Here is my question: I used to live in Japan, and to this day I still enjoy practicing Japanese with native speakers. If I wanted to move from San Francisco to a part of the country that has more Japanese people, where should I move?

Step 1: Find the Table for the Data

Data in the census bureau is stored in tables. One way to find the table for a particular metric is to use the function ?acs.lookup from the acs package. (Note that to run this code you will need to get and install a census API key; I explain how to do that here).

> library(acs) > acs.lookup(keyword = "Japanese", endyear = 2013) An object of class "acs.lookup" endyear= 2013 ; span= 5 results: variable.code table.number table.name variable.name 1 B02006_009 B02006 Asian Alone By Selected Groups Japanese 2 B16001_069 B16001 Language Spoken at Home by Ability to Speak English for the Population 5+ Yrs Japanese: 3 B16001_070 B16001 Language Spoken at Home by Ability to Speak English for the Population 5+ Yrs Japanese: Speak English 'very well' 4 B16001_071 B16001 Language Spoken at Home by Ability to Speak English for the Population 5+ Yrs Japanese: Speak English less than 'very well'

The Census Bureau has two “Japanese” tables: the first relates to race and the second to language. For simplicity, let’s focus on race (B02006). The “_009” at the end indicates the column of the table; each column tabulates a different Asian nationality.

Step 2: Get the Data

There are a few ways to get the data from that table into R. One way is to use the function ?acs.fetch in the acs package. If your end result is to map the data with the choroplethr package, however, you might find it easier to use the function ?get_acs_data in the choroplethr package:

> library(choroplethr) > l = get_acs_data("B02006", "county", column_idx=9)

What’s returned is a list with 2 elements. The first element is a data frame with the (region, value) pairs. The second element is the title of the column:

str(l) List of 2 \$ df :'data.frame': 3143 obs. of 2 variables: ..\$ region: num [1:3143] 1001 1003 1005 1007 1009 ... ..\$ value : num [1:3143] 10 25 0 0 0 0 0 103 2 19 ... \$ title: chr "Asian Alone By Selected Groups: Japanese" Step 3: Analyze the Data

The first way to analyze the data is to simply look at the data frame:

> df = l[[1]] > head(df) region value 1 1001 10 2 1003 25 3 1005 0 4 1007 0 5 1009 0 6 1011 0

People who have taken my course will recognize the regions as FIPS County Codes. We can use a boxplot to look at the distribution of values:

boxplot(df\$value)

I draw two conclusions from this chart: 1) the median is very low and 2) there are two very large outliers.

To find out the names of the outliers we need to convert the FIPS Codes to English. We can do that by merging df with the data frame ?county.regions.

> data(county.regions) > head(county.regions) region county.fips.character county.name state.name state.fips.character state.abb 1 1001 01001 autauga alabama 01 AL 36 1003 01003 baldwin alabama 01 AL 55 1005 01005 barbour alabama 01 AL 15 1007 01007 bibb alabama 01 AL 2 1009 01009 blount alabama 01 AL 16 1011 01011 bullock alabama 01 AL > df2 = merge(df, county.regions) > df2 = df2[order(-df2\$value), ] > head(df2) region value county.fips.character county.name state.name state.fips.character state.abb 548 15003 150984 15003 honolulu hawaii 15 HI 205 6037 103180 06037 los angeles california 06 CA 216 6059 33211 06059 orange california 06 CA 229 6085 28144 06085 santa clara california 06 CA 2971 53033 21493 53033 king washington 53 WA 223 6073 18592 06073 san diego california 06 CA

So the outliers are Honolulu county and Los Angeles county. San Francisco isn’t even in the top 6. So if I ever decide to give up my career in technology for a career focused on Japanese, I should move to Honolulu!

It’s also easy to create a choropleth map of the values. This allows us to see the geographic distribution of the values.

library(choroplethrMaps) county_choropleth(df, title = "2012 County Estimates:nNumber of Japanese per County")

According to this map, by living on the west coast I am already in a part of the country with a high concentration of Japanese people.

Conclusion

A final note to my Japanese friends: どう思いますか？アメリカで一番興味がある場所はホノルルとロサンゼルスですか？口コミしてください！

jQuery.post( "http://www.arilamstein.com/", { 'fca_eoi_track_form_id': 400 } );

LEARN TO MAP CENSUS DATA Subscribe and get my free email course: Mapping Census Data in R!

100% Privacy. We don’t spam.

The post How to Search for Census Data from R appeared first on AriLamstein.com.

To leave a comment for the author, please follow the link and comment on their blog: AriLamstein.com » R. R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

Categories: Methodology Blogs

### simmer 3.0.0 is on CRAN

Mon, 2015-11-16 07:00

(This article was first published on FishyOperations, and kindly contributed to R-bloggers)

I’m very pleased to announce the first CRAN release of simmer. (https://cran.rstudio.com/web/packages/simmer/). This release has been realised thanks to the efforts made by Iñaki.

To reiterate a bit, simmer is a discrete-event simulation (DES) package for R. It is the first R package that focuses on creating a robust DES framework for R. It provides a framework somewhat similar as e.g SimPy and SimJulia but is arguably written from a bit more of an applied angle.

R might not be the most efficient language to implement a DES framework due to its method of memory allocation. Therefore, simmer implements a C++ backend by making use of Rcpp.

This release follows the philosophy and workflow of the pre-release version but adds a more robust event-based C++ backend and a more flexible frontend. These are the main improvements:

• Enhanced programmability The timeout activity is more than just a delay. It admits a user-defined function, which can be as complex as needed in order to interact with the simulation model. The old v2.0 was no more than a queueing network simulator. This feature makes simmer a flexible and generic DES framework. Moreover, we have finally got rid of the infamous add_skip_event function to implement a more flexible and user-friendly branching method.
• Robustness The event-based core design is rigorous and simple, which makes simmer faster and less error-prone, at the same level of other state-of-the-art DES frameworks.
Much better performance. Instead of creating n arrivals beforehand, this release leverages the concept of generator of arrivals, which is faster and more flexible. At the same time, the concept of trajectory as a chain of activities is implemented entirely in C++ internally. Our tests show that simmer is even faster than SimPy when it comes to simulate queueing networks.
• Replication In the pre-release, replication was implemented inside simmer. This no longer makes sense since, with the current design, it is more than straightforward to replicate and even parallelize the execution of replicas using standard R tools.

A nice getting started manual is provided in the vignette: https://cran.r-project.org/web/packages/simmer/vignettes/introduction.html.

For feedback, questions and bug reports please create an issue at the GitHub repository: https://github.com/Bart6114/simmer. Here you can also find the latest development version.

Let’s close with a simple example:

library(simmer) coffee_prep_duration<-function(){ coffees<-list(espresso=runif(1, min=5, max=10), macchiato=runif(1, min=6, max=12), americano=runif(1, min=3, max=8)) sample(coffees, 1)[[1]] } coffee_trajectory<- create_trajectory("I want coffee!") %>% seize("barista", 1) %>% timeout(coffee_prep_duration) %>% release("barista", 1) env<- simmer("coffee shop") %>% add_resource("barista", 2) %>% add_generator("caffeine addict", coffee_trajectory, function() abs(rnorm(1, 5, 2))) %>% run(until=480) plot_evolution_arrival_times(env, type = "flow_time")

Have fun with it! (and let us know what we can improve further)

To leave a comment for the author, please follow the link and comment on their blog: FishyOperations. R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

Categories: Methodology Blogs

### Polar Circles

Mon, 2015-11-16 06:00

(This article was first published on Ripples, and kindly contributed to R-bloggers)

You cannot find peace by avoiding life (Virginia Woolf)

Combining polar coordinates, RColorBrewer palettes, ggplot2 and a simple trigonometric function to define the width of the tiles is easy to produce nice circular plots like these:

Do you want to try? Here you have the code:

library(ggplot2) library(dplyr) library(RColorBrewer) n=500 m=50 w=sapply(seq(from=-3.5*pi, to=3.5*pi, length.out=n), function(x) {abs(sin(x))}) x=c(1) for (i in 2:n) {x[i]=x[i-1]+1/2*(w[i-1]+w[i])} expand.grid(x=x, y=1:m) %>% mutate(w=rep(w, m))-> df opt=theme(legend.position="none", panel.background = element_rect(fill="white"), panel.grid=element_blank(), axis.ticks=element_blank(), axis.title=element_blank(), axis.text=element_blank()) ggplot(df, aes(x=x,y=y))+geom_tile(aes(fill=x, width=w))+ scale_fill_gradient(low=brewer.pal(9, "Greens")[1], high=brewer.pal(9, "Greens")[9])+ coord_polar(start = runif(1, min = 0, max = 2*pi))+opt ggplot(df, aes(x=x,y=y))+geom_tile(aes(fill=w, width=w))+ scale_fill_gradient(low=brewer.pal(9, "Reds")[1], high=brewer.pal(9, "Reds")[9])+ coord_polar(start = runif(1, min = 0, max = 2*pi))+opt ggplot(df, aes(x=x,y=y))+geom_tile(aes(fill=y, width=w))+ scale_fill_gradient(low=brewer.pal(9, "Purples")[1], high=brewer.pal(9, "Purples")[9])+ coord_polar(start = runif(1, min = 0, max = 2*pi))+opt ggplot(df, aes(x=x,y=y))+geom_tile(aes(fill=w*y, width=w))+ scale_fill_gradient(low=brewer.pal(9, "Blues")[9], high=brewer.pal(9, "Blues")[1])+ coord_polar(start = runif(1, min = 0, max = 2*pi))+opt

To leave a comment for the author, please follow the link and comment on their blog: Ripples. R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

Categories: Methodology Blogs

### Why Submitting Ideas To The R Consortium Is A Good Idea

Mon, 2015-11-16 05:16

(This article was first published on Mango Solutions, and kindly contributed to R-bloggers)

By Steph Locke

Announced in June 2015, the R Consortium aims to support the R ecosystem and community and help grow the adoption of R. To deliver on these aims, the R Consortium has an Infrastructure Steering Committee (ISC). The ISC are responsible for directing technical focus and overseeing projects to deliver improvements.

You’ll note in all this, that there’s no mention of a central vision for R beyond supporting it. That’s because the R Consortium is looking to support the projects that the community thinks will help it grow and thrive.

Every six months, the R Consortium will grant awards to proposals for projects that they feel best help the community. The ISC are responsible for receiving, evaluating, and selecting projects.

The first proposal selected by the ISC is R-Hub. R-Hub will provide a vital pre-CRAN build and check process, designed to give package developers a facility to test their proposed builds before it goes to CRAN, making the process quicker for developers and easier for the CRAN volunteers.

The next proposals to be accepted, will need to be submitted by January 10, 2016 and the announcements will be made mid-February. On the R Consortium site, the ISC have provided guidance and information on the process.

To help get you started, Mango’s Steph Locke has built a template proposal that you can customise for your own proposal. This template is intended to be extensive but will need customising, whether that’s reducing the number of sections to reflect a very modest proposal, or to add specific topical sections.

At Mango, we actively encourage our consultants to contribute to the community and they’re working on a number of proposals for consideration; one of the most recent is the satRday – user-group driven conferences, facilitated by the R Consortium.

We strongly encourage all our readers to think about the things they wish could be better and come up with some ideas about how they could help solve the improve the situation. We will help raise awareness and give feedback on any proposals pre-submission so if you have an idea, let us know by either tweeting us directly @MangoTheCat or using the hash-tag #ISCProposal!

To leave a comment for the author, please follow the link and comment on their blog: Mango Solutions. R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

Categories: Methodology Blogs

### Aggregate – A Powerful Tool for Data Frame in R

Mon, 2015-11-16 02:47

(This article was first published on DataScience+, and kindly contributed to R-bloggers)

This post gives a short review of the aggregate function as used for data.frames and presents some interesting uses: from the trivial but handy to the most complicated problems I have solved with aggregate.

Aggregate (data.frame): Technical Overview

Aggregate is a function in base R which can, as the name suggests, aggregate the inputted data.frame d.f by applying a function specified by the FUN parameter to each column of sub-data.frames defined by the by input parameter.

The by parameter has to be a list. However, since data.frame‘s are handled as (named) lists of columns, one or more columns of a data.frame can also be passed as the by parameter. Interestingly, if these columns are of the same data.frame as the one inputted as x, those columns are not passed on to the FUN function.

The function to apply has to be able to accept a vector (since it will be called with parts of a column of a data.frame as input).

The sub-data.frames defined by the by input parameter can be thought of as logical indexing:

d.f <- data.frame(rating = c("AAA", "A", "A", "AAA", "BB", "BB", "AAA", "A")) i <- 1 by <- d.f\$rating sub.data.frame <- d.f[by == unique(by)[i], ]

and do this for every i between 1 and length(unique(by)). Note that the by variable doesn’t have to agree with one (or more) column of the data.frame but could be anything. Hence, one can reproduce the aggregate functionality by a for cycle running the cycle variable over the unique values of the variable passed as by and an sapply applying the function passed as FUN to each column of the data.frame sub.data.frame. Such a workaround however would be very difficult to document, as it would be unclear what (and why) this code is actually doing.

Aggregate always returns a data.frame as a result. This data.frame will contain the (now unique) values from the input parameter by as the first column and then columns containing the results of the call to the function in the FUN parameter applied to the parts of the columns of the inputted data.frame. It is interesting to note that if the function FUN returns multiple values, the class of the columns of the result data.frame will be list or something a list can be cast to (see the last example below).

It is important to note that the function call is applied to nameless vectors rather than named columns of a data.frame and hence referring to the names of the data.frame will not work, nor will column references such as s.d.f[,1].

Basic Examples

The most basic uses of aggregate involve base functions such as mean and sd. It is indeed one of the most common uses of aggregate to compare the mean or other properties of sample groups.

Recently I reproduced calculations from an Excel sheet. Most formulae were subtotals and grand totals. The Excel sheet was not very comfortably organized for this purpose: sums over rows, columns and totals of those sums were used. In R, I have changed the data to a star schema representation (when all metadata are represented row-wise and every value gets its own row) using reshape2 package and melt then used aggregate along different variables to get the different totals. The less variables you use in by the more aggregated the end-result: the grand total along a dimension is simply using that dimension as “by”, while subtotals can be achieved using multiple variables as by. The FUN in this case was of course sum.

One handy use of aggregate and a base function is getting the number of appearances of the various values:

values <- data.frame(value = c("a", "a", "a", "a", "a", "b", "b", "b", "c", "c", "c", "c")) nr.of.appearances <- aggregate(x = values, by = list(unique.values = values\$value), FUN = length)

My favourite use of aggregate with a base function is getting the last day of each month in a series of dates. To do so, one can use the following code (assuming your dates are stored in a “YYYY-MM-DD” format as strings or as Date):

dates <- data.frame(date = as.Date("2001-01-01", format = "%Y-%m-%d") + 0:729) dates date 1 2001-01-01 2 2001-01-02 3 2001-01-03 4 2001-01-04 ..... last.day <- aggregate(x = dates["date"], by = list(month = substr(dates\$date, 1, 7)), FUN = max) last.day month date 1 2001-01 2001-01-31 2 2001-02 2001-02-28 3 2001-03 2001-03-31 4 2001-04 2001-04-30 .....

This came in very handy when working with banking information where the last day of the month depended on banking holidays as well as weekends.

More advanced uses of aggregate depend on writing your own function, e.g. anonymous functions passed on as the FUN parameter. To do so, one can use the syntax

# do not run the syntax aggregate(x = d.f, by = by.list, FUN = function(s.d.f){y <- s.d.f; return(y)}

The possible uses range from calling complex portfolio risk metrics for the homogeneous risk groups of a portfolio via fitting a distribution to categories of samples to anything you can image, really.

Here is an example with a “complex” portfolio risk metric (exposure to different counterparties in different asset classes):

assets <- data.frame(asset.class = c("equity", "equity","equity", "option","option","option", "bond", "bond"), rating = c("AAA", "A", "A", "AAA", "BB", "BB", "AAA", "A"), counterparty.a = c(runif(3), rnorm(5)), counterparty.b = c(runif(3), rnorm(5)), counterparty.c = c(runif(3), rnorm(5))) assets asset.class rating counterparty.a counterparty.b counterparty.c 1 equity AAA 0.9026004 0.6029417 0.8629453 2 equity A 0.8834034 0.5809589 0.4654721 3 equity A 0.1007586 0.9368537 0.3090811 4 option AAA -1.0508915 0.7171532 0.2224984 .....

Here is the use of aggregate() function.

exposures <- aggregate(x = assets, by = assets, FUN = function(market.values){ sum(pmax(market.values, 0)) }) exposures asset.class rating counterparty.a counterparty.b counterparty.c 1 bond A 1.0038714 0.6382029 2.2822936 2 equity A 0.9841620 1.5178126 0.7745532 3 bond AAA 0.0000000 0.0000000 0.0000000 4 equity AAA 0.9026004 0.6029417 0.8629453 .....

Next up: fitting a Gaussian distribution to observations by categories:

library(MASS) categories <- data.frame(category = c("a", "a", "a", "a", "a", "b", "b", "b", "b", "b", "c", "c", "c", "c")) observations <- data.frame(observation = c(rnorm(5, mean = 3, sd = 0.2), rnorm(5, mean = -2, sd = 0.4), rnorm(4, mean = 0, sd = 1)))

Below we use the aggregate() function to find the mean and standard deviation by categories.

distr.estimate <- aggregate(x = observations, by = categories, FUN = function(observations){ fitdistr(observations, densfun = "normal")\$estimate }) distr.estimate category observation.mean observation.sd 1 a 3.0606926 0.1779962 2 b -2.1446040 0.1658481 3 c -0.1881841 0.5613013

This last example showcases several interesting properties. First, the data.frame to aggregate and the list of by variables don’t have to be the same. While this is implied in other places of the post, this is an explicit example of such a setup. Secondly, the function passed as FUN is not only an anonymous function, it is curried from a function with more than one input parameter. A function of a single input variable observations has been created from the two-input variable function fitdistr: fixing one of the input variables by setting densfun = "normal". Thirdly, rather than returning the full return value of the fitdistr function, the return value is restricted to the element estimate from the return value. And last but not least, the return value of the anonymous function passed to FUN consists of two variables and not only one. Interestingly, aggregate casts the return value from list to a matrix and names the elements for us. However, these names can’t be used to reference the columns of the matrix. You can however reference them as follows:

distr.estimate\$observation[1,][["mean"]] [1] 3.016988 Closing Words

I hope that you have found the above useful. Now that you are more familiar with aggregate, it is time for the truth: everything above and much more can be done with data.table, and with a much faster performance. However, data.table has a complex syntax and one really has to understand how things work under the hood, while aggregate is simple and insightful. Until you are comfortable with both the logic of aggregation and the syntax of data.table, it is a worthy investment to first write the code using aggregate and then optimize it by rewriting it using data.table.

For those of you who are interested, a dedicated post is coming where the above is redone with data.table, along with some additional use cases specific to data.table.

To leave a comment for the author, please follow the link and comment on their blog: DataScience+. R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

Categories: Methodology Blogs

### Interactive Data Science with R in Apache Zeppelin Notebook

Mon, 2015-11-16 02:42

(This article was first published on SparkIQ Labs Blog » R, and kindly contributed to R-bloggers)

Introduction

The objective of this blog post is to help you get started with Apache Zeppelin notebook for your R data science requirements. Zeppelin is a web-based notebook that enables interactive data analytics. You can make beautiful data-driven, interactive and collaborative documents with Scala(with Apache Spark), Python(with Apache Spark), SparkSQL, Hive, Markdown, Shell and more.

However, the latest official release, version 0.5.0, does not yet support the R programming language. Fortunately NFLabs, the company driving this open source project, pointed me this pull request that provides an R Interpreter. An Interpreter is a plug-in which enables zeppelin users to use a specific language/data-processing-backend. For example to use scala code in Zeppelin, you need a spark interpreter. So, if you are impatient like I am for R-integration into Zeppelin, this tutorial will show you how to setup Zeppelin for use with R by building from source.

Prerequisites
• We will launch Zeppelin through Bash shell on Linux. If you are using Windows OS I recommend that you install and use the Cygwin terminal (It provides functionality similar to a Linux distribution on Windows).
• Make sure Java 1.7 and Maven 3.2.x are installed on your host machine and their environment variables are set.
Build Zeppelin from Source

In my case I have downloaded and unzipped the folder onto my Desktop

Step 2: Build Zeppelin

Run the following code in your terminal to build zeppelin on your host machine in local mode. If you are installing on a cluster then add these options found in the Zeppelin documentation.

\$ cd Desktop/Apache/incubator-zeppelin-rinterpreter

\$ mvn clean package -DskipTests

This will take around 6 minutes to build zeppelin, Spark and all interpreters including R, Markdown, Hive, Shell, and others. (as shown in the image below).

Step 3: Start Zeppelin

Run the following command to start zeppelin.

\$ ./bin/zeppelin-daemon.sh start

Go to localhost on your web browser and listen on port 8080. (i.e. http://localhost:8080). At this point you are ready to start creating interactive notebooks with code and graphs in Zeppelin.

Interactive Data Science

Step 1: Create a Notebook

Click the dropdown arrow next to the “Notebook” page and click “Create new note”.

Give your notebook a name or you can use the assigned default name. I named mine “Base R in Apache Zeppelin”.

To use R, use the “%spark.r” or “%spark.knitr” tags as shown in the images below. First let’s use markdown to write some instruction text.

Now let’s install some packages that we may need for our analysis.

Now let’s read in our data set. We shall use the “flights” dataset which shows flights departing New York in 2013.

Now let’s do some data manipulation using dplyr (with the pipe operator)

You can also use bar graphs and pie charts to visualize some descriptive statistics from your data.

Now let’s do some data exploration with ggplot2

Now let’s do some statistical machine learning using the caret package.

Final Remarks

Zeppelin allows you to create interactive documents with beautiful graphs using multiple programming languages. The objective of this post was to help you setup Zeppelin for use with the R programming language. Hopefully the Project Management Committees (PMC) of this wonderful open source project can release the next version with an R interpreter. It will surely make it easier to launch Zeppelin faster without having to build from source.

Also it’s worth mentioning that there is another R interpreter for Zeppelin produced by the folks at Data Layer. You can find instructions on how to use it here: https://github.com/datalayer/zeppelin-R-interpreter.

Try out both interpreters and share your experiences in the comments section below.

As a follow-up to this post, We shall see how to use Apache Spark (especially SparkR) within Zeppelin in the next blog post.

Filed under: Apache Spark, Data Science, Machine Learning, R, SparkR, Zeppelin Tagged: Data Science, R, Zeppelin

To leave a comment for the author, please follow the link and comment on their blog: SparkIQ Labs Blog » R. R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

Categories: Methodology Blogs

Sun, 2015-11-15 23:39

(This article was first published on Quintuitive » R, and kindly contributed to R-bloggers)

Markets are very smart in absorbing and reflecting information. If you think otherwise, try making money by trading. If you are new to it, make sure you don’t bet the house.

In other words, markets are efficient. At least most of the time. So then why people trade? The general believe is that there are windows during which prices of certain assets are inefficient. Thus, there are opportunities to make money. Is the presence of autocorrelation one such opportunity? Let’s find out.

To keep things simple, we will use the standard R acf function to compute the autocorrelation.

require(quantmod) spy=getSymbols("SPY",from="1900-01-01",auto.assign=F) # Note: We use adjusted close - it's unrealistic to expect anything from actual close in the # presence of dividends spy.rets = ROC(Ad(spy),na.pad=F) aa = acf(tail(spy.rets,500),main="ACF computed over the last 500 days") head(aa) # [1] 1.000000000 0.051116484 -0.037469705 -0.010014871 -0.126667484 -0.004113005

This will produce the following chart:

It shows the autocorrelation coefficients at different lags. The first lag is the correlation of the series with itself (lag 0), and, it’s always 1. The second value (0.051116484) is the correlation of the series with the series lagged by one.

The two dashed lines are the confidence intervals for the lags. They are of special interest since we are going to use them to decide when there is significant autocorrelation. How are they computed? To find the answer, I had to look at the acf’s code:

# xx is the series, conf.level is the confidence level - think 0.95 for instance conf = qnorm((1 + conf.level)/2)/sqrt(sum(!is.na(xx)))

The above is nothing else but computing confidence intervals for normal (0,1) distribution. It still puzzles me (I couldn’t find an answer quickly when I thought about it) why the correlation coefficients at different lags are distributed normally in (0,1), but that’s irrelevant.

The last question is how to trade extreme autocorrelation – do we bet that the autocorrelation persists, or do we bet that the autocorrelation fades? There are two variables here: the sign of the correlation and the sign of the last day return. A table is helpful.

Correlation Sign Return Sign Persisting Signal Fading Signal 1 1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 -1 1 -1

My gut tells me to go with the fading – markets are efficient, especially this one. So, let’s run a quick backtest, putting all this mumbo-jumbo into code.

high.acf = function(xx,conf.level=0.95,lag=1) { aa = acf(xx,plot=F) conf = qnorm((1 + conf.level)/2)/sqrt(sum(!is.na(xx))) if(abs(aa\$acf[lag+1,1,1]) > conf) sign(aa\$acf[lag+1,1,1]) else 0 } backtest.acf = function(rets, n=21, conf.level=0.95, lag=1, fade=F, dates="2004/2013") { aa = na.trim(rollapplyr(rets, width=n+lag, FUN=high.acf, conf.level=conf.level, lag=lag)) bb = merge(rets, aa, all=F) ind = sign(bb[,1]*bb[,2]) if(fade) ind = -ind cc = merge(rets, lag.xts(ind), all=F) dd = cc[,1]*cc[,2] strat = dd[dates] n.win = NROW(which(as.numeric(strat) > 0, arr.ind=T)) n.trades = NROW(which(as.numeric(strat) != 0, arr.ind=T)) str = paste(round(n.win/n.trades*100,2), "% [", n.win, "/", n.trades, "]", sep="") print(str) } # About 3 (trading) months of history backtest.acf(spy.rets, dates="2004/2013", fade=T, n=63) # [1] "54.88% [45/82]"

The percentage looks ok, but the sample is small (only 82 opportunities over 10 years). Is it possible that there is some opportunity here – too early to tell. Certainly worth looking further though, IMO.

The post Trading Autocorrelation? appeared first on Quintuitive.

To leave a comment for the author, please follow the link and comment on their blog: Quintuitive » R. R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

Categories: Methodology Blogs

### Rcpp 0.12.2: More refinements

Sun, 2015-11-15 22:43

(This article was first published on Thinking inside the box , and kindly contributed to R-bloggers)

The second update in the 0.12.* series of Rcpp is now on the CRAN network for GNU R. As usual, I will also push a Debian package. This follows the 0.12.0 release from late July which started to add some serious new features, and builds upon the 0.12.1 release in September. It also marks the sixth release this year where we managed to keep a steady bi-montly release frequency.

Rcpp has become the most popular way of enhancing GNU R with C or C++ code. As of today, 512 packages on CRAN depend on Rcpp for making analytical code go faster and further. That is up by more than fifty package from the last release in September (and we recently blogged about crossing 500 dependents).

This release once again features pull requests from two new contributors with Nathan Russell and Tianqi Chen joining in. As shown below, other recent contributors (such as such as Dan) are keeping at it too. Keep’em coming! Luke Tierney also email about a code smell he spotted and which we took care of. A big Thank You! to everybody helping with code, bug reports or documentation. See below for a detailed list of changes extracted from the NEWS file.

Changes in Rcpp version 0.12.2 (2015-11-14)
• Changes in Rcpp API:

• Correct return type in product of matrix dimensions (PR #374 by Florian)

• Before creating a single String object from a SEXP, ensure that it is from a vector of length one (PR #376 by Dirk, fixing #375).

• No longer use STRING_ELT as a left-hand side, thanks to a heads-up by Luke Tierney (PR #378 by Dirk, fixing #377).

• Rcpp Module objects are now checked more carefully (PR #381 by Tianqi, fixing #380)

• An overflow in Matrix column indexing was corrected (PR #390 by Qiang, fixing a bug reported by Allessandro on the list)

• Nullable types can now be assigned R_NilValue in function signatures. (PR #395 by Dan, fixing issue #394)

• operator<<() now always shows decimal points (PR #396 by Dan)

• Matrix classes now have a transpose() function (PR #397 by Dirk fixing #383)

• operator<<() for complex types was added (PRs #398 by Qiang and #399 by Dirk, fixing #187)

• Changes in Rcpp Attributes:

• Enable export of C++ interface for functions that return void.

• Changes in Rcpp Sugar:

• Added new Sugar function cummin(), cummax(), cumprod() (PR #389 by Nathan Russell fixing #388)

• Enabled sugar math operations for subsets; e.g. x[y] + x[z]. (PR #393 by Kevin and Qiang, implementing #392)

• Changes in Rcpp Documentation:

• The NEWS file now links to GitHub issue tickets and pull requests.

• The Rcpp.bib file with bibliographic references was updated.

Thanks to CRANberries, you can also look at a diff to the previous release As always, even fuller details are on the Rcpp Changelog page and the Rcpp page which also leads to the downloads page, the browseable doxygen docs and zip files of doxygen output for the standard formats. A local directory has source and documentation too. Questions, comments etc should go to the rcpp-devel mailing list off the R-Forge page.

This post by Dirk Eddelbuettel originated on his Thinking inside the box blog. Please report excessive re-aggregation in third-party for-profit settings.

To leave a comment for the author, please follow the link and comment on their blog: Thinking inside the box . R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

Categories: Methodology Blogs

### Partools, Recommender Systems and More

Sun, 2015-11-15 18:07

Recently I attended a talk by Stanford’s Art Owen, presenting work done with his student, Katelyn Gao. This talk touched on a number of my interests, both mathematical and computational. What particularly struck me was that Art and Katelyn are applying a very old — many would say very boring — method to a very modern, trendy application: recommender systems.  (See also a recent paper by P.O. Perry.)

The typical example of a recommender system involves movie review data, such as that in Netflix and MovieLens, with various users rating various movies. We have a matrix, with rows representing users and columns for movies, but it is a very sparse matrix, because most users don’t rate most movies. So, we’d like to guess the missing elements, a process sometimes termed the matrix completion problem.

There are many methods for this, such as matrix factorization, as exemplified in the R package NMF. I much prefer the approach taken by Art and Katelyn, as it is more statistical and thus potentially can provide more insight into the underlying processes.

I was inspired by their research to start a project that takes a new approach to random effects models in general. I’ll explain that later, but let’s first address the major computational problems that can occur with random effects models when applied to big data. Note carefully that not only is this a problem with long run times, but also with memory limitations. For example, Pennell and Dunson wrote, in Fitting Semiparametric Random Effects Models to Large Data Sets, Biostatistics, 2007, “SAS PROC MIXED ran out of memory when we attempted to fit a model with random smoking effects.”

One of the features of my partools package is what I call Software Alchemy, with the term alluding to the fact that the procedure turns a statistical algorithm into a statistically equivalent version that is easily parallelized. The idea is simple: Break i.i.d. data into chunks, apply the original algorithm to each chunk, then average the results. It’s easy to prove that the statistical accuracy of this method matches that of the original one, and one not only gets a speedup but also may solve the memory problem, since each chunk is smaller. See my paper on this, to appear soon in the Journal of Statistical Software, but also outlined in my book on parallel computation in data science.  (Chunking has been tried in random effects contexts before; see references in the Perry paper above.)

Here is an example of all this, applied to the MovieLens data, 1-million records version. I’ll use Art and Katelyn’s model here, which writes the rating by user i of movie j as

Yij = μ + ai + bj + εij

where μ is an unknown constant and ai + bj + εij are independent, mean-0 variables (not necessarily normal) with variances σa2, σb2 and σε2. One can also add covariates, but we’ll keep things simple here. In addition to estimating the variances, one can estimate the ai and bj, and thus predict the missing Yij. (See the b component in the output of lmer(), to be discussed below. It lists the estimated ai and bj, in that order.)

Art and Katelyn (as well as Perry) use the Method of Moments rather than the standard methods for random effects models, so that (a) they reduce computation, by deriving explicit solutions and (b) are not tied to assumptions of normality. I use MM in my project too, but here I’ll use the famous R package for random effects models, lme4.  I ran the experiments here on a quad-core machine, using an R parallel cluster of that size.

I first ran without chunking:
> system.time(lmerout <- lmer(V3 ~ (1|V1) +  (1|V2),ml1m))
user system elapsed
1105.129 3.692 1108.946
> getME(lmerout,'theta')
V1.(Intercept) V2.(Intercept)
0.3994240 0.6732281

Not good — almost 20 minutes! So, I turned to partools.The main Software Alchemy function is cabase(). For a given statistical algorithm, one needs to specify the function that implements the algorithm and a function to extract the estimated parameters from the output of the algorithm. There is already a wrapper for lm() in the package, so I just modified it for lmer():

calmer <- function (cls, lmerargs) { ovf <- function(u) { tmp <- paste("lmer(", lmerargs, ")", collapse = "") docmd(tmp) } estf <- function(lmerout) { getME(lmerout,'theta') } cabase(cls, ovf, estf) }

I had earlier used partools function distribsplit() to distribute my data frame m1m to the cluster nodes, in chunks of the same name. So, the call to cabase() does the estimation on each chunk, collects the results, and averages them. Here is what happened:

> system.time(print(calmer(cls,'V3 ~ (1|V1) + (1|V2),ml1m'))) \$thts V2.(Intercept) V1.(Intercept) [1,] 0.6935344 0.3921150 [2,] 0.6623718 0.3979066 [3,] 0.6505459 0.3967739 [4,] 0.6825985 0.4103681 \$tht V2.(Intercept) V1.(Intercept) 0.6722626 0.3992909 user system elapsed 0.812 0.080 189.744

We got an almost 6-fold improvement! And note that we used only 4 cores; this is the superlinear effect that I discuss in my Software Alchemy paper. Just as important, we got essentially the same estimates as from the non-chunked computation. (Please note: The chunked version can be shown to have the same asymptotic variance as the original one; neither is an “approximation” to the other.)

But wait a minute. Didn’t I say that Software Alchemy is intended for i.i.d. data? I did. In the paper I said that the theory could easily be extended to the case of independent but non identically-distributed data, but it would be difficult to state general conditions for this.

Random effects models generally do NOT assume the data are identically distributed. In Art and Katelyn’s paper, for instance, they define indicator variables zij representing that user i rated movie j, and treat these as constants. So, for example, the row counts — numbers of movies rated by each user — are treated as constants.

It’s my contention that such quantities should be treated as random, just like other user behavior. (We haven’t brought in variables such as movie genre here, but if we did, a similar statement would hold.) And once we do that, then voila!, we now have i.i.d. data, and Software Alchemy applies.

In my paper on arXiv,  I also show why it doesn’t matter: The MM estimators are going to turn out essentially indentical, regardless of assuming fixed or random quantities of this type.

Moreover, the advantage of doing things this way goes beyond simply being able to use Software Alchemy.  Those quantities may be of statistical interest in their own right. For instance, consider a study of children within families. In the classical approach, the number of kids in a family in a random effects model would be treated as fixed. But research has shown that family size is negatively correlated with income, and the number of kids should be an important variable — a random variable — in the analysis.

Now, back to MM. The analysis will compute a lot of quantities involving within-user variances and so on. If so, it will help to group together rows or columns of our input data frame, e.g. group by user and group by movie. How might we do this using par tools?

The R function split() comes to mind, but there is an issue of combining the outputs of split(),  which are R lists, from the various cluster nodes. Note in particular that some user, for example, may have data only at one of the nodes. But this is exactly the situation that the partools function addlists() is designed to handle!

That helps, but still we have to worry about overhead incurred by shipping so much data back and forth between the worker nodes and manager node in the cluster. Also, what about memory? Again, we should try (may need) to avoid holding the entire  data set in memory at one time.  The  partools function filesplit() helps in that regard, but this continues to be a computational challenge. It may be that the only good solutions involve C.

Categories: Methodology Blogs

### Using htmlwidgets with knitr and Jekyll

Sun, 2015-11-15 12:39

(This article was first published on Brendan Rocks >> R, and kindly contributed to R-bloggers)

A few weeks ago I gave a talk at BARUG (and wrote a post) about blogging with the excellent knitr-jekyll repo. Yihui’s system is fantastic, but it does have one drawback: None of those fancy new htmlwidgets packages seem to work…

A few people have run into this. I recently figured out how to fix it for this blog (which required a bit of time reading through the rmarkdown source), so I thought I’d write it up in case it helps anyone else, or my future-self.

TL;DR

You can add a line to build.R which calls a small wrapper function I cobbled together (brocks::htmlwidgts_deps), add a snippet of liquid syntax to ./_layouts/post.html, and you’re away.

What’s going on?

Often, when you ‘knit’ an .Rmd file to html, (perhaps without knowing it) you’re doing it via the rmarkdown package, which adds its own invisible magic to the process. Behind the scenes, rmarkdown uses knitr to convert the file to markdown format, and then uses pandoc to convert the markdown to HTML.

While knitr executes R code and embeds results, htmlwidgets packages (such as leaflet, DiagrammR, threejs, and metricsgraphics) also have js and css dependencies. These are handled by rmarkdown’s second step, and so don’t get included when using knitr alone.

The rmarkdown invisible magic works as follows:

• It parses the .Rmd for special dependencies objects, linking to the js/css source (by calling knitr::knit_meta)
• It then (by default) combines their source-code into a huge data:uri blob, which it writes to a temp-file
• This is injected into the the final HTML file, by passing it to pandoc’s --include-in-header argument
A fix: htmlwdigets_deps

Happily, including bits of HTML in other bits of HTML is one of Jekyll’s strengths, and it’s possible to high-jack the internals of rmarkdown to do something appropriate. I did this with a little function htmlwdigets_deps, which:

• Copies the js/css dependencies from the R packages, into a dedicated assets folder within in your blog

• Writes a little HTML file, containing the links to the source code above

With a small tweak to the post.html file, Jekyll’s liquid templating system can be used to pull in that little HTML file, if htmlwidgets are detected in your post.

If you’re using knitr-jekyll, all that’s needed to make everything work as you’d expect, is to call the function from your build.R file, like so:

local({ # Your existing configurations... # See https://github.com/yihui/knitr-jekyll/blob/gh-pages/build.R brocks::htmlwidgets_deps(a) })

(The parameter a refers to the input file — if you’re using a build file anything like Yihui’s example, this will work fine.)

If you’d like to have a look at the internals of htmlwidgets_deps yourself, it’s in my personal package up on GitHub. Long story short, it hi-jacks rmarkdown:::html_dependencies_as_string. The rest of this post walks through what it actually does.

1. Copying dependencies to your site

To keep things transparent, the dependency source files are kept in their own folder (./htmlwidgets_deps). If it doesn’t exist, it’ll be created. This behaviour is different to the rmarkdown default of compressing everything into huge in-line data:uri blobs. While that works great for keeping everything in one big self-contained file (e.g. to email to someone), it makes for a very slow web page. For a blog, having separate files is preferable, as it allows the browser to load files asynchronously, reducing the load time.

After compiling your sites, if you’ve used htmlwidgets you’ll have an extra directory within your blog, containing the source for all the dependencies, a bit like this:

- _includes - _layouts - _posts - _sass - _site - _source - js/ - css/ - htmlwidgets_deps/ - d3-3.5.3/ - LICENCE - bower.json - d3.js - d3.min.js - jquery-1.11.1 - AUTHORS.txt - jquery.min.js - ... - ... 2. Writing the extra HTML

Once all the dependencies are ready to be served from your site, you still need to add HTML pointers to your blog post, so that it knows where to find them. htmlwidgets_deps automates this, by adding a file for each htmlwidgets post to the ./_includes directory (which is where Jekyll goes to look for HTML files to include). For each post which requires it, the extra HTML file will be generated in the htmlwidgets sub-directory, like this:

- _includes/ - htmlwidgets/ - my-new-htmlwidgets-post.html - footer.html - head.html - header.html - _layouts/ ...

The file itself if pretty simple. Here’s an example:

<script src="{{ "/htmlwidgets_deps/htmlwidgets-0.5/htmlwidgets.js" | prepend: site.baseurl }}"></script> <script src="{{ "/htmlwidgets_deps/jquery-2.1.3/dist/jquery.min.js" | prepend: site.baseurl }}"></script> <script src="{{ "/htmlwidgets_deps/d3-3.5.3/d3.min.js" | prepend: site.baseurl }}"></script> <link href="{{ "/htmlwidgets_deps/metrics-graphics-2.1.0/dist/metricsgraphics.css" | prepend: site.baseurl }}" rel="stylesheet" /> <script src="{{ "/htmlwidgets_deps/metrics-graphics-2.1.0/dist/metricsgraphics.min.js" | prepend: site.baseurl }}"></script> <script src="{{ "/htmlwidgets_deps/metricsgraphics-binding-0.8.5/metricsgraphics.js" | prepend: site.baseurl }}"></script>

The HTML comes pre-wrapped in the usual liquid syntax.

3. Including the extra HTML

Now you have a little file to include, you just need to get it into the HTML of the blog post. Jekyll’s templating system liquid is all about doing this.

Because htmlwdigets_deps gives the dependency file the same name as your .Rmd input (and thus the post), it’s quite easy to write a short {% include %} statement, based on the name of the page itself. However, things get tricky if the file doesn’t exist. By default, htmlwdigets_deps only produces files when necessary (e.g. when you are actually using htmlwidgets). To handle this, I used a plugin, providing the file_exists function.

Adding the following the bottom of ./_layouts/default.html did the trick. You could also use ./_layouts/post.html if you wanted to. It’s a good idea to put it towards to the bottom, otherwise the page won’t load until all the htmlwdigets dependencies are loaded, which could make things feel rather slow.

<!-- htmlwidgets dependencies --> {% assign dep_file = page.url | replace_regex:'/\$','.html' | prepend : 'htmlwidgets' %} {% assign dep_file_inc = dep_file | prepend : '_includes/' %} {% capture hw_used %}{% file_exists {{ dep_file_inc }} %}{% endcapture %} {% if hw_used == "true" %} {% include {{dep_file}} %} {% endif %} With GitHub Pages

The solution above proves a little tricky if you’re using GitHub pages, as this doesn’t allow plugins. While I’m sure an expert with the liquid templating engine could come up with a brilliant solution to this, in lieu, I present a filthy untested hack.

By setting the htmlwdigets_deps parameter always = TRUE, a dependencies file will always be produced, even if there’s no htmlwidgets detected (the file will be empty). This means that you can do-away with the logic part (and the plugin), and simply add the lines:

<!-- htmlwidgets dependencies --> {% assign dep_file = page.url | replace_regex:'/\$','.html' | prepend : 'htmlwidgets' %} {% include {{dep_file}} %}

The disadvantage is that you’ll end up with some empty HTML files in ./_includes/htmlwidgets/, which may or may not bother you. If you’re only going to be using htmlwidgets for blog posts (and not the rest of your site) I’d recommend doing this for the ./_layouts/post.html file, (as opposed to default.html) so that other pages don’t have trouble finding dependencies they don’t need.

If you give this a crack, let me know!

How to do the same

In summary:

• Add the snippet of liquid syntax to one of your layout files

• Add the following line to your build.R file, just below the call to knitr::knit

brocks::htmlwidgets_deps(a)

And you should be done!

Showing Off

After all that, it would be a shame not to show off some interactive visualisations. Here are some of the htmlwidgets packages I’ve had the chance to muck about with so far.

MetricsGraphics

MetricsGraphics.js is a JavaScript API, built on top of d3.js, which allows you to produce a lot of common plots very quickly (without having to start from scratch each time). There’s a few libs like this, but MetricsGraphics is especially pleasing. Huge thanks to Ali Almossawi and Mozilla, and also to Bob Rudis for the R interface.

library(metricsgraphics) plots <- lapply(1:4, function(x) { mjs_plot(rbeta(1000, x, x), width = 300, height = 300, linked = TRUE) %>% mjs_histogram(bar_margin = 2) %>% mjs_labs(x_label = sprintf("Plot %d", x)) }) mjs_grid(plots)

leaflet

leaflet.js allows you to create beautiful, mobile-friendly maps (based on OpenStreetMap data), incredibly easily. Hat tip to Vladimir Agafonkin, and Joe Cheng et al for the R interface!

Here’s the Pride of Spitalfields, which I occasionally pine for, from beneath the palm trees of sunny California.

library(leaflet) m <- leaflet() %>% addTiles() %>% # Add default OpenStreetMap map tiles addMarkers(lng = -0.07125, lat = 51.51895, popup = "Reasonably Priced Stella Artois") m

threejs

three.js is a gobsmackingly brilliant library for creating animated, interactive 3D graphics from within a Web browser. Here’s an interactive 3D globe with the world’s populations mapped as, erm, light-sabers. Probably not as informative as a base graphics plot, but it is much more Bond villianish. Drag it around and have a zoom!

library("threejs") library("maps") ## ## # ATTENTION: maps v3.0 has an updated 'world' map. # ## # Many country borders and names have changed since 1990. # ## # Type '?world' or 'news(package="maps")'. See README_v3. # data(world.cities, package = "maps") cities <- world.cities[order(world.cities\$pop,decreasing = TRUE)[1:1000],] value <- 100 * cities\$pop / max(cities\$pop) # Set up a data color map and plot col <- rainbow(10, start = 2.8 / 6, end = 3.4 / 6) col <- col[floor(length(col) * (100 - value) / 100) + 1] globejs(lat = cities\$lat, long = cities\$long, value = value, color = col, atmosphere = TRUE)

Kudos to Ricardo Cabello/mrdoob for three.js, and Bryan W. Lewis for the R package.

Wrapping up

So, there we go. I hope this might be useful to someone. If you do have a go at using this, let me know how you get on!

To leave a comment for the author, please follow the link and comment on their blog: Brendan Rocks >> R. R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

Categories: Methodology Blogs

### Wind in Netherlands

Sun, 2015-11-15 09:20

(This article was first published on Wiekvoet, and kindly contributed to R-bloggers)

In climate change discussions, everybody talks about temperature. But weather is much more than that. There is at least rain and wind as directly experienced quality, and air pressure as measurable quantity. In the Netherlands, some observation stations have more than a century of daily data on these things. The data may be broken in the sense that equipment and location can have changed. To quote: “These time series are inhomogeneous because of station relocations and changes in observation techniques. As a result, these series are not suitable for trend analysis. For climate change studies we refer to the homogenized series of monthly temperatures of De Bilt link or the Central Netherlands Temperature link.” Since I am not looking at temperature but wind, I will keep to this station’s data. Data Data are from daily observations from KNMI. I have chosen station De Kooy. For those less familiar with Dutch geography, this is close to Den Helder, in the tip North West of Netherlands. This means pretty close to the North Sea, Wadden Sea and Lake IJssel. Wind should be relatively unhindered there. The data themselves are daily observations. For wind there are: DDVEC     Vector mean wind direction in degrees
(360=north, 90=east, 180=south, 270=west, 0=calm/variable)
FHVEC     Vector mean windspeed (in 0.1 m/s)
FG             Daily mean windspeed (in 0.1 m/s)
FHX          Maximum hourly mean windspeed (in 0.1 m/s)
FHXH       Hourly division in which FHX was measured
FHN          Minimum hourly mean windspeed (in 0.1 m/s)
FHNH       Hourly division in which FHN was measured
FXX          Maximum wind gust (in 0.1 m/s)
FXXH       Hourly division in which FXX was measured
The header of the data downloaded contains this, and much more information. I am sure there are good reasons to do speed in 0.1 m/s, but personally I find m/s more easy.
The two first variables are ‘vector means’. It is obvious that one cannot simply average directions. Luckily there is the circular package, which does understand direction.
Thus the data reading script becomes:
r2 <- r1[grep(‘^#’,r1):length(r1)]
explain <- r1[1:(grep(‘^#’,r1)-1)]
# explain
r2 <- gsub(‘#’,”,r2)
library(dplyr)
library(circular)
methods(sd)
r4 <- mutate(r3,
Date=as.Date(format(YYYYMMDD),format=’%Y%m%d’),
year=floor(YYYYMMDD/1e4),
month=factor(format(Date,’%B’),levels=month.name),
rDDVEC=as.circular(DDVEC,units=’degrees’,template=’geographics’),
# Vector mean wind direction in degrees
# (360=north, 90=east, 180=south, 270=west, 0=calm/variable)
rFHVEC=FHVEC/10, # Vector mean windspeed (in 0.1 m/s)
rFG=FG/10,   # Daily mean windspeed (in 0.1 m/s)
rFHX=FHX/10, # Maximum hourly mean windspeed (in 0.1 m/s)
rFHN=FHN/10, # Minimum hourly mean windspeed (in 0.1 m/s)
rFXX=FXX/10 # Maximum wind gust (in 0.1 m/s)
) %>%
select(.,YYYYMMDD,Date,year,month,rDDVEC,rFHVEC,
rFG,rFHX,rFHN,rFXX) Plots

Plot of mean wind speed shows several effects. There is an equipment change just before year 2000. At the beginning of the curve the values are lowest, while in the sixties there is a bit more wind, as was n the nineties. I wonder about that. Is that equipment? I can imagine that hundred years ago there was lesser equipment giving such a change, but fifty or twenty years ago? Finally, close to the end of the war there is missing data.
ggplot(data=r4,aes(y=rFG,x=Date))+
geom_smooth()+
geom_point(alpha=.03) +
ylab(‘Mean wind speed x (m/s)’)+
xlab(‘Year’)

A second plot is by month. This shows somewhat different patterns. There is still most wind in the middle of last century. However, September and October have the most wind just before 1950, while November and December have most wind after 1950. Such a pattern cannot be attributed to changes in equipment. It would seem there is some kind of change in wind speeds then. r5 <- group_by(r4,month,year) %>%     summarise(.,mFG=mean(rFG),mFHX=max(rFHX),mFXX=max(rFXX)) ggplot(data=r5,aes(y=mFG,x=year)) +     geom_smooth(method=’loess’) +     geom_point(alpha=.5)+     facet_wrap(~ month)   Wind direction In the Netherlands there is a clear connection between wind and the remainder of the weather. Most of the wind is from the SW (south west, I will be using N, E, S, W to abbreviate directions from here on). N, NW, W and SW winds take humidity from the North Sea and Atlantic Ocean, which in turn will bring rain. In winter, the SW wind will also bring warmth, there will be no frost with W and SW wind. In contrast, N, NE and E will bring cold. A winter wind from Siberia will bring skating fever. In summer, the nice and sunny weather is associated with S to E winds the E wind in May is associated with nice spring weather. SE is by far the least common direction.  The circular package has a both density and plot functions. Combining these gets the following directions for the oldest part of the data.  par(mfrow=c(3,4),mar=c(0,0,3,0)) lapply(month.name,function(x) {       xx <- r4\$rDDVEC[r4\$year<1921 & r4\$month==x]       xx <- xx[!is.na(xx)]       density(xx,bw=50)  %>%            plot(main=x,xlab=”,ylab=”,shrink=1.2)       1     }) title(“1906-1920″, line = -1, outer = TRUE) I would be hard pressed to see significant differences between old and recent data. The densities are slightly different, but not really impressive. Note the lack of E wind in summer, indicating that recent summers have been not been very spectacular.  par(mfrow=c(3,4),mar=c(0,0,2,0)) lapply(month.name,function(x) {       xx <- r4\$rDDVEC[r4\$year>=2000 & r4\$month==x]       xx <- xx[!is.na(xx)]       density(xx,bw=50)  %>%            plot(main=x,xlab=”,ylab=”,shrink=1.2)       1     }) title(“2000-now”, line = -1, outer = TRUE)

To leave a comment for the author, please follow the link and comment on their blog: Wiekvoet. R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

Categories: Methodology Blogs

### Bioenergetics in R Workshop

Sun, 2015-11-15 01:00

(This article was first published on fishR Blog, and kindly contributed to R-bloggers)

It was just brought to my attention that there will be a workshop at the upcoming Midwest Fish and Wildlife Conference (Grand Rapids, MI) on the Bioenergetics 4.0 shiny app. The announcement from here (where there is a registration link) is below (I added the links):

Instructors:

• Dr. David Deslauriers, Post Doctoral Fellow, Department of Biological Sciences, University of Manitoba
• Dr. Steven R. Chipps, Unit Leader, USGS South Dakota Cooperative Fish & Wildlife Research Unit, Department of Natural Resource Management, South Dakota State University

Bioenergetics models are widely used as a tool in fisheries management and research. Although Fish Bioenergetics 3.0 (Hanson et al. 1997) remains a popular software package, it is now over 18 years old and is incompatible with many new operating systems. Moreover, since Fish Bioenergetics 3.0 was released, the number of published fish bioenergetics models has increased from 33 to 115 models. This workshop will introduce Fish Bioenergetics 4.0, an R-based platform that consists of a graphical user interface application (Shiny by RStudio). Instructors will provide an overview of bioenergetics concepts and applications, and introduce attendees to the new modeling platform. Example exercises and group projects will be covered to aid in navigating the software and to answer basic and applied questions in fish ecology.

To leave a comment for the author, please follow the link and comment on their blog: fishR Blog. R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

Categories: Methodology Blogs

### Correlation and Linear Regression

Sat, 2015-11-14 08:44

(This article was first published on DataScience+, and kindly contributed to R-bloggers)

Before going into complex model building, looking at data relation is a sensible step to understand how your different variable interact together. Correlation look at trends shared between two variables, and regression look at causal relation between a predictor (independent variable) and a response (dependent) variable.

Correlation

As mentioned above correlation look at global movement shared between two variables, for example when one variable increases and the other increases as well, then these two variables are said to be positively correlated. The other way round when a variable increase and the other decrease then these two variables are negatively correlated. In the case of no correlation no pattern will be seen between the two variable.

Let’s look at some code before introducing correlation measure:

x<-sample(1:20,20)+rnorm(10,sd=2) y<-x+rnorm(10,sd=3) z<-(sample(1:20,20)/2)+rnorm(20,sd=5) df<-data.frame(x,y,z) plot(df[,1:3])

Here is the plot:

From the plot we get we see that when we plot the variable y with x, the points form some kind of line, when the value of x get bigger the value of y get somehow proportionally bigger too, we can suspect a positive correlation between x and y.

The measure of this correlation is called the coefficient of correlation and can calculated in different ways, the most usual measure is the Pearson coefficient, it is the covariance of the two variable divided by the product of their variance, it is scaled between 1 (for a perfect positive correlation) to -1 (for a perfect negative correlation), 0 would be complete randomness. We can get the Pearson coefficient of correlation using the function cor():

cor(df,method="pearson") x y z x 1.0000000 0.8736874 -0.2485967 y 0.8736874 1.0000000 -0.2376243 z -0.2485967 -0.2376243 1.0000000 cor(df[,1:3],method="spearman") x y z x 1.0000000 0.8887218 -0.3323308 y 0.8887218 1.0000000 -0.2992481 z -0.3323308 -0.2992481 1.0000000

From these outputs our suspicion is confirmed x and y have a high positive correlation, but as always in statistics we can test if this coefficient is significant. Using parametric assumptions (Pearson, dividing the coefficient by its standard error, giving a value that follow a t-distribution) or when data violate parametric assumptions using Spearman rank coefficient.

cor.test(df\$x,df\$y,method="pearson") Pearson's product-moment correlation data: df\$x and df\$y t = 7.6194, df = 18, p-value = 4.872e-07 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.7029411 0.9492172 sample estimates: cor 0.8736874 cor.test(df\$x,df\$y,method="spearman") Spearman's rank correlation rho data: df\$x and df\$y S = 148, p-value < 2.2e-16 alternative hypothesis: true rho is not equal to 0 sample estimates: rho 0.8887218

An extension of the Pearson coefficient of correlation is when we square it we obtain the amount of variation in y explained by x (this is not true for the spearman rank based coefficient where squaring it has no statistical meanings). In our case we have around 75% of the variance in y that is explained by x.

However such results do not allow any causal explanation of the effect of x on y, indeed x could act on y in various way that are not always direct, all we can say from the correlation is that these two variables are linked somehow, to really explain and measure causal effect of x on y we need to use regression method, which will come next.

Linear Regression

Regression is different from correlation because it try to put variables into equation and thus explain causal relationship between them, for example the most simple linear equation is written : Y=aX+b, so for every variation of unit in X, Y value change by aX. Because we are trying to explain natural processes by equations that represent only part of the whole picture we are actually building a model that’s why linear regression are also called linear modelling.

In R we can build and test the significance of linear models.

m1<-lm(mpg~cyl,data=mtcars) summary(m1) Call: lm(formula = mpg ~ cyl, data = mtcars) Residuals: Min 1Q Median 3Q Max -4.9814 -2.1185 0.2217 1.0717 7.5186 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 37.8846 2.0738 18.27 < 2e-16 *** cyl -2.8758 0.3224 -8.92 6.11e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3.206 on 30 degrees of freedom Multiple R-squared: 0.7262, Adjusted R-squared: 0.7171 F-statistic: 79.56 on 1 and 30 DF, p-value: 6.113e-10

The basic function to build linear model (linear regression) in R is to use the lm() function, you provide to it a formula in the form of y~x and optionally a data argument.

Using the summary() function we get all information about our model: the formula called, the distribution of the residuals (the error of our models), the value of the coefficient and their significance plus an information on the overall model performance with the adjusted R-squared (0,71 in our case) that represent the amount of variation in y explained by x, so 71% of the variation in ‘mpg’ can be explain by the variable ‘cyl’.

Before shouting ‘Eureka’ we should first check that the models assumptions are met, indeed linear models make a few assumptions on your data, the first one is that your data are normally distributed, the second one is that the variance in y is homogeneous over all x values (sometimes called homoscedasticity) and independence which means that a y value at a certain x value should not influence other y values.

There is a marvelous built-in methods to check all this with linear models:

par(mfrow=c(2,2)) plot(m1)

The par() argument is just to put all graphs in one window, the plot function is the real one.

Here is the plot:

The graphs on the first columns look at variance homogeneity among other things, normally you should see no pattern in the dots but just a random clouds of points. In this example this is clearly not the case since we see that the spreads of dots increase with higher values of cyl, our homogeneity assumptions is violated we can go back at the beginning and build new models this one cannot be interpreted… Sorry m1 you looked so great…

For the record the graph on the top right check the normality assumptions, if your data are normally distributed the point should fall (more or less) in a straight line, in this case the data are normal.

The final graph show how each y influence the model, each points is removed at a time and the new model is compared to the one with the point, if the point is very influential then it will have a high leverage value. Points with too high leverage value should be removed from the dataset to remove their outlying effect on the model.

Transforming the data

There are a few basics mathematical transformations that can be applied to non normal or heterogeneous data, usually it is a trial and error process;

mtcars\$Mmpg<-log(mtcars\$mpg) plot(Mmpg~cyl,mtcars)

Here is the plot we get:

In our case this looks ok, but we can still remove the two outliers in ‘cyl’ categorie 8;

n<-rownames(mtcars)[mtcars\$Mmpg!=min(mtcars\$Mmpg[mtcars\$cyl==8])] mtcars2<-subset(mtcars,rownames(mtcars)%in%n)

The first line ask for row names in ‘mtcars’ (rownames(mtcars)), but only return the one where the value of the variable ‘Mmpg’ is not equal (!=) to the minimum value of the variable ‘Mmpg’ which fall in the category of 8 cylinders. Then the list ‘n’ contain all these rownames and the next step is to make a new data frame that only contain rows with rownames present in the list ‘n’.

In this stage of transforming and removing outliers from the data you should use and abuse of plots to help you through the process.

Now let’s look back at our bivariate linear regression model from this new dataset;

model<-lm(Mmpg~cyl,mtcars2) summary(model) Call: lm(formula = Mmpg ~ cyl, data = mtcars2) Residuals: Min 1Q Median 3Q Max -0.19859 -0.08576 -0.01887 0.05354 0.26143 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.77183 0.08328 45.292 < 2e-16 *** cyl -0.12746 0.01319 -9.664 2.04e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.1264 on 28 degrees of freedom Multiple R-squared: 0.7693, Adjusted R-squared: 0.7611 F-statistic: 93.39 on 1 and 28 DF, p-value: 2.036e-10 plot(model)

Here is the plot for the model:

Again we have highly significant intercept and slope, the model explain 76% of the variance in log(mpg) and is overall significant. Now we biologist are trained to love and worship ANOVA table, in R there are several way to do it (as always an easy and straightforward way and another with more possibilities for tuning);

anova(model) Analysis of Variance Table Response: Mmpg Df Sum Sq Mean Sq F value Pr(>F) cyl 1 1.49252 1.49252 93.393 2.036e-10 *** Residuals 28 0.44747 0.01598 library(car) Le chargement a nécessité le package : MASS Le chargement a nécessité le package : nnet Anova(model) Anova Table (Type II tests) Response: Mmpg Sum Sq Df F value Pr(>F) cyl 1.49252 1 93.393 2.036e-10 *** Residuals 0.44747 28

The second function Anova() allow you to define which type of sum-of-square you want to calculate (here is a nice explanation of their different assumptions) and also to correct for variance heterogeneity;

Anova(model,white.adjust=TRUE) Analysis of Deviance Table (Type II tests) Response: Mmpg Df F Pr(>F) cyl 1 69.328 4.649e-09 *** Residuals 28

You would have noticed that the p-value is a bit higher. This function is very useful for unbalanced dataset (which is our case) but need to take care when formulating the model especially when there is more than one predictor variables since the type II sum of square assume that there is no interaction between the predictors.

To sum up, correlation is a nice first step to data exploration before going into more serious analysis and to select variable that might be of interest (anyway it always produce sexy and easy to interpret graphs which will make your supervisor happy), then the next step is to model the variable relationship and the most basic models are bivariate linear regression that put the relation between the response variable and the predictor variable into equation and testing this using the summary and anova() function. Since linear regression make several assumptions on the data before interpreting the results of the model you should use the function plot and look if the data are normally distributed, that the variance is homogeneous (no pattern in the residuals~fitted values plot) and when necessary remove outliers.

Next step will be using multiple predictors and looking at generalized linear models.

To leave a comment for the author, please follow the link and comment on their blog: DataScience+. R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

Categories: Methodology Blogs

### Linear model with time series random component

Sat, 2015-11-14 06:00

(This article was first published on Peter's stats stuff - R, and kindly contributed to R-bloggers)

What do auto-correlated residuals do to your linear model?

For training purposes I wanted to illustrate the dangers of ignoring time series characteristics of the random part of a classical linear regression, and I came up with this animation to do it:

I like this, because it shows how easy it is to fit something that looks to be a good fit but actually misses important parts of reality. The red lines show where the fitted model is, based on a small window of the data – from 5 to 200 points. The black line shows the true data generating process. From very early on the model fit to the simple cross-sectional has converged to pretty close to the black line. However, the model fit to the data with time series errors spends a long time greatly overestimating the value of one of the parameters in the model, and not until there are 120 observations has it converged to anywhere near the true process.

At the very least, it shows that you need many more – four times as many in this case, but unfortunately that’s not a magic number that will always work – observations from a time series to reliably estimate the structural part of a model. Even if we’d explicitly modelled the time series part of the data on the right of the animation, we’d still have that problem.

By including the residual plots below the scatter plots we get a nice picture of a warning sign in this basic (and should be fundamental and universal) diagnostic plot. In this particular case the pattern is obvious; when working with real data you should check with partial autocorrelation function plots too.

Simulating data

The animation illustrates the results of simulating and contrasting two fairly extreme cases:

• cross section data, generated exactly from a model of y = a + b.x + e, e ~ N(0, 1). This is the textbook case introduced in any basic statistics course;
• time series data, generated with exactly the same model except the error term, in addition to be normally distributed with mean of zero and standard deviation of 1, has a high autocorrelation.

I chose to make the intercept of my model (a in the above formulation) 1, and the slope (b) equal to 0.3. Here’s what the first 200 observations of the response variable looks like:

In fact, I’ve over-simplified things by leaving x in both datasets as independent and identically distributed white noise. In reality, if y has a time series random component, x probably will have too. But I wanted to illustrate how a single violation of our assumptions can lead to problems, rather than create a fully realistic case (which obviously would show up even more problems).

The data were generated as follows. To illustrate a point and make it a realistic test, I generate a much larger “population” time series, and the mean of zero and standard deviation of 1 applies only to that larger population. The first 200 points is all we see.

library(dplyr) library(tidyr) library(ggplot2) library(gridExtra) library(showtext) #-------------set up------------- # Fonts and themes: font.add.google("Poppins", "myfont") showtext.auto() theme_set(theme_light(base_family = "myfont")) # sample and population size: n <- 200 popn <- n * 10 #----------simulate data--------- set.seed(123) # Linear model with a time series random element, n * 10 in length: df1 <- data.frame(x = rnorm(popn)) %>% mutate(y = 1 + 0.3 * x + scale(arima.sim(list(ar = 0.99), popn)), ind = 1:popn, type = "TimeSeries") # cut back to just the first n points: df1 <- df1[1:n, ] # Same linear model, with i.i.d. white noise random element: df2 <- data.frame(x = rnorm(n)) %>% mutate(y = 1 + 0.3 * x + rnorm(n), ind = 1:n, type = "CrossSection") # draw the time series response: p0 <- df1 %>% ggplot(aes(x = ind, y = y)) + geom_line() + labs(x = "Time") + ggtitle("Simulated response variable from linear modelnwith time series random element")

Creating the animation is straightforward graphics. I make use of ggplot2’s faceting feature to cut down on some code, drawing the top two connected scatterplot images with one chunk and the bottom two residuals with another. Each frame is saved as an individual PNG image, and ImageMagick ties it all together into an animated GIF as easily as usual.

df_both <- rbind(df1, df2) for(i in 5:n){ # I name the images i + 1000 so alphabetical order is also numeric png(paste0(i + 1000, ".png"), 700, 600, res = 100) df1_tmp <- df1[1:i, ] df2_tmp <- df2[1:i, ] residuals1 <- data.frame(res = residuals(lm(y ~ x, data = df1_tmp)), ind = 1:i, type = "TimeSeries") residuals2 <- data.frame(res = residuals(lm(y ~ x, data = df2_tmp)), ind = 1:i, type = "CrossSection") # connected scatter plots: p1 <- ggplot(df_both[c(1:i, (n + 1) : (n + i)), ], aes(x, y, colour = ind)) + facet_wrap(~type, ncol = 2) + geom_path() + geom_point() + geom_abline(intercept = 1, slope = 0.3) + geom_smooth(method = "lm", se = FALSE, size = 2, colour = "red") + theme(legend.position = "none") + xlim(range(df_both\$x)) + ylim(range(df_both\$y)) + ggtitle(paste("Connected scatterplot showing regression on first", i, "points")) # Residuals plots p2 <- residuals1 %>% rbind(residuals2) %>% mutate(type = factor(type, levels = c("CrossSection", "TimeSeries"))) %>% ggplot(aes(x = ind, y = res)) + scale_x_continuous(limits = c(0, n)) + facet_wrap(~type) + geom_line() + geom_point() + ggtitle("Residuals from regression so far") + labs(x = "Time", y = "Residuals") grid.arrange(p1, p2) dev.off() } # combine them into an animated GIF system('"C:\Program Files\ImageMagick-6.9.1-Q16\convert" -loop 0 -delay 10 *.png "timeseries.gif"')

To leave a comment for the author, please follow the link and comment on their blog: Peter's stats stuff - R. R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

Categories: Methodology Blogs

### James Bond movies

Sat, 2015-11-14 06:00

(This article was first published on Opiate for the masses, and kindly contributed to R-bloggers)

James Bond: Do you expect me to talk?
Auric Goldfinger: No, Mr. Bond, I expect you to die!

James Bond

I’m a big James Bond fan, so naturally I went to watch the new Bond movie Spectre which – spoiler alert! – is pretty bad. It also got me to reminice about the good Bond films of the past. My personal candidate for worst Bond film is Die Another Day, but what does the “objective” opinion say on this hotly debated topic? Does my taste conform to the Internet’s taste?

The Economist newspaper did some data analysis when the last Bond (Skyfall) came out, stacking up the different Bond actors on killing, drinking martinis, and love conquests (“Booze, bonks and bodies”). They updated the data for the UK release of Spectre, with Daniel Craig jumping in ranking a lot, mainly from his many kills. The Economist also recently did a comparison of box office opening weekends of the Bond films.

As an economist, I’d have to argue that box office success is one objective measure of film quality; if the movie is bad, people don’t go to the cinemas to watch it, and after all, the market is always right, right?

Others might argue that a critical review scale can assess the quality of a film. This is always debateable: who decides what and how things are grouped into a quality scale? Nevertheless, metacritic and other ratings are used a lot in different industries, because it’s the best you can get.

Data

I decided to pull some data on the Bond movies. Luckily, wikipedia has a nice article, listing all bond films with their respective budget, box office returns, and several critic scales. I’ll use the excellent and easy to use rvest package to pull the data.

## pull oldid of the wikipedia page to ensure reproduceability bond.url <- "https://en.wikipedia.org/w/index.php?title=List_of_James_Bond_films&oldid=688916363" ## read the page into R bond.wiki <- read_html(bond.url) ## film data bond.films <- bond.wiki %>% html_nodes("table") %>% # first table in the page .[[1]] %>% # fill, because wiki uses multi-cell formatting :( html_table(fill = TRUE)

rvest really is that easy to use. It was the first time I used it, and I have to say, I like it a lot. It’ll make pulling data from internet webpages much, much easier in the future for me!

The only problem is that the table in the wiki article has a lot of extra information (like footnotes) that we now need to clean to get a nice, usable dataframe. It’s relatively straightforward, I’ve written a couple of custom functions doing mainly some regexp cleaning. These are the f. functions in the code below. If you’re interested in the details you can check out the full code, including the function code on github.

# lots of cleaning to work to be done now... bond.films %<>% # make clean column names (replace space with dot) setNames(make.names(names(bond.films), unique = TRUE)) %>% # remove first and last line (inner-table headers) head(-1) %>% tail(-1) %>% mutate( # clean wiki-related data import problems Title = f.clean.titles(Title), Bond.actor = f.clean.double.names(Bond.actor), Director = f.clean.double.names(Director), # NA coerced here are what we want -- ignore the warnings Box.office = as.numeric(f.clean.footnote.marks(Box.office)), Budget = as.numeric(f.retrieve.initial.numbers(Budget)), Salary.of.Bond.actor = as.numeric(f.retrieve.initial.numbers(Salary.of.Bond.actor)), # rename the 2005 adjusted price values Box.office.2005.adj = as.numeric(Box.office.1), Budget.2005.adj = as.numeric(Budget.1), Salary.of.Bond.actor.2005.adj = as.numeric(f.retrieve.initial.numbers(Salary.of.Bond.actor.1)), RoI = Box.office.2005.adj/Budget.2005.adj ) %>% # remove old uninformative uniqueness naming select(-contains("1"))

The make.namescommand was new for me, and I am very impressed. It makes setting up proper R-usable variable names a breeze, and in this case was also helpful with uniqueness: The wiki article page uses multi-cell formatting to identify columns, which got lost in my transformations. It’s a bother people don’t conform to clean data standards everywhere!

I then pull and clean the second table as well, which has information on several awards and critic ratings, and merge both for the final dataset to use:

## ratings on films bond.ratings <- bond.wiki %>% html_nodes("table") %>% # second table in the page .[[2]] %>% # fill, because of multi-cell formatting :( html_table(fill = TRUE) bond.ratings %<>% setNames(make.names(names(bond.ratings), unique = TRUE)) %>% mutate( # rename for later merging Title = f.clean.titles(Film), Actor = f.clean.double.names(Actor) ) %>% # reorder to original order select(Title, Year:Awards) %>% # split the RT into rating and number of reviewers separate(Rotten.Tomatoes, c("Rotten.Tomatoes.rating", "Rotten.Tomatoes.reviews"), sep="%") %>% mutate( Rotten.Tomatoes.reviews = sub("^ \((\d*) reviews\).*", "\1", Rotten.Tomatoes.reviews) ) ## final data set bond.dta <- bond.films %>% merge(bond.ratings) %>% select(-Actor) %>% arrange(Year)

The Rotten Tomatoes rating is the one I will be using, by virtue of the fact that it’s available for all films.

Graphing

With all the data pulled, let’s take a look at it more closely.

A little messy, but you can already get the general idea. There was a sharp decline in box office earnings in the late 1980s, and the older Bond movies (with Connery as Bond in particular) have better ratings.

We’ll clean this graph a bit more, and also add the Bond actors. For that, we need to generate a dataset of Bond actors and their time of service.

actor.grp <- bond.dta %>% group_by(Bond.actor) %>% summarise( yearmin = min(Year), yearmax = max(Year) ) %>% ungroup() %>% arrange(yearmin) %>% # George Lazenby only had one film, which makes geom_rect "invisible" mutate( yearmax = ifelse(yearmin==yearmax, yearmax+1, yearmax) )

With this information, we can build a large ggplot graph with all the information in one place!

ggplot() + # place geom_rect first, then other geoms will "write over" the rectangles # write a background rectangle for each Bond actor years of service geom_rect(data = actor.grp, aes(xmin = yearmin, xmax = yearmax, ymin = -Inf, ymax = Inf, fill = Bond.actor), alpha = 0.3)+ # write actor names on rectangles geom_text(data = actor.grp, aes(x = yearmin, # place text rather at the top of the y-axis y = max(bond.dta\$Box.office.2005.adj, na.rm = TRUE), label = Bond.actor, # colour is already mapped to RT.rating (continuous) # fill=Bond.actor, angle = 90, hjust = 1, vjust = 1 ), alpha = 0.6, size = 5)+ # film names geom_text(data = bond.dta, aes(x = Year, y = 0, label = Title, angle = 90, hjust = 0, vjust = 0.5), size=4)+ # film data geom_point(data = bond.dta, aes(x = Year, y = Box.office.2005.adj, size = Budget.2005.adj, colour = as.numeric(Rotten.Tomatoes.rating)))+ # Rotten Tomatoes rating gradient scale_colour_continuous(low="red", high="green", name = "Rotten Tomatoes rating")+ # increase minimum point size for readability scale_size_continuous(name = "Budget (2005 mil. dollars)", range = c(3, 10))+ theme_bw()+ theme(plot.title = element_text(lineheight=.8, face="bold"))+ # remove actor names from legend guides(fill=FALSE)+ labs(title = "Box office results, budgets, and ratings of James Bond filmsn", x="", y="Box office earnings (in 2005 mil. dollars)") # export to size that fits everything into graph, use golden ratio ggsave(file="bond-full.png", width = 30, height = 30/((1+sqrt(5))/2), units = "cm")

Results

It turns out that my personal least favourite Bond isn’t that bad by the scales provided above: it was relatively successful at the Box office, and also wasn’t rated that badly.

The honour of the worst-rated Bond film goes to A View to a Kill (Roger Moore), and the film with the worst box office results is License to Kill (Timothy Dalton). Perhaps it’s the naming policy? It seems you do not make a killing with Bond movies if they have the word “kill” in the title – these two are the only ones.

Roger Moore presided over a continuous decline in popularity in the 1980s, and Timothy Dalton could not stop that trend. This lead to a long pause in Bond films until the franchise was resurrected in 1995 with Pierce Brosnan. Also interesting is the fact that with Brosnan, Bond film budgets noticeably increased in size. Only Moonraker comes close to being as expensive as the later Bonds, and that was set in space.

The early Connery Bond films were the most profitable – easily topping the current Craig films, by up to a factor of 22 (Dr. No yielded 64 times its costs, versus Quantum of Solace which made 2.8). Given the quality of Spectre, I fully expect the current Bond to not be a big success (ratings are already very bad in general). This would lead to a new Bond coming up, if history repeats itself – and Craig has already said he no longer wants to play Bond.

Code and data for this analysis is availabe on github, as always.

James Bond movies was originally published by Kirill Pomogajko at Opiate for the masses on November 14, 2015.

To leave a comment for the author, please follow the link and comment on their blog: Opiate for the masses. R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

Categories: Methodology Blogs

### What it means to be a US Veteran Today

Fri, 2015-11-13 21:56

(This article was first published on Econometrics by Simulation, and kindly contributed to R-bloggers)

Six easy graphs that tell a big story:

1. You represent a much small portion of the American people than veterans in the 1980s.

2. You currently have the highest risk of being classified as poor for any time period since 1980. Since 2005, the rate of poverty among veterans has nearly doubled. That makes you still better off than the general population but not by much.

3. You are now less likely to own your own home than any time since 1980.

4. On a positive note, you are more likely to have completed high-school than at any other time in history.

5. On a not so positive note, as a veteran in 2010+ you are much less likely to have completed four or more years of college than those with no military service. Unfortunately the current world is not friendly to those without a college degree.

6. And to top it off. The strain of being a veteran has negatively affected your marriage. Starting in the 1990s and getting worse over time, the likelihood of being separated or divorced from your spouse is significantly higher than that of the no-military-service population.

So what is the takeaway? Vote for a president you know is going to support your issues.  Source:

In this quick analysis, I look at the census records of 470 thousand random adults between the ages of 18 and 65 sampled each of the years (1980,1990,2000,2005,2010,2013). The source of the data is from IPUMS-USA.

In order control the effect of disproportionate representations of ages across years each year sample has been reduced so that a constant proportions of all ages have been represented for each year.

PUMS-USA, University of Minnesota, www.ipums.org.

To leave a comment for the author, please follow the link and comment on their blog: Econometrics by Simulation. R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

Categories: Methodology Blogs

### Blog Post at Pluralsight

Fri, 2015-11-13 18:57

(This article was first published on R-Chart, and kindly contributed to R-bloggers)

Final post in the three part series is now up at Pluralsight.  The series is geared towards business users so if you have some friends that you are encouraging to set aside their spreadsheets and take R for a spin – send them this way!

Previous posts:

To leave a comment for the author, please follow the link and comment on their blog: R-Chart. R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

Categories: Methodology Blogs

Fri, 2015-11-13 12:49

(This article was first published on Mango Solutions, and kindly contributed to R-bloggers)

By Paulin Shek

Working at Mango is generally busy, fun, but also at times, quite surreal. Lunchtime conversation amongst the consultants can get quite animated, especially when Andy Nicholls, the Head of Consultancy, finds a topic that he has a strong opinion on.

Today, it was chocolate.

In Andy’s mind, chocolate falls into two categories: “Chocolate I would eat after swimming” and “Chocolate I would not eat after swimming”.

I severely questioned his expertise, so now it is up to me to prove him wrong. So off to the supermarket I went, for my all-important “research”.

Andy’s classification of these was as follows – * indicates chocolate bars that Andy was uncertain about.

Now for my part: I collected data from the chocolate bars, collecting fields such as Weight, Calories, Protein etc. which were readily found on the wrapper. I also made flags for whether a chocolate bar has nuts, wafer, caramel etc.

## Weight Calories Fat Carb Protein Salt HasNuts HasWafer ## KitKat Chunky 40.0 516 25.6 65.1 5.4 0.18 0 1 ## Boost 48.5 515 28.5 58.5 5.9 0.30 0 0 ## Dairy Milk 45.0 530 30.5 56.5 7.5 0.23 0 0 ## Galaxy 42.0 546 32.4 56.0 6.7 0.25 0 0 ## Twix 50.0 495 24.0 64.6 4.5 0.44 0 0 ## Picnic 48.4 485 22.5 61.0 7.7 0.53 1 1 ## HasCaramel HasHoneycomb HasNougat ## KitKat Chunky 0 0 0 ## Boost 1 0 0 ## Dairy Milk 0 0 0 ## Galaxy 0 0 0 ## Twix 1 0 0 ## Picnic 1 0 0

Already I was realising that I did not have enough chocolate for the number of fields that I’d collected! Also, did I really want Kinder Buenos with hasNuts=TRUE? I had originally expected this field to be some indicator of the “bulkiness” of the chocolate bar, and Kinder Buenos did not fit with this expectation. In the end, I decided to change the field of hasPeanuts, which excluded Kinder Buenos but kept things like Snickers and Picnic bars. I was very quickly realising the heuristic nature of clustering.

normalise <- function(x){ (x - mean(x))/var(x)} choc <- apply(chocolates, 2, normalise)

Next, I ran the standard k-means clustering algorithm from the stats package. I decided to try 2, 3 and 4 clusters, because I was not convinced that Andy’s binary classification made sense, but also due to the small sample, I couldn’t try too many either.

We could interpret the orange coloured cells to denote “Swim bars”, and this stays consistent regardless of how many clusters we set the k-means algorithm to.

It is interesting to see that Wispas and Wispa Golds are never in the same category, which is not the case in Andy’s categorisation. However, it makes intuitive sense to others- another consultant Aimee said “Well, if I wanted to eat a Wispa, a Wispa Gold would not be a suitable substitute”. So there we go.

Taking K=3, we can see that the Crunchie has been identified as an outlier, but then it is hard to make sense of the blue group from 4 clusters!

I decided that perhaps normalising every column was not the best idea- the binary values definitely needed to be adjusted, but I was less certain about things like salt. I tried running the cluster analysis again, starting again from the original data, but adjusting the binary variables (HasCaramel etc.) to be much bigger, multiplying them by a scalar that’s on par with the other values in the data.

This time K=2 looks less reasonable, with Picnic and Snickers bars not grouped with the Mars bar. K=3 looks a bit better- the white group could almost be a “plain chocolate” group.

I think I’ve found the perfect classification with K=4 though. The blue group are clearly the “Wafer” group, all the chocolate only bars are picked out by the purple group, and them Crunchie and Wispa Golds are put together as the “super sweet” group.

So, I’ve found the perfect classification, (which is clearly better than Andy’s!). I can stop and reward myself with one of my many chocolate bars now!

At a glance of the data, I could already tell that the 1, 0 values were going to be too small and hence get ignored by the clustering. So, I decided to normalise all the values.

To leave a comment for the author, please follow the link and comment on their blog: Mango Solutions. R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

Categories: Methodology Blogs

### Annotables: R data package for annotating/converting Gene IDs

Fri, 2015-11-13 10:54

(This article was first published on Getting Genetics Done, and kindly contributed to R-bloggers)

I work with gene lists on a nearly daily basis. Lists of genes near ChIP-seq peaks, lists of genes closest to a GWAS hit, lists of differentially expressed genes or transcripts from an RNA-seq experiment, lists of genes involved in certain pathways, etc. And lots of times I’ll need to convert these gene IDs from one identifier to another. There’s no shortage of tools to do this. I use Ensembl Biomart. But I do this so often that I got tired of hammering Ensembl’s servers whenever I wanted to convert from Ensembl to Entrez gene IDs for pathway mapping, get the chromosomal location for some BEDTools-y kinds of genomic arithmetic, or get the gene symbol and full description for reporting. So I used Biomart to retrieve the data that I use most often, cleaned up the column names, and saved this data as an R data package called annotables. This package has basic annotation information from Ensembl release 82 for:
• Human (grch38)
• Mouse (grcm38)
• Rat (rnor6)
• Chicken (galgal4)
• Worm (wbcel235)
• Fly (bdgp6)
Where each table contains:
• ensgene: Ensembl gene ID
• entrez: Entrez gene ID
• symbol: Gene symbol
• chr: Chromosome
• start: Start
• end: End
• strand: Strand
• biotype: Protein coding, pseudogene, mitochondrial tRNA, etc.
• description: Full gene name/description.
Additionally, there are tables for human and mouse (grch38_gt and grcm38_gt, respectively) that link ensembl gene IDs to ensembl transcript IDs. Usage The package isn’t on CRAN, so you’ll need devtools to install it. # If you haven't already installed devtools...
install.packages("devtools")

# Use devtools to install the package
devtools::install_github("stephenturner/annotables")
It isn’t necessary to load dplyr, but the tables are tbl_df and will print nicely if you have dplyr loaded. library(dplyr)
library(annotables)
Look at the human genes table (note the description column gets cut off because the table becomes too wide to print nicely): grch38
## Source: local data frame [66,531 x 9]
##
## ensgene entrez symbol chr start end strand biotype
## (chr) (int) (chr) (chr) (int) (int) (int) (chr)
## 1 ENSG00000210049 NA MT-TF MT 577 647 1 Mt_tRNA
## 2 ENSG00000211459 NA MT-RNR1 MT 648 1601 1 Mt_rRNA
## 3 ENSG00000210077 NA MT-TV MT 1602 1670 1 Mt_tRNA
## 4 ENSG00000210082 NA MT-RNR2 MT 1671 3229 1 Mt_rRNA
## 5 ENSG00000209082 NA MT-TL1 MT 3230 3304 1 Mt_tRNA
## 6 ENSG00000198888 4535 MT-ND1 MT 3307 4262 1 protein_coding
## 7 ENSG00000210100 NA MT-TI MT 4263 4331 1 Mt_tRNA
## 8 ENSG00000210107 NA MT-TQ MT 4329 4400 -1 Mt_tRNA
## 9 ENSG00000210112 NA MT-TM MT 4402 4469 1 Mt_tRNA
## 10 ENSG00000198763 4536 MT-ND2 MT 4470 5511 1 protein_coding
## .. ... ... ... ... ... ... ... ...
## Variables not shown: description (chr)
Look at the human genes-to-transcripts table: grch38_gt
## Source: local data frame [216,133 x 2]
##
## ensgene enstxp
## (chr) (chr)
## 1 ENSG00000210049 ENST00000387314
## 2 ENSG00000211459 ENST00000389680
## 3 ENSG00000210077 ENST00000387342
## 4 ENSG00000210082 ENST00000387347
## 5 ENSG00000209082 ENST00000386347
## 6 ENSG00000198888 ENST00000361390
## 7 ENSG00000210100 ENST00000387365
## 8 ENSG00000210107 ENST00000387372
## 9 ENSG00000210112 ENST00000387377
## 10 ENSG00000198763 ENST00000361453
## .. ... ...
Tables are tbl_df, pipe-able with dplyr: grch38 %>%
filter(biotype=="protein_coding" & chr=="1") %>%
select(ensgene, symbol, chr, start, end, description) %>%
pander::pandoc.table(split.table=100, justify="llllll", style="rmarkdown")
ensgene symbol chr start end ENSG00000158014 SLC30A2 1 26037252 26046133 ENSG00000173673 HES3 1 6244192 6245578 ENSG00000243749 ZMYM6NB 1 34981535 34985353 ENSG00000189410 SH2D5 1 20719732 20732837 ENSG00000116863 ADPRHL2 1 36088875 36093932 ENSG00000188643 S100A16 1 153606886 153613145 Table: Table continues below description solute carrier family 30 (zinc transporter), member 2 [Source:HGNC Symbol;Acc:HGNC:11013] hes family bHLH transcription factor 3 [Source:HGNC Symbol;Acc:HGNC:26226] ZMYM6 neighbor [Source:HGNC Symbol;Acc:HGNC:40021] SH2 domain containing 5 [Source:HGNC Symbol;Acc:HGNC:28819] ADP-ribosylhydrolase like 2 [Source:HGNC Symbol;Acc:HGNC:21304] S100 calcium binding protein A16 [Source:HGNC Symbol;Acc:HGNC:20441] Example with RNA-seq data Here’s an example with RNA-seq data. Specifically, DESeq2 results from the airway package, made tidy with biobroom: # Load libraries (install with Bioconductor if you don't have them)
library(DESeq2)
library(airway)

# Load the data and do the RNA-seq data analysis
data(airway)
airway = DESeqDataSet(airway, design = ~cell + dex)
airway = DESeq(airway)
res = results(airway)

# tidy results with biobroom
library(biobroom)
res_tidy = tidy.DESeqResults(res)
## Source: local data frame [6 x 7]
##
## gene baseMean estimate stderror statistic
## (chr) (dbl) (dbl) (dbl) (dbl)
## 1 ENSG00000000003 708.6021697 0.37424998 0.09873107 3.7906000
## 2 ENSG00000000005 0.0000000 NA NA NA
## 3 ENSG00000000419 520.2979006 -0.20215550 0.10929899 -1.8495642
## 4 ENSG00000000457 237.1630368 -0.03624826 0.13684258 -0.2648902
## 5 ENSG00000000460 57.9326331 0.08523370 0.24654400 0.3457140
## 6 ENSG00000000938 0.3180984 0.11555962 0.14630523 0.7898530
## Variables not shown: p.value (dbl), p.adjusted (dbl)
Now, make a table with the results (unfortunately, it’ll be split in this display, but you can write this to file to see all the columns in a single row): res_tidy %>%
inner_join(grch38, by=c("gene"="ensgene")) %>%
select(gene, estimate, p.adjusted, symbol, description) %>%
pander::pandoc.table(split.table=100, justify="lrrll", style="rmarkdown")
gene estimate p.adjusted symbol ENSG00000152583 -4.316 4.753e-134 SPARCL1 ENSG00000165995 -3.189 1.44e-133 CACNB2 ENSG00000101347 -3.618 6.619e-125 SAMHD1 ENSG00000120129 -2.871 6.619e-125 DUSP1 ENSG00000189221 -3.231 9.468e-119 MAOA ENSG00000211445 -3.553 3.94e-107 GPX3 ENSG00000157214 -1.949 8.74e-102 STEAP2 ENSG00000162614 -2.003 3.052e-98 NEXN ENSG00000125148 -2.167 1.783e-92 MT2A ENSG00000154734 -2.286 4.522e-86 ADAMTS1 ENSG00000139132 -2.181 2.501e-83 FGD4 ENSG00000162493 -1.858 4.215e-83 PDPN ENSG00000162692 3.453 3.563e-82 VCAM1 ENSG00000179094 -3.044 1.199e-81 PER1 ENSG00000134243 -2.149 2.73e-81 SORT1 ENSG00000163884 -4.079 1.073e-80 KLF15 ENSG00000178695 2.446 6.275e-75 KCTD12 ENSG00000146250 2.64 1.143e-69 PRSS35 ENSG00000198624 -2.784 1.707e-69 CCDC69 ENSG00000148848 1.783 1.762e-69 ADAM12 Table: Table continues below description SPARC-like 1 (hevin) [Source:HGNC Symbol;Acc:HGNC:11220] calcium channel, voltage-dependent, beta 2 subunit [Source:HGNC Symbol;Acc:HGNC:1402] SAM domain and HD domain 1 [Source:HGNC Symbol;Acc:HGNC:15925] dual specificity phosphatase 1 [Source:HGNC Symbol;Acc:HGNC:3064] monoamine oxidase A [Source:HGNC Symbol;Acc:HGNC:6833] glutathione peroxidase 3 [Source:HGNC Symbol;Acc:HGNC:4555] STEAP family member 2, metalloreductase [Source:HGNC Symbol;Acc:HGNC:17885] nexilin (F actin binding protein) [Source:HGNC Symbol;Acc:HGNC:29557] metallothionein 2A [Source:HGNC Symbol;Acc:HGNC:7406] ADAM metallopeptidase with thrombospondin type 1 motif, 1 [Source:HGNC Symbol;Acc:HGNC:217] FYVE, RhoGEF and PH domain containing 4 [Source:HGNC Symbol;Acc:HGNC:19125] podoplanin [Source:HGNC Symbol;Acc:HGNC:29602] vascular cell adhesion molecule 1 [Source:HGNC Symbol;Acc:HGNC:12663] period circadian clock 1 [Source:HGNC Symbol;Acc:HGNC:8845] sortilin 1 [Source:HGNC Symbol;Acc:HGNC:11186] Kruppel-like factor 15 [Source:HGNC Symbol;Acc:HGNC:14536] potassium channel tetramerization domain containing 12 [Source:HGNC Symbol;Acc:HGNC:14678] protease, serine, 35 [Source:HGNC Symbol;Acc:HGNC:21387] coiled-coil domain containing 69 [Source:HGNC Symbol;Acc:HGNC:24487] ADAM metallopeptidase domain 12 [Source:HGNC Symbol;Acc:HGNC:190] Explore! This data can also be used for toying around with dplyr verbs and generally getting a sense of what’s in here. First, tet some help. ls("package:annotables")
?grch38
Let’s join the transcript table to the gene table. gt = grch38_gt %>%
inner_join(grch38, by="ensgene")
Now, let’s filter to get only protein-coding genes, group by the ensembl gene ID, summarize to count how many transcripts are in each gene, inner join that result back to the original gene list, so we can select out only the gene, number of transcripts, symbol, and description, mutate the description column so that it isn’t so wide that it’ll break the display, arrange the returned data descending by the number of transcripts per gene, head to get the top 10 results, and optionally, pipe that to further utilities to output a nice HTML table. gt %>%
filter(biotype=="protein_coding") %>%
group_by(ensgene) %>%
summarize(ntxps=n_distinct(enstxp)) %>%
inner_join(grch38, by="ensgene") %>%
select(ensgene, ntxps, symbol, description) %>%
mutate(description=substr(description, 1, 20)) %>%
arrange(desc(ntxps)) %>%
pander::pandoc.table(split.table=100, justify="lrll", style="rmarkdown")
ensgene ntxps symbol description ENSG00000165795 77 NDRG2 NDRG family member 2 ENSG00000205336 77 ADGRG1 adhesion G protein-c ENSG00000196628 75 TCF4 transcription factor ENSG00000161249 68 DMKN dermokine [Source:HG ENSG00000154556 64 SORBS2 sorbin and SH3 domai ENSG00000166444 62 ST5 suppression of tumor ENSG00000204580 58 DDR1 discoidin domain rec ENSG00000087460 57 GNAS GNAS complex locus [ ENSG00000169398 57 PTK2 protein tyrosine kin ENSG00000104529 56 EEF1D eukaryotic translati Let’s look up DMKN (dermkine) in Ensembl. Search Ensembl for ENSG00000161249, or use this direct link. You can browse the table or graphic to see the splicing complexity in this gene. Or, let’s do something different. Let’s group the data by what type of gene it is (e.g., protein coding, pseudogene, etc), get the number of genes in each category, and plot the top 20. library(ggplot2)
grch38 %>%
group_by(biotype) %>%
summarize(n=n_distinct(ensgene)) %>%
arrange(desc(n)) %>%