---
title: "Introduction to the the R Project for Statistical Computing and the S language"
author: "D G Rossiter"
date: "`r Sys.Date()`"
output:
  html_document:
    fig_align: center
    fig_height: 4
    fig_width: 4
    number_section: yes
    theme: spacelab
    toc: yes
    toc_float: yes
  word_document:
    toc: yes
editor_options:
  chunk_output_type: inline
---

```{r global_options, include=FALSE}
knitr::opts_chunk$set(fig.align='center', warning=FALSE, message=FALSE)
options(show.signif.stars=FALSE)
```

This tutorial shows the main features of the S language as implemented in the R Environment for Statistical Computing (https://www.r-project.org/). 

R is fully described in: Venables, W N, D M Smith, and R Development Core Team. 2019. An Introduction to R; Notes on R: A Programming Environment for Data Analysis and Graphics. Version 3.6.2 (2019-12-12). Vienna: R Foundation for Statistical Computing. https://cran.r-project.org/manuals.html.

A simple introduction to R for statistical analysis is: Dalgaard, Peter. 2008. Introductory Statistics with R. Second. Springer. http://link.springer.com/book/10.1007%2F978-0-387-79054-1.

You should work through the code examples. At the end of each section is an *exercise*. Please create a new R Markdown document and complete these exercises; you may want to use some of the code from this document.

# R as an expression language

S code can be executed at the R console command line, or from a R script, or from an R Markdown code chunk, as a calculator. 

The simplest use of S is to evaluate mathematical *expressions*.

```{r}
2*pi/360
```

The ratio of a circle's radius squared to its area $\pi$ is a built-in constant `pi`. 

S has the usual scalar arithmetic operators +, -, *, /, ^ and some less-common ones like %% (modulus) and %/% (integer division). Expressions are evaluated in accordance with the usual operator precedence; parentheses may be used to change the precedence or make it explicit.

```{r}
3 / 2^2 + 2 * pi
((3 / 2)^2 + 2) * pi
```

Exercise 1: Write an expression to compute the number of seconds in a 365-day year, and execute the expression. 

# Assignment and workspace objects

Results of S expressions can be saved as *objects* in the *R workspace*. There are two (equivalent) assignment operators, the *left-hand side* is the name of the workspace object, which is assigned the value of the expression on the *right-hand side*.

Object names are case-sensitive and may contain special characters, but no spaces.

```{r}
rad.deg <- 2*pi/360
Rad2Deg = 2*pi/360
```

By default nothing is printed; but all of these give the same output:

```{r}
(rad.deg <- 2*pi/360)
rad.deg
print(rad.deg)
```

Workspace objects are listed with the `ls` "list" function and deleted with the `rm` "remove" function:

```{r}
ls()
rm(rad.deg)
ls()
```

The `character(0)` result means that the list of objects in the workspace is a zero-length character vector, i.e., null.

Exercise 2: Define a workspace object which contains the number of seconds in 365-day year, and display the results.

# Functions and methods

Most work in S is done with *functions* or *methods*. These are parts of expressions containing:

1. Method or function name; arguments to the function are between parentheses ( ) following the function name
2. Argument list:
  + Required
  + Optional, with defaults
  + positional and/or named

Functions _return_ the results.  These can be directly printed, assigned to workspace objects, or directly used as arguments to other functions or as parts of expressions.

For example, the `c` "catenate", "make a chain" function builds a vector; these can be saved as workspace objects. The arguments of this function are the items to be joined into one vector:

```{r}
(primes.lt.20 <- c(2, 3, 5, 7, 11, 13, 17))
(colours <- c("indigo", "violet", "red", "orange", "yellow", "green", "blue"))
```

Base R defined many common mathematical functions:

```{r}
sum(primes.lt.20); min(primes.lt.20); max(primes.lt.20); median(primes.lt.20)
sin(pi/4); cos(pi/4); tan(pi/4); atan(1)
log(128); log2(128); exp(4.85203)
```

Note that multiple expressions can be written on one line, separated by `;`.

Also many vector manipulation functions:

```{r}
rev(colours)  # reverse the order of elements in a vector
sort(colours) # sort in ascending order
sort(colours, decreasing=TRUE) # sort in decending order
```

Exercise 3: Find the function name for base-10 logarithms, and compute the base-10 logarithm of 10, 100, and 1000 (use the `??` function at the console to search).

The results of one function can be used as part of an expression:

```{r}
sqrt(2)
sqrt(2)/2
```

The results of one function can be used as an argument to another:

```{r}
asin(sqrt(2)/2)
rev(sort(colours)) # equivalent to `sort(colours, decreasing=TRUE)`
```

An example of a function with required and optional arguments is `rnorm` "random sample from a normal distribution". The only required argument is the number of items to be returned.

```{r}
rnorm(12)
```

Optional arguments are the parameters of this distribution, i.e., its mean and standard deviation:

```{r}
rnorm(n=12, mean=180)
rnorm(n=12, mean=180, sd=10)
```

How do we know all this? From the help system.

```{r}
help(rnorm)
```

(You can also access the help for this function by enterung `?rnorm` at the console).

This shows:

1. Title and package where found
2. Description
3. Usage: how to call the function
4. Arguments (what each one means, defaults)
5. Details of the algorithm
6. Value returned
7. Source of code
8. References to the statistical or numerical methods
9. See Also (related commands)
10. Examples of use and output

Exercise 4: What are the arguments of the `rbinom` (random numbers following the binomial distribution) function? Are any default or must all be specified? What is the value returned?

Exercise 5: Display the vector of the number of successes in 24 trials with probability of success `0.2` (20%), this simulation carried out 128 times. 


Since S is an expression language, the results of one function can be passed as an argument to another.

Here, the results of the `rnorm` function are passed to the `sort` function.

```{r}
sort(rnorm(n=12, mean=180, sd=10), decreasing=TRUE)
```

Note the use here of the _optional_ argument `decreasing` of the `sort` function.

# Including computations in the text

An important feature of R Markdown is the ability to include results of computations in the text, using the syntax `` `r ` `` followed by any R expression, which can include any workspace object defined in previous code chunks.

For example, suppose we compute the sample mean of a simulated normal distribution, rounded to two decimal places:

```{r}
print(my.mean <- round(mean(rnorm(n=12, mean=180, sd=10)), 2))
```

We then could report it like this: the sample mean is `r my.mean`. Notice how in the compiled document the in-line expression beginning with `` `r ` `` is replaced by its value from the computation.

Exercise 6: Summarize the  result of `rbinom` (previous exercise) with the `table` function.   What is the range of results, i.e., the minimum and maximum values? Which is the most likely result? For these, write text which includes the computed results.  This is necessary because the results change with each random sampling.

(Hint: convert the result of `table` to a `data.frame`; see `?table`; you may want to come back here after reading about `data.frame` in a later section.)

# Vectorized operations

A major feature of S is that most operations are _vectorized_, that is, they work element-wise on vectors. 

Above we built a vector with the `c` function:

```{r}
(primes.lt.20 <- c(2, 3, 5, 7, 11, 13, 17))
```

Elements can be selected from a vector with the `[]` matrix selection operator:

```{r}
primes.lt.20[1]
primes.lt.20[1:3]
primes.lt.20[c(1,3,5)]
```

For example, create a sequence of integers 1..10 and then add a normally-distributed error to it. These are added in parallel:

```{r}
(s <- seq(1, 10))
(r <- rnorm(10, 0, 0.1))
s[1] + r[1]
s + r
```

If one of the arguments is shorter than the other, it is _recycled_:

```{r}
primes.lt.20/2
```

Many functions operate element-wise on vectors:

```{r}
seq(0, 2*pi, by=pi/6)
round(sin(seq(0, 2*pi, by=pi/6)),4)
```

Others summarize a vector:

```{r}
sum(primes.lt.20)
```

Exercise 7: Create and display a vector representing latitudes in degrees from $0^\circ$ (equator) to $+90^\circ$ (north pole), in intervals of $5^\circ$. Compute and display their cosines -- recall, the trig functions in R expect arguments in radians. Find and display the maximum cosine.

# Packages

R is built from a set of packages, several of which are loaded by default in all R installations and required for normal operation. The `search` function shows the loaded packages and the order in which they are searched for function and object names:

```{r}
search()
```

There are thousands of contributed packages; these must first be installed on each R system from a CRAN repository^[https://cran.r-project.org/] using the `install.packages` function, and then loaded into the search space with the `library` or `require` function. Installation only needs to be done once per system, but loading one time each session. It is common practice to also load any packages on which the installed package depends.

Although packages are usually installed via the RStudio "Packages" tab, they can also be installed directly  with `install.packages`.

Here we load the `MASS` package that implements statistical methods from the textbook: Venables, W N, and B D Ripley. 2002. Modern Applied Statistics with S. Fourth edition. New York: Springer-Verlag.

```{r eval=FALSE}
install.packages("MASS", dependencies=TRUE)
```

Once installed, they must be _loaded_ into the workspace. Then their functions are available for use.

```{r}
library(MASS)
search()
help(package=MASS)
```

Exercise 8: Check if the `gstat` package is installed on your system. If not, install it. Load it into the workspace. Display its help and find the `variogram` function. What is its description?

# Classes

## Fundamental classes

All objects in S have a class. This controls what can be done with the object. For example, it does not make sense to use arithmetic operations on character strings.

Some basic types are `logical` (TRUE/FALSE), `numeric`, and `character`.

The `class` function returns the data type of an object:

```{r}
class(TRUE); class(316); class(3.1415926); class("blue")
```

Note both integers and real numbers are class `numeric`.

The type `list` can combine any objects, each with its own class:

```{r}
class(l <- list(1,FALSE,"red"))
```

Elements of a list are extracted with the `[[]]` operator:

```{r}
class(l[[2]])
```

Functions are also a class:

```{r}
class(sort)
```

The classes `logical`, `numeric`, and `character` are all vectors with one or more elements.

```{r}
class(c(TRUE, FALSE, TRUE))
class(primes.lt.20)
class(colours)
```

Exercise 9: Display the classes of the built-in constant `pi` and of the built-in constant `letters`.


## Derived classes

These are basic classes with some additional attributes.

Examples:

1. an `array` is a vector with a `dim` ``dimensions'' attribute
2. a `matrix` is a 2-D `array`

```{r}
primes.lt.20
dim(primes.lt.20)
my.a <- array(primes.lt.20)
class(my.a)
dim(my.a)
(my.matrix <- matrix(data=rep(primes.lt.20,2), nrow=length(primes.lt.20), byrow=FALSE))
class(my.matrix); dim(my.matrix)
my.matrix[,2] <- my.matrix[,2]^2
my.matrix
```

Certain operations and functions are restricted to work on certain classes. For example, a `matrix` may be transposed with the `t` "transpose" function:

```{r}
t(my.matrix)
```

Some functions automtically "promote" the class of their arguments, if possible:

```{r}
class(t(primes.lt.20))
```

## Classes defined by functions

Functions can define their own classes. For example, the `rlm` "robust linear models" function of the `MASS` package returns an object with several classes defined by both `rlm` and the function `lm` "linear models" on which it depends:

```{r}
data(birthwt, package="MASS")
class(model <- rlm(low ~ ., birthwt))
```

Model formulas and built-in datasets will be explained below.

Exercise 10: What is the class of the object returned by the `variogram` function? (Hint: see the heading "Value" in the help text.)

# Example datasets

R has many example datasets in the `datasets` package, which is loaded by default in all R installations. There are also datasets in most contributed packages. These are intended to show the features of various functions,

(You can also load your own data, of course; see below.)

The `data` function with the optional `package` argument lists the datasets available in the named packages:

```{r}
data(package="datasets")
data(package="ggplot2")
```

These can be loaded into the workspace with the `data` function. If the package which provides the dataset is loaded, there is no need to specify it. 


```{r}
data(diamonds, package="ggplot2")
data("birthwt") # in MASS, this was loaded above, so no need to specify
```

A dataset (actually, any S object) can be summarized with the `summary` function, which gives output appropriate to the object class. An object's internal structure is revealed with the `str` "structure" function:

```{r}
is(diamonds)
summary(diamonds)
str(diamonds)
```

Exercise 11: List the datasets in the `gstat` package.

Exercise 12: Load, summarize, and show the structure of the `oxford` dataset.

# Data frames

The "data frame" is the most common class for statistical and database work.
An object of class `data.frame` is a `matrix` with column (field) `colnames`  and (optional) `row.names`.

Rows are generally observations or individuals (database "cases") and columns are attributes (database "fields").

```{r}
(m.d <- as.data.frame(my.matrix))
class(m.d)
str(m.d)
colnames(m.d)
colnames(m.d) <- c("primes", "primes.sq")
rownames(m.d)
rownames(m.d) <- letters[1:length(rownames(m.d))]
rownames(m.d)
str(m.d)
```

There are a variety of ways to extract parts of a dataframe. We illustrate this with the `trees` sample dataset.

The most common way to extract a field as a vector is with the `$` operator:

```{r}
data(trees)
str(trees)
# fields
trees$Volume 
```

However, since a dataframe is just a matrix with named rows and columns, matrix selection operator `[]` also can be used:

```{r}
trees[,1]
trees[,"Volume"]
# cases (records, tuples)
trees[1,]
# subsets of cases and selected fields
trees[1,1]
trees[1:3, c(1,3)]
```

Exercise 13: load the `women` sample dataset. How many observations (cases) and how many attributes (fields) for each case? What are the column (field) and row names? What is the height of the first-listed woman?

## Factors

R makes a strong distinction between _continuous numeric_ variables and _categorical_ variables. The latter are called R `factors` and may be ordered (a natural ordering of class levels) or unordered.

```{r}
is.factor(diamonds$color)
is.ordered(diamonds$color)
levels(diamonds$color)
is.factor(diamonds$carat)
```

The identification as a factor or not has major implications for statistical models (see below).

Exercise 14: List the factors in the `oxford` dataset.

# Missing values

S has a special `NA` value for missing values. For example, in a data frame there may be no value for a certain field for one or more cases (observations). The missing value is handled correctly by functions.

Suppose there was no volume measurement for the first tree in the `trees` dataset:

```{r}
trees[1,]
trees[1,"Volume"] <- NA
trees[1,]
summary(trees$Volume)
```

This is different from "not a number", `NaN`, which is the result of an illegal mathematical operation, and `Inf`, which is the result of dividing by zero. If these are used in further mathematical operations, the error wil propagate.

```{r}
(x <- sqrt(-1))
x^2
(x <- 100/0)
x+1
```


# Logical expressions

The usual logical S operators `<`, `>`, `<=`, `>=`, `==`, `!=` return the logical values  `TRUE` and `FALSE`. The `which` function then selects the indices of the `TRUE` elements of a vector. 

```{r}
trees$Height > 80
(ix <- which(trees$Height > 80))
dim(trees)
(tall.trees <- trees[ix, ])
dim(tall.trees)
```

Exercise 15: Identify the thin trees, defined as those with height/girth ratio more than 1 s.d. above the mean. You will have to define a new field in the dataframe with this ratio, and then use the `mean` and `sd` summary functions, along with a logical expression.

# Combination: functions, logical expressions, simulation

Here is an example of using logical expressions, along with simulation: "“What are the chances that three people on an elevator with 13 buttons press three consecutive floors?"

```{r}
mean(replicate(2^12, all(diff(sort(sample.int(n=13, size=3, replace=TRUE))) == 1)))
```

Let's break this down "inside-out", but only using a small number of simulations:

```{r}
sample.int(n=13, size=3, replace=TRUE) # three floors selected, can be repeats
sort(sample.int(n=13, size=3, replace=TRUE)) # same, but in order
diff(sort(sample.int(n=13, size=3, replace=TRUE))) # delta between each choice
diff(sort(sample.int(n=13, size=3, replace=TRUE))) == 1 # are they consecutive 1 by 1?
all(diff(sort(sample.int(n=13, size=3, replace=TRUE))) == 1) # are they all consecutive?
replicate(12, all(diff(sort(sample.int(n=13, size=3, replace=TRUE))) == 1)) # repeat the test many times
mean(replicate(2^12, all(diff(sort(sample.int(n=13, size=3, replace=TRUE))) == 1))) # what is the average of the replicated tests?
```

Challenge: run this several times and see how variable is the simulation result.

# Graphics

R has several graphics systems, beginning with _base graphics_ which are installed by default. The `plot` method is generic, and changes its behaviour according to its arguments.

The simplest plot is the _scatterplot_ of two variables, which uses the `~` formula operator (see 'Statistical Models', below) to symbolize "y-axis" vs. "x-axis".


```{r}
plot(Height ~ Girth, data=trees)
# plot(trees$Height ~ trees$Girth)
```

Base graphics also has functions for histograms, boxplots, dot charts, and conditioning plots. See `help(package=graphics)`.

```{r out.width=c('30%', '30%', '30%'), fig.show='hold'}
hist(trees$Volume); rug(trees$Volume)
boxplot(diamonds$carat)
dotchart(trees$Volume)
```

Another graphics system, `ggplot2`, is becoming ever more popular. This views graphics as being built in layers, each with functions, starting with `ggplot2` to initialize the graphic, then each layer added with a `+` operator.

For example, here is a histogram of the tree volumes with a "rug" plot to show the actual values:

```{r}
library(ggplot2)
g.h <- ggplot(data=trees) +
        geom_histogram(mapping = aes(x=Volume), binwidth = 5,
                        colour="blue") +
        geom_rug(mapping = aes(x=Volume))
print(g.h)
```

See [here](https://ggplot2.tidyverse.org/index.html) for more on `ggplot2`.

Exercise 16: Display a histogram of the diamond prices in the `diamonds` dataset.

# Statistical models

S specifies statistical models in symbolic form with *model formulae*, which are arguments to many statistical methods, for example `lm` (linear models, in the base packages) and `rlm` (robust linear models, in the `MASS` package).

These have a:

* _Left-hand side_: (mathematically) dependent variable(s)
* _Formula operator_ `~`, read as "depends on", "is modelled by"
* _Right-hand side_: (mathematically) indepdendent variable(s)

The simplest use is in simple linear regression with the `lm` function. For example, to model the linear relation of tree `Volume` from its `Height`, the formula is `Volume ~ Height`.

```{r}
model.vh <- lm(Volume ~ Height, data=trees)
# equivalent to: model <- lm(trees$Volume ~ trees$Height)
summary(model.vh)
```

Exercise 17: Write a model to predict tree height from tree girth. How much of the height can be predicted from the girth? 

The right-hand side can have _formula operators_, which look like arithmetic opertors but are interpreted in terms of the model structure:

* `+` additive effect
* `*` interaction effect
* `/` nested model
* `-` remove terms from the model

For arithmetic within model formulas, the `I` "identity" function must be used.

To illustrate the use of opertors, here is a model of tree volume `Volume` modelled by the square root of tree height `sqrt(Height)`, interacting with the tree radius `I(Girth/(2*pi))`, with no intercept (`-1`):

```{r}
model <- lm(Volume ~ sqrt(Height)*I(Girth/(2*pi)) - 1, data=trees)
summary(model)
```

_This model may or may not make sense physically!_

Exercise 18: Write a model to predict tree volume as a linear function of tree height and tree girth, with no interaction.

Variables in statistical models may be categorical, i.e., R `factors`; these are specified the same way in models but interpreted differently. For example in the linear model `lm` function, a factor is converted into _contrasts_, typically so-called "dummy" variables:

```{r}
summary(m.c <- lm(price ~ carat, data=diamonds))
head(model.matrix(m.c))
summary(m.f <- lm(price ~ clarity - 1, data=diamonds))
head(model.matrix(m.f))
```

# Control structures

Within an expression the `ifelse` function can be used to select two outcomes, based on a logical expression.

For example, to select a plotting colour based on whether one or the other of two variables in a scatterplot is greater:

```{r fig.height=4, fig.width=4}
# two vectors of 100 random uniform variates
n <- 128
x <- runif(n); y <- runif(n)
plot(y ~ x, asp=1, col=ifelse(y > x, "red", "green"), pch=20)
abline(0, 1, lty=2); grid()
```

S has ALGOL-like _control structures_:

* `if` ... `else`
* `for`
* `while`, `repeat`
* `break`, `next`

Note that _vectorized functions_ or the `apply` family of functions can often be used where other programming languages would use `for` loops.

Here is an example of the `while` control structure. For some simulation we want to draw a sample from the normal distribution but make sure that there is an extreme value, so we repeat the sampling until we get what we want:

```{r}
n <- 1
while (max(abs(sample <- rnorm(80))) < 3) {
  print("No extreme")
  n <- n + 1
}
print(paste("Sample required",n,"draws"))
```

The `for` loop works as in other languages:

```{r}
v <- vector(mode="numeric", length=11)
for (power in 0:10) {
  v[power+1] <- 2^power
}
print(v)
```

But of course this is much easier with a vectorized operation:

```{r}
0:10
(v <- 2^(0:10))
```

## The `apply` family of functions

These apply a function across either the rows or columns of a matrix, including data frames.

For example, to find the maximum value of all the variables in a data frame we can use the `sapply` "apply a function and simplify the result" function:

```{r}
sapply(trees, FUN=max)
```

This applies the `max` built-in function along the fields of the data frame.

The `by` function performs a function on subsets defined by a factor (categorical varibable):

```{r}
by(diamonds[,c(1,3:5)], diamonds$cut, summary)
``` 

# User-defined functions

A function is defined with the `function` function (!), specifying the arguments, body of the function, and the return values:

For example, here is a function to compute the geometric mean of a vector $\overline{v_h} = \left[ \prod_{i=1\ldots n} \; v_i  \right]^{1/n}$:

```{r}
hm <- function(v) {
  return(exp(sum(log(v))/length(v)))
}
class(hm)
hm(1:99); mean(1:99)
```

Of course, a function should check for valid inputs. This example shows the use of the `if`, `else if`, `else` control structure, as well as the functions `is.numeric`, `any`, and `length`, and the `return` function to return a value from a function.

```{r}
g.mean <- function(v) { 
    if (!is.numeric(v)) { 
      print("Argument must be numeric"); return(NULL) 
      } 
    else if (any(v <= 0)) { 
      print("All elements must be positive"); return(NULL) 
      } 
    else return(exp(sum(log(v))/length(v))) 
} 
g.mean(letters) 
g.mean(c(-1, -2, 1, 2)) 
g.mean(1:99) 
```

Exercise 19: Write a function to restrict the values of a vector to the range $0 \ldots 1$. Any values $< 0$ should be replaced with $0$, and any values $>1$ should be replaced with $1$.  Test the function on a vector with elements from $-1.2$ to $+1.2$ in increments of $0.1$ -- see the `seq` "sequence" function.

# Import/Export

R can exchange data in any conceivable format.

Reference: "R Data Import/Export", an R manual installed with R; available under the Help menu as `R Help`.`

The most common interchange format for "flat" files is the Comma-Separated Values (CSV). Here we first download a sample CSV file and then import to R:

```{r}
download.file(
  url="https://data.ny.gov/api/views/vpv5-zd4k/rows.csv?accessType=DOWNLOAD",
  destfile="enplanements.csv")
file.show("enplanements.csv")  # examine in another window
airports <- read.csv("enplanements.csv")
str(airports)
```

This file is described [here](https://data.ny.gov/Transportation/Annual-Enplanements-in-NYS-Airports-Beginning-1997/vpv5-zd4k)

The more general `read.table` function can read tables in many formats.

# The Tidyverse

In the past few years "an opinionated collection of R packages designed for data science" called the `tidyverse` has been developed by Hadley Wickham and colleagues.
The `ggplot2` package is part of this; other important packages are `dplyr` for data manipulation,`readr` for data import, and `stringr` for string manipulation.

The complete tidyverse can be installed with:

```{r eval=FALSE}
install.packages("tidyverse")
```

The `dplyr` package provides a flexible grammar of data manipulation,
including the "pipe" operator `%>%`, which works from left to right, unlike the conventional `<-` assignment operator.

```{r}
library(dplyr)
```

The typical use is to do some data manipulation by filtering, as well as to summarize.
For example, to subset the airports to just the upstate hubs and summarize the enplanements by year:

```{r}
levels(airports$Airport.Group)
airports %>% filter(Airport.Group=="UPSTATE HUB") %>% 
  group_by(Year) %>% 
  summarize(Enplanements=sum(Enplanements))
```

Another typical use is to define new variables, possibly as combinations of others.
For example, load the Meuse River geochemical dataset, select only the most frequently-flooded areas, rename the columns for the heavy metals, compute new columns with their base-10 logarithms, compute a ratio, and finally save as a new object:

```{r}
library(sp)
data(meuse)
names(meuse)
meuse %>% 
  filter(ffreq==1) %>%
  select(-ffreq) %>%  # column is no longer  needed, it is constant
  rename(Cd=cadmium, Cu=copper, Pb=lead, Zn=zinc) %>%
  mutate(lCd=log10(Cd), lCu=log10(Cu), lPb=log10(Pb), lZn=log10(Zn)) %>%
  mutate(ratio.lCu.lPb=lCu/lPb) ->
  meuse.new
meuse.new[1:6,]
```

Note the use here of the right-hand assignment operator `->` to save the result into a new object (we could also overwrite the old object `meuse` if we had no use for its original form).

Bonus Exercise : Use `tidyverse` functions and pipes on the `trees` dataset, to select the trees (use the `filter` function) with a volume greater than the median volume (use the `median` function), compute the ratio of girth to height as a new variable (use the `mutate` function), and sort by this (use the `arrange` function) from thin to thick trees.

```{r}
names(trees)
mv <- median(trees$Volume)
trees %>% 
  filter(Volume > median(Volume)) %>%
  mutate(thickness=round(Girth/Height,3)) %>%
  arrange(thickness)
```