Notice that variables like Genus, Biome, Species... are all called factors. This means they are categorical variables. R will automatically assign categorical status to your data, but it can sometimes make mistakes. For example, suppose BIOME were actually not a factor - let's force R to do that temporarily:

Now we can convert it back:

Take a look at ?factor and ?levels for more information about this. You can define ordered and unordered factors for any set of levels. Checking levels is always good - you may have had a typo in your datasheet that will add an extra level to your analysis and reduce power!

Let's try learning some simple things about our data - what, for example, is the average log.LMA (log10 leaf mass per leaf area). Here we need the $ operator to access a column of the data frame. We'll use the mean command to get an average.

That didn't work, because some rows of the dataset have missing values for log.LMA. We can find them and return the rows of the data that have NA values. The which commands returns the indices of the TRUE entries of a vector.

let's calculate the mean on only the non-NA values of log LMA now, by first getting a subset of the data. We use [] to access different rows and columns of the dataset. The ! operator means NOT (i.e. converts true to false, false to true)

Now we have an answer - but there is an easier way too. Read the help for the mean command, and you'll find that it can automatically remove missing values.

We can also get standard deviation, median, etc. the same way

maybe more broadly we are curious about the coverage of the dataset - does it include as many deciduous as evergreen species? We can count these using the table function:

how many rows does our dataset have?

We can also count the number of items in each individual column/row vector using 'length':

This dataset has a column which assigns a categorical variable (BIOME) to each point. Suppose we wanted to find the mean logLMA in each biome. Here we applied the 'mean' function to log.LMA for each unique value found in BIOME.

Suppose we only wanted to study oak (Quercus) species in temperate forests. Let's subset the data to find only these points. First, what is the temperate forest biome called?

Now, let's use the == (equality test) operator along with the & (logical and) operator to get the subset

These temperate oaks can have either deciduous or evergreen leaves. Let's look at the LMA of the evergreen ones using the [] and $ operators together. First we will choose all rows that have evergreen leaves using the [,] operator (leaving the space before or after the comma means to take all values, i.e. [,x] is all rows satisfying x while [x,] is all columns satisfying x.

We can also do more complicated things to different categories within a dataframe. Suppose we wanted to calculate the ratio of deciduous to evergreen leaf lifespans within each biome. This is not too hard - we can use the 'by' command and then write a custom function that will return the right ratio. We can define a new function just like defining a variable:

Now we will use the 'by' command to apply this function to each unique BIOME of the data. We then have to call 'unlist' to transform the result from a list (another data type) to a data frame (much easier to work with)

Now let's sort in decreasing order the output to see which biome has the highest ratio of leaf lifespans So in tropical forests leaf lifespans are most equivalent between deciduous and evergreen species!

Suppose we wanted something simpler, like the name of the species with the highest value of leaf lifespan:

That returned the entire row - let's output the full name of the species and the actual lifespan (un-logged) of that species' leaves: 288 months - 24 years is really long! Check it out at http://www.esu.edu/~milewski/intro_biol_two/lab_3_seed_plts/Araucaria.html

Let's try using the [,] and $ operators more. These are very useful for subsetting and manipulating your data! should give the same output

all rows in the Genus column

Use a logical test to select rows, and take two columns of the result

Only get a single column this time

Sometimes it is useful to be able to sort a dataframe by multiple columns. Let's sort by biome, then by decreasing LL. We use the 'order' command to generate a list of indices that we then use as row IDs: a lot of output, but you can see it's ordered first by BIOME, then by LL!

We can also get a random subset of rows - here, choose 10 random values from the integers ranging from 1 to nrow(glopnet) Try running these two previous commands again - you should get different values this time.

We can get the mean log leaf lifespan from a random subset of the data:

Is there a significant difference in leaf lifespan between deciduous and evergreen species? First, let's get the right data:

Now let's try a two-sided t-test at the alpha = 0.05 level for no difference in means:

Of course we really should have tested to see if the data met the assumptions of a t test, the most important of which being normality. We can test normality using a Shapiro-Wilk test: Not looking so good!

Let's do a Mann-Whitney (Wilcoxon Rank Sum) non-parametric test instead, that doesn't assume data are normally distributed And the difference is still significant!

We can also test for correlations between variables. For example, perhaps leaf lifespan is significantly correlated with leaf nitrogen content?

We can also directly access different parts of the test result - for example, the p value:

You can find documentation on many more of these tests by googling. We're going to focus instead on linear models, which encompass many of these tests and also include things such as ANOVA, ANCOVA, MANOVA.

The most important command to learn is lm, which generates linear models. These models are assumed to represent data points with normally distributed errors, homoscedastic variances, independent data points, random sampling... etc. Your data may not match these assumptions (for example, count data (e.g. number of marmots seen today) will probably have Poisson-distributed errors, and presence/absence data may require a logistic regression) and you may have multiple non-independent points in your dataset. But don't worry! R can handle these things too (scroll down). Let's just explore the simple cases first.

Let's see if there is a relationship between photosynthetic rate (log.Aarea) and nitrogen content (log.Nmass). This might be expected because photosynthesis depends on the Rubisco protein, which is very nitrogen-rich. Notice how we specify a formula with the dependent variable on the left of the ~ and the independent variable on the right. This prints the coefficients of the regression A much more useful output. Looks like a significant positive relationship! Try str(summary(m_an)) to see what else we can extract from the model summary.

It's very important to verify that the data meet all the assumptions of a linear model. We can check residuals, leverage, etc: This will generate several plots you can step through with the mouse click or the return key on the console. Probably the simplest diagnostic is the normal quantile-quantile plot - the quantiles of your data should match the quantiles of a scaled and translated normal distribution, generating a straight line. If you don't want to see all the linear model diagnostics, just try which looks pretty good.

This last model included by default an intercept term (in the language of linear models, Y ~ B0 + B1 X). You might not want an intercept term, so try instead adding '-1' to the formula:

Creating a linear model object lets us access a number of useful additional variables like
the model coefficients the residuals at each data point confidence intervals the effects

We can also specify much more complex formulas. You should know the I(), :, and * notation. * adds terms that include linear and interaction terms. : adds only interaction terms. I adds terms 'as-is', which are numerically determined -before- running the linear model (for example, adding a quadratic term to the model. Let's try a few: You can try removing the I from the above model - notice it won't keep a quadratic term for you!

We can also perform an ANOVA to look at sums of squares and see which of the interactions and terms are overall significant in a model: This approach is especially useful if we are using categorical variables. In R, you can treat categorical variables just like continuous ones. Let's see if there is an interaction between BIOME and leaf nitrogen that influences Aarea: yes, a significant interaction was found!
You do need to be careful about something here - this is type I ANOVA that works by sequentially adding terms to look at the effect on the residual sum of squares. So the order in which you specify the variables will potentially impact your conclusions.

You need to think about type II or type III ANOVA if you don't think this is appropriate. Unbalanced designs will also require other types of ANOVAs. Fortunately this is possible, but we need to bring in some extra help. You will want to install the extra 'car' package for R by going up to the Packages & Data Menu, and selecting Package Installer. There are hundreds of packages which extend R's functionality. All of these are freely available, and are documented to varying extents. Some work better than others - google before using! 'car' is very well-respected. Once you have done the one-time download of the package, you need to load it into this R session:

We will now use the Anova command in 'car' (note the new capitalization) to try a type III sum of squares: compare this to anova(m_threeway) - notice the different results! You can try ?Anova to understand more what is going on.

You may also want to explore the 'aov' function for performaning analyses of variances. However this requires a balanced data set - I recommend using 'lme' instead (see later). 'lm' can also handle contrasts if you are into those.

We can also compare between models using the standard 'anova' function (not Anova in 'car') There is a very tricky thing going on here - how do we decide how to compare models? You can do a F test, Chi-square test, or a Mallow's Cp test. Think carefully about what you want.

Probably a better way to compare models is to use the Akaike Information Criterion (a penalized log-likelihood). For a given dependent variable and dataset, we can assess which set of predictors (i.e. which model) is best. Lower AIC values are better, but there is no absolute scale - AIC only allows comparison between candidate models. We can also programmatically find the minimium AIC (i.e. the best model)

We can also programmatically find a good model using AIC-based model selection, stepwise deleting terms from a model until one is found where no further deletions improve the AIC: which shows the 'best' model by an AIC criterion

And it's also easy to use a linear model for future prediction. Suppose we wanted to make a prediction at some new data points. We need to define a new data frame with all of the required predictors:

Ordinary least-squares may not always be the most appropriate type of regression. While OLS estimators are maximum likelihood estimators under many conditions, they are not ideal under all conditions. For example, OLS is very sensitive to outlying points, so that minimizing squared residuals may not work well for smaller high-variance samples. In this case, you may want to instead use robust regression in the MASS package (you will need to install it). Here you can see for the relatively unsampled non-Angiosperm taxa, the robust-regression method estimates a higher intercept and lower slope. This estimate may be more useful if the data have known but unremovable outliers.

OLS regression (Type I regression) is useful if you want to compare the likelihoods of different models, or make predictions. However if you are more interested in obtaing a true estimate of regression coefficients (e.g. slope) and/or your data have uncertainty in both x and y, you may want to use Type II regression instead. Type II regression will minimize (optionally rescaled) perpendicular, rather than just vertical residuals. This is called major axis or standard/reduced major axis regression. Fortunately there is a package for that (you will need to install it first)

the smatr package doesn't let you enter formulas or use a data= argument - so we have to manually specify y and x variables. see ?line.cis for details. Notice that this method estimates a higher slope than the OLS method did.

We can also do slope/intercept tests on Type II regressions. Suppose we hypothesize that the true slope of the A/N relationship is 1, because photosynthesis is proportional to nitrogen (maybe) We get p=0, so the true slope is not 1!

Let's move onto more complex types of linear models. Suppose our data were not all independent points but rather was psuedoreplicated. To deal with this, and more generally with unbalanced designs, we need to start working with mixed models. Let's first load the nlme package use the lme command for linear mixed models. We can now specify fixed and nested random effects. Suppose (wrongly) that we nee random effects of species gropued within genus when considering the relationship between Aarea and Nmass:

You can also specify random slopes, random intercepts at different levels of hierarchy... but the syntax gets very complex very quickly, and documentation is poor. You can consider also using the 'lmer' command in the 'lme4' library. Mixed models are controversial and poorly documented in R unfortunately.

Let's instead consider generalized linear models, which allow non-normally distributed errors. For example, suppose we were considering Aarea instead of logAarea - Aarea varies from 0 to infinity, so might be best described as Gamma-distributed. We can account for this quite easily: Notice how it calculates a dispersion for the gamma distribution and deals with missing values automatically...
You can learn more about other more complex models (mostly exponential family distributions) with ?family. Notice the options for probit, logit, logistic regression if your data merit these options. GLM is a very powerful way to analyze your data!

Another very useful command is gls in the nlme package, for generalized least squares. This function lets you break the assumption of uncorrelated data - for example, dealing with phylogenetic correlation, temporal autocorrelation, etc. Here's an example of using GLS to account for phylogenetic structure (an analysis equivalent to creating phylogenetic independent contrasts). We have data on leaf area shrinkage (after drying) and LMA for multiple species, at the leaf level. First let's calculate species level means:

Now let's read in a phylogenetic tree that contains these species using a phylogenetic analysis package (ape, please install it!)

Let's compare the model with and without phylogenetic independent contrasts (Brownian evolution correlation structure) In both models, the shrinkage effect is significant, but is slightly stronger (57 vs 56) when phylogeny is taken into account.) So in this case the relationship is probably not strongly affected by evolutionary history. You may also want to use gls if you have spatially or temporally autocorrelated data.

R has packages that can do many other types of modeling (classification/regression trees, generalized additive models, generalized linear mixed models, LOESS, quantile regression, whatever you like - to find packages to do what you want, look on CRAN (http://cran.r-project.org/web/packages/)

One last very important function is principal components analysis (PCA). Note that R also has packages for things like non-metric multidimensional scaling, canonical correlation analysis, redundancy analysis... you will again want to look at CRAN for these more complex ordination techniques. But for simple multidimensional descriptive statistics on normally distributed data, PCA is probably your best bet...) Suppose we wanted to understand variation in multiple leaf traits like LMA, Aarea, Nmass, and LL. How many axes of variation are there, and what traits correlate with each other? We can specify a one-sided formula (no dependent variable). If our data are not all in the same units we should probably look at correlations rather than covariances between them to assess the PCA axes. Our data are also already log-transformed to meet normality assumptions, but you may need to do this to your data before running a PCA. This is incredible - the first principal component explains 79.7% of the total variation in the dataset. So all these leaf traits are very correlated with each other. Let's see in which direction we find these correlations: So the first, most important axis is dominated by a tradeoff between high LMA, low Amass, low Nmass, and high LL. This is the global leaf economics spectrum right here!

R is capable of making publication-quality figures. The problem is that it often takes a lot of work to set up all the margins, sizes, fonts, labels, colors, etc. just as you would like. However it is worth the trouble to learn how!

In this tutorial we are going to cover just the most basic types of graphs. There are many packages out there which can make nearly any graph you like. To find out how, I recommend browsing the R graph gallery (http://addictedtor.free.fr/graphiques/thumbs.php) - will show you examples of graph types and demo code for how to generate them.

Let's first make a new window to plot in. If you don't specify this command before plotting, R will open a default-sized window. But often if you are making figures for publication (e.g. Nature requires figures that are exactly 3.503 inches wide, good luck doing that in Powerpoint) you will appreciate this extra control.

Now let's plot the relationship between LMA and LL: Of course we can have much more control over this plot - let's relabel the axes, add a title, and keep 1:1 aspect ratio, for starters We can also label the points based on any criteria we like. Here we will show the Angiosperms as black circles(pch=16) and gymnosperms as open triangles (pch=2). See if you can understand how this bit of code works! We can also change point color and size - how about we try angiosperms in red, gymnosperms in blue: We can also do more complex coloring - let's instead color the points with a rainbow color map that varies with each species' nitrogen content. Here we will define a rainbow color map (?rainbow) and take different index values from the cut function to get our plotting colors. As you can see leaves that have higher LMA and LL also have lower Nmass.

We should probably add a legend to the plot. Let's do this for the simpler gymno/angio plot:

We can also directly put text on the graph. Let's draw a graph with just the gymnosperm points, labeled by their Latin name in a small font and rotated 45 degrees: Now we can even overlay the angiosperm points on this graph - R lets us continually add to graphs. You might also want to investigate the lines and segments functions. segments is good at drawing error bars, for example.

We can also put fancier Greek symbols and subscripts/superscripts on the graph - look at ?plotmath for instructions. For example, let's subscript the Log10 in the axes:

We can also put multiple figures on one plot, and control margins and other such options. Try looking at ?par: If you ever want to control figure margins, font size, font face, etc... look in par first!

Another very useful type of plot is the boxplot, which will show the distribution of your data by a categorical variable. You may be tempted to show a bar plot, but bar plots don't show uncertainty ranges - and if you don't know an uncertainty range, you might not want to be making a graph yet anyway!

We can also easily add multiple combinations of categorical variables - here all the combinations of C3/C4 vs Angio/notAngio

What if we want to see lots of correlations between lots of variables at once? Fortunately this is a built-in function: With this function you can assess data distribution, correlations, and so on - very useful for initial exploratory work.

There is a more useful version of this function in the 'car' library: Here we also see histograms of data distributions and linear fits. You can explore the 'car' library for many more nice plotting functions. ggplot2 and lattice also have some nice plotting functionality if you want to explore them. I will instead show you how to generate this functionality using base functions.

Maybe we want to draw a line of best-fit through the data. Suppose we are interested in the relationship between Amass and Nmass:

We can now draw the fit on the graph with a double line width Strictly speaking we might not want to draw this line, because the regression is not valid at points beyond the range of the predictors. So let's draw a better line by first generating a sequence of evenly spaced points in the x-range, then predicting their y-values:

Generating x values and predicting y-values works for nonlinear models too, whereas abline only works for linear models. Suppose we try to fit a shifted-scaled square-root (k0*sqrt(x-k1)) model to the data. This uses the nonlinear least-squares command, which can also be useful for other data-fitting procedures. Now you can really see the curvature!

let's finally add a locally-smoothed line to the graph using LOESS. This can be useful for very complex datasets for which you don't have an a-priori model to fit.

We can also make histograms. Suppose we wanted to visualize the distribution of LMA. We can draw a histogram using the hist function. You can choose the number of 'breaks' or bins in the histogram, as well as do many fancier plots too - look at the help!

Histograms are rarely a good thing to display - the story they tell depends very heavily on the number of bins/breaks the user chooses. A better option is a kernel density estimate: You can choose smoothing bandwidths to generate a better looking density estimate, but the default options are usually quite good.
Let's superimpose the distributions of the angios and gymnos on the whole plot:

And now for a nice trick, let's graph the LMA distributions by each biome type - in three lines of code!

It's also often useful to visualize the results of a principal components analysis to assess pairwise axes of variation. First set dimension names so the biplot looks nice: Now do the pca and biplot (defaulting to showing PCA1 and PCA2)

Last up - phylogenetic analyses. If you haven't already loaded it for an earlier part of the tutorial, we need the 'ape' and 'picante' packages, which you can install if needed. Now plot species-mean LMA on the tree, using a 'fan' layout

We can also test for phylogenetic signal in the data using Blomberg's K statistic (you can also find tests for Pagel's lambda, but it seems to be less used these days): Don't worry about the warning. You can see here we get K=0.45, so this trait evolves significantly more rapidly (P = 0.001) than expected under a Brownian evolution model.

You can also calculate phylogenetic independent contrasts if you don't like using GLS with correlational structure (see above). Now we can look for a phylogeny-independent relationship between LMA and shrinkage using regression through the origin: A relationship still exists (but is weak).

Here are some suggestions on other packages you may want to investigate. These are all available on the Comprehensive R Archive Network (CRAN) (which can be accessed from the Package Installer within R).

• ape and picante have many more tools for phylogenetic analyses
• MASS has many nice utility functions and data manipulation/summary functions
• nlme and lme4 have tools for mixed modeling
• sem has structural equation modeling (path analysis) but a more modern tool is lavaan
• vegan is good for ecological analyses including community structure, rank-abundance, diversity metrics, etc.
• spatial / map and related packages have tools for plotting and mapping data, and working with spatial datasets (e.g. you can use gls with spatial autocorrelation)
• lattice and ggplot2 and car have many additional high-level plotting functions that will produce nice results
• rpart will do classification and regression trees
• igraph or sna will do network analyses (good for animal studies)

To really make the most of R, you will want to learn how to

• work with loops and control structures (if, for, while)
• become good at mass-processing data (tapply, sapply, apply, by)
• write your own functions (e.g. x <- function(a,b,c) to handle complex but often-repeated tasks
• extract parts of objects to programmatically do analyses or insert numbers onto graphs (e.g. get a R2 value from a lm summary object and add it to a plot)

Additional resources include

• http://www.rseek.org - a search engine devoted only to R and R mailing lists
• http://stat.ethz.ch/mailman/listinfo/r-help - a mailing list for help with R. Read the archives first - nearly everyone on this list is very grumpy
• http://cran.r-project.org/ - for finding extra packages
• http://cran.r-project.org/other-docs.html - other documentation and guides for R
• http://books.google.com/books/about/The_R_book.html?id=8D4HVx0apZQC - a hardcopy book that is a very extensive guide to many packages and common tasks