## Quantum Forest

### notes in a shoebox

#### Category: tutorials

The Swarm Lab presents a nice comparison of R and Python code for a simple (read ‘one could do it in Excel’) problem. The example works, syringe but I was surprised by how wordy the R code was and decided to check if one could easily produce a shorter version.

The beginning is pretty much the same, this although I’ll use ggplot2 rather than lattice, because it will be a lot easier (and shorter) to get the desired appearance for the plots:

The whole example relies on only three variables and—as I am not great at typing—I tend to work with shorter variable names. I directly changed the names for variables 1 to 3:

It is a lot easier to compare the regression lines if we change the shape of the data set from wide to long, where there is one variable for year, one for tax type, and one for the actual tax rate. It would be possible to use one of Hadley’s packages to get a simpler syntax for this, but I decided to stick to the minimum set of requirements:

And now we are ready to produce the plots. The first one can be a rough cut to see if we get the right elements:

First cut of the taxes per year plot.

Yes, this one has the points, lines, linear regression and 95% confidence intervals for the mean predicted responses, but we still need to get rid of the grey background and get black labels (theme_bw()), set the right axis labels and ticks (scale_x... scale_y...) and set the right color palette for points and lines (scale_colour_manual) and filling the confidence intervals (scale_colour_fill) like so:

Way closer to the desired plot, still much shorter.

One can still change font sizes to match the original plots, reposition the legend, change the aspect ratio while saving the png graphs (all simple statements) but you get the idea. If now we move to fitting the regression lines:

This gives the regression coefficients for Corporate (3.45 – 1.564e-04 year) and Individual ([3.45 + 4.41] + [-1.564e-04 + 1.822e-04] year or 7.84 + 2.58e-05 year). As a bonus you get the comparison between regression lines.

In R as a second language I pointed out that ‘brevity reduces the time between thinking and implementation, so we can move on and keep on trying new ideas’. Some times it seriously does.

Disappeared for a while collecting frequent flyer points. In the process I ‘discovered’ that I live pharm +Canterbury&hl=en&ll=-43.580391,172.617188&spn=171.212657,59.0625&sll=-43.068888,171.914063&sspn=171.285878,59.0625&oq=chris&hnear=Christchurch,+Canterbury&t=m&z=2″>in the middle of nowhere, as it took me 36 hours to reach my conference destination (Estoril, Portugal) through Christchurch, Sydney, Bangkok, Dubai, Madrid and Lisbon.

Where was I? Showing how split-plots look like under the bonnet (hood for you US readers). Yates presented a nice diagram of his oats data set in the paper, so we have the spatial location of each data point which permits us playing with within-trial spatial trends.

Rather than mucking around with typing coordinates we can rely on Kevin Wright’s version of the oats dataset contained in the agridat package. Kevin is a man of mistery, a James Bond of statisticians—so he keeps a low profile—with a keen interest in experimental design and analyses. This chap has put a really nice collection of data sets WITH suggested coding for the analyses, including nlme, lme4, asreml, MCMCglmm and a few other bits and pieces. Recommended!

Plants (triffids excepted) do not move, which means that environmental trends within a trial (things like fertility, water availability, temperature, etc) can affect experimental units in a way that varies with space and which induces correlation of the residuals. Kind of we could be violating the independence assumption if you haven’t got the hint yet.

Gratuitous picture: Detail of Mosteiro dos Jerónimos, Belém, Lisboa, Portugal (Photo: Luis).

There are a few ways to model environmental trends (AR processes, simple polynomials, splines, etc) that can be accounted for either through the G matrix (as random effects) or the R matrix. See previous post for explanation of the bits and pieces. We will use here a very popular approach, which is to consider two separable (so we can estimate the bloody things) autoregressive processes, one for rows and one for columns, to model spatial association. In addition, we will have a spatial residual. In summary, the residuals have moved from $$oldsymbol{R} = sigma^2_e oldsymbol{I}$$ to $$oldsymbol{R} = sigma^2_s oldsymbol{R}_{col} otimes oldsymbol{R}_{row}$$. I previously showed the general form of this autoregressive matrices in this post, and you can see the $$oldsymbol{R}_{col}$$ matrix below. In some cases we can also add an independent residual (the so-called nugget) to the residual matrix.

We will first fit a split-plot model considering spatial residuals using asreml because, let’s face it, there is no other package that will give you the flexibility:

So we have to build an autoregressive correlation matrix for rows, one for columns and multiply the whole thing for a spatial variance. Then we can add an independent residual (the nugget, if we want—and can estimate—one). Peter Dalgaard has neat code for building the autocorrelation matrix. And going back to the code in the previous post:

Which are the same results one gets from ASReml-R. QED.

P.S. Many thanks to Kevin Wright for answering my questions about agridat.

A commenter on this blog reminded me of one of the frustrating aspects faced by newbies, dosage not only to R but to any other programming environment (I am thinking of typical students doing stats for the first time). The statement “R is a language” sounds perfectly harmless if you have previous exposure to programming. However, visit this if you come from a zero-programming background the question is What do you really mean?

R—and many other statistical systems for that matter—is often approached from either one of two extremes:

1. Explaining R as a programming language, as most of the R documentation and books (like The Art of R Programming, quite good by the way) do.
2. The other one is from a hodgepodge of statistical analyses, introducing the language as a bonus, best represented by Crowley’s The R Book (which I find close to unreadable). In contrast, Modern Applied Statistics with S by Ripley and Venables is much better even when it doesn’t mention R in the title.

If you are new to both statistics and R I like the level of the Quick-R website as a starting point, which was expanded into a book (R in Action). It uses the second approach listed above, so if you come from a programming background the book will most likely be disappointing. Nevertheless, if you come from a newbie point of view both the website and book are great resources. In spite of this, Quick-R assumes that the reader is familiar with statistics and starts with “R is an elegant and comprehensive statistical and graphical programming language.

## A simpler starting point

I would like to start from an even simpler point, ignoring for a moment programming and think about languages like English, Spanish, etc. In languages we have things (nouns) and actions (verbs). We perform actions on things: we measure a tree, draw a plot, make some assumptions, estimate coefficients, etc. In R we use functions to perform actions on objects (things in the previous explanations). Our data sets are objects that we read, write, fit a model (and create objects with results), etc. “R is a language” means that we have a grammar that is designed to deal with data from a statistical point of view.

A simple sentence “Luis writes Quantum Forest” has two objects (Luis and Quantum Forest) and one function (writes). Now lets look at some simple objects in R; for example, a number, a string of characters and a collection of numbers (the latter using the function c() to keep the numbers together):

Up to this point we have pretty boring software, but things start becoming more interesting when we can assign objects to names, so we can keep acting on those objects (using functions). In this blog I use = instead of <- to assign things (objects) to names outside function calls. This is considered "bad form" in the R world, but to me is much more readable§. (Inside function calls the arguments should always be referred with an = sign, as we'll see in a future post). Anyway, if you are feeling in a conformist mood replace the = by <- and the code will work equally well.

R is case sensitive, meaning that upper- and lower-case letters are considered different. Thus, we can assign different objects to variables named sex, Sex and SEX and R will keep track of them as separate entities (A Kama Sutra of names!). Once objects are assigned to a variable R stops printing the object back to the user. However, it is possible to type the object name, press enter and get the content stored in the name. For example:

The lowly = sign is a function as well. For example, both a and b are assigned the same bunch of numbers:

Even referring to an object by its name calls a function! print(), which is why we get [1] 23.000000 44.000000 3.141593 when typing b and hitting enter in R.

## Grouping objects

Robert Kabacoff has a nice invited post explaining data frames. Here I will present a very rough explanation with a toy example.

Objects can be collected in other objects and assigned a name. In data analysis we tend to collect several variables (for example tree height and stem diameter, people's age and income, etc). It is convenient to keep variables referring to the same units (trees, persons) together in a table. If you have used Excel, a rough approximation would be a spreadsheet. Our toy example could be like:

The last line combines two objects (x and y) in an R table using the function data.frame() and then it assigns the name "toy" to that table (using the function =). From now on we can refer to that data table when using other functions as in:

Incidentally, anything following a # is a comment, which helps users document their work. Use them liberally.

Fitting a linear regression will produce lots of different outputs: estimated regression coefficients, fitted values, residuals, etc. Thus, it is very handy to assign the results of the regression to a name (in this case "fm1") for further manipulation. For example:

The last line of code extracts the residuals of fm1 (using the function resid()) and the fitted values (using the function fitted()), which are then combined using the function plot().

In summary: in this introduction we relied on the R language to manipulate objects using functions. Assigning names to objects (and to the results of applying functions) we can continue processing data and improving our understanding of the problem under study.

## Footnotes

R is an implementation (a dialect) of the S language. But remember, a language is a dialect with an army and a navy.

Natural languages tend to be more complex and will have pronouns, articles, adjectives, etc. Let's ignore that for the moment.

§Languages change; for example, I speak Spanish—a bad form of Latin—together with hundreds of millions of people. Who speaks Latin today?