Working with the data on gliding snakes’ undulation rates we will calculate the mean, median, range, interquartile range, variance, and standard deviation. We will also construct some plots in order to see a few things in action.

```
snakeData = read.csv(url("http://whitlockschluter.zoology.ubc.ca/wp-content/data/chapter03/chap03e1GlidingSnakes.csv"))
head(snakeData)
```

```
## snake undulationRateHz
## 1 1 0.9
## 2 2 1.2
## 3 3 1.2
## 4 4 1.3
## 5 5 1.4
## 6 6 1.4
```

Drawing a histogram of these data will allow us to see how undulation rates are distributed.

`hist(snakeData$undulationRateHz, right = FALSE, las = 1, col = "gray", breaks = seq(0.8,2.2,by=0.2), xlab = "Undulation rate (Hz)", ylab = "Frequency", main = "Distribution of Undulation Rates (in Hz)", ylim=c(0,3))`

The code below constructs a similar plot albeit with the `ggplot2`

package

```
library(ggplot2)
ggplot(snakeData, aes(undulationRateHz)) + geom_histogram(binwidth=0.2) + xlab("Undulation rate (Hz)")
```

## Calculating Central Tendency & Variability Calculating the mean, median, variance, and standard deviation is pretty straightforward:

`mean(snakeData$undulationRate, na.rm=TRUE)`

`## [1] 1.375`

`median(snakeData$undulationRate, na.rm=TRUE)`

`## [1] 1.35`

`var(snakeData$undulationRate, na.rm=TRUE)`

`## [1] 0.105`

`sd(snakeData$undulationRate, na.rm=TRUE)`

`## [1] 0.324037`

**Note:** If you have missing values in your data-set then you should run `mean(snakeData$undulationRateHz, na.rm=TRUE)`

etc. for all commands. As a matter of fact R has several commands that will not run unless you specify `na.rm=TRUE`

. Further, if we wanted to we could have named each estimate, for example, `mean.ur = mean(snakeData$undulationRate)`

and so on.

We have to be a bit careful when calculating the quartiles (\(Q_1\) and \(Q_3\)) and the interquartile range. This is because there are nine ways of calculating these in R. We will use `type=5`

in this class since that is what is used by Whitlock & Schluter. The relevant commands are shown below:

`quantile(snakeData$undulationRateHz, probs=c(0.25, 0.75), type=5, na.rm=TRUE)`

```
## 25% 75%
## 1.2 1.5
```

`IQR(snakeData$undulationRateHz, type=5, na.rm=TRUE)`

`## [1] 0.3`

**Note:** \(IQR = Q_3 - Q_1 = 1.5 - 1.2 = 0.3\) . Specifying `probs=c(0.25, 0.75)`

specifies that we want \(Q_1\) and \(Q_3\). If you wanted the median then you’d have to specify a probs value of `0.50`

.

The Coefficient of Variation tells us something about how a variable is distributed about its mean. Given that it is a ratio, the higher is the coefficient of variation the more dispersed are the values of the variable. Why is it a useful measure if we have others to choose from? Because it is not influenced by the units of measurement. For example, if you are measuring heights in inches and calculate the variance and standard deviation, then convert height into centimeters and recompute the variance and standard deviation, the two sets of values will differ a good bit (see below):

```
height.in = c(65, 72, 77, 60, 58)
mean(height.in) # No missing values so na.rm=TRUE not specified
```

`## [1] 66.4`

`sd(height.in) # No missing values so na.rm=TRUE not specified`

`## [1] 8.018728`

```
height.cm = height.in * 2.54
mean(height.cm) # No missing values so na.rm=TRUE not specified
```

`## [1] 168.656`

`sd(height.cm) # No missing values so na.rm=TRUE not specified`

`## [1] 20.36757`

Now suddenly it seems as if there is more variability when height is measured in centimeters, but this is because of the units of measurement being used! If we now calculate the coefficient of variation for height in inches and then for height in centimeters we get identical estimates of how much variation there is in the heights of these individuals:

```
cov.in = (sd(height.in) / mean(height.in)) * 100
cov.in
```

`## [1] 12.0764`

```
cov.cm = (sd(height.cm) / mean(height.cm)) * 100
cov.cm
```

`## [1] 12.0764`

Now calculating the coefficient of variation for the gliding snakes data-set:

`100 * sd(snakeData$undulationRate)/mean(snakeData$undulationRate)`

`## [1] 23.56633`

For this exercise we’ll replicate the in-text example that uses the number of convictions data-set (see below):

```
convictionsFreq = read.csv(url("http://whitlockschluter.zoology.ubc.ca/wp-content/data/chapter03/chap03t1_2ConvictionsFreq.csv"))
head(convictionsFreq)
```

```
## convictions frequency
## 1 0 265
## 2 1 49
## 3 2 21
## 4 3 19
## 5 4 10
## 6 5 10
```

Notice the format of the data … it is given to us as a frequency table. Not a problem so long as we understand what we are looking at. In brief, 265 boys had no convictions, 49 had 1 conviction, and so on. We can expand the data into its long (or raw) form as shown below:

```
convictions = rep(convictionsFreq$convictions, convictionsFreq$frequency)
head(convictions)
```

`## [1] 0 0 0 0 0 0`

Now the calculations can commence:

`mean(convictions, na.rm=TRUE)`

`## [1] 1.126582`

`median(convictions, na.rm=TRUE)`

`## [1] 0`

`sd(convictions, na.rm=TRUE)`

`## [1] 2.456562`

`var(convictions, na.rm=TRUE)`

`## [1] 6.034698`

`quantile(convictions, probs=c(0.25, 0.75), type=5, na.rm=TRUE)`

```
## 25% 75%
## 0 1
```

`IQR(convictions, type=5, na.rm=TRUE)`

`## [1] 1`

`( sd(convictions) / mean(convictions) ) * 100`

`## [1] 218.0544`

`tab.c = table(convictions); tab.c`

```
## convictions
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
## 265 49 21 19 10 10 2 2 4 2 1 4 3 1 2
```

`barplot(tab.c, ylim=c(0,300), xlab="Number of Convictions", ylab="Frequency")`

## The Speeds of Male Tidarren spiders (pre- and post-amputation)

```
spiderData = read.csv(url("http://whitlockschluter.zoology.ubc.ca/wp-content/data/chapter03/chap03e2SpiderAmputation.csv"))
head(spiderData)
```

```
## spider speed treatment
## 1 1 1.25 before
## 2 2 2.94 before
## 3 3 2.38 before
## 4 4 3.09 before
## 5 5 3.41 before
## 6 6 3.00 before
```

**Note:** The two variables of interest are `speed`

and `treatment`

. We can start by constructing histograms of speeds before and after amputation. We will then look at a different graphic … the `boxplot`

, an extremely useful plot because it conveys a lot of information very simply. First the histograms

```
library(lattice)
histogram(~ speed | treatment, data=spiderData, layout=c(1,2), type=c("count"))
```

Now the box-plots:

`boxplot(speed ~ treatment, data=spiderData, horizontal=TRUE, ylab="Amputation Treatment", xlab="Speed (cm/s)", col="gray")`

For each group the boxplot maps a few statistics. First, the dark line inside each box is the Median for the group. Clearly the median speed is higher for the post-amputation group than for the pre-amputation group. The edges of the gray boxes reflect the first and third quartiles for each group; \(Q_1\) is to the left of the Median and \(Q_3\) is to the right of the Median. The “whiskers” go out to the minimum and maximum values.

If a distribution is symmetric then the box should be equally wide either side of the Median and the whiskers should be of similar length on either side of the box. If you see the plot looking stretched out on one side then that tells us it is skewed towards this side. For instance, see the stretch to the right for the “after”" group; these data are positively skewed (aka skewed right). For the “before” group we see the stretch to the left; this group’s speeds are negatively skewed (i.e. skewed left).

We also see a lone dot at the lower extreme of the plot for the “before” group. This is an **outlier** … any data point that is smaller than \(Q_1 - (1.5 \times IQR)\) or larger than \(Q_3 + (1.5 \times IQR)\). An outlier is essentially an extremely unusual value.

We will now extract and save the before-amputation data using the `subset`

switch. We will then repeat this exercise for the after-amputation data.

```
speedBefore = subset(spiderData, treatment == "before")
speedBefore
```

```
## spider speed treatment
## 1 1 1.25 before
## 2 2 2.94 before
## 3 3 2.38 before
## 4 4 3.09 before
## 5 5 3.41 before
## 6 6 3.00 before
## 7 7 2.31 before
## 8 8 2.93 before
## 9 9 2.98 before
## 10 10 3.55 before
## 11 11 2.84 before
## 12 12 1.64 before
## 13 13 3.22 before
## 14 14 2.87 before
## 15 15 2.37 before
## 16 16 1.91 before
```

```
speedAfter = subset(spiderData, treatment == "after")
speedAfter
```

```
## spider speed treatment
## 17 1 2.40 after
## 18 2 3.50 after
## 19 3 4.49 after
## 20 4 3.17 after
## 21 5 5.26 after
## 22 6 3.22 after
## 23 7 2.32 after
## 24 8 3.31 after
## 25 9 3.70 after
## 26 10 4.70 after
## 27 11 4.94 after
## 28 12 5.06 after
## 29 13 3.22 after
## 30 14 3.52 after
## 31 15 5.45 after
## 32 16 3.40 after
```

Now we can calculate the mean, median, variance, standard deviation, quartiles, and the interquartile range for each group.

`mean(speedBefore$speed, na.rm=TRUE); mean(speedAfter$speed, na.rm=TRUE)`

`## [1] 2.668125`

`## [1] 3.85375`

`median(speedBefore$speed, na.rm=TRUE); median(speedAfter$speed, na.rm=TRUE)`

`## [1] 2.9`

`## [1] 3.51`

`var(speedBefore$speed, na.rm=TRUE); var(speedAfter$speed, na.rm=TRUE)`

`## [1] 0.4115896`

`## [1] 0.9853183`

`sd(speedBefore$speed, na.rm=TRUE); sd(speedAfter$speed, na.rm=TRUE)`

`## [1] 0.6415525`

`## [1] 0.992632`

`quantile(speedBefore$speed, probs=c(0.25, 0.75), type=5, na.rm=TRUE)`

```
## 25% 75%
## 2.340 3.045
```

`quantile(speedAfter$speed, probs=c(0.25, 0.75), type=5, na.rm=TRUE)`

```
## 25% 75%
## 3.22 4.82
```

`IQR(speedBefore$speed, type=5, na.rm=TRUE)`

`## [1] 0.705`

`IQR(speedAfter$speed, type=5, na.rm=TRUE)`

`## [1] 1.6`

Just for completeness’ sake, multiply the IQR of each group by 1.5

`1.5 * IQR(speedBefore$speed, type=5, na.rm=TRUE)`

`## [1] 1.0575`

`1.5 * IQR(speedAfter$speed, type=5, na.rm=TRUE)`

`## [1] 2.4`

Now subtract 1.0575 from \(Q_1\) of the “before” group

`2.340 - 1.0575`

`## [1] 1.2825`

This yields a value of \(1.2825\) … Spider 1 in the “before” group has a speed of \(1.25\) and hence is flagged as an outlier in the boxplot.

```
par(mfrow=c(1,2))
plot( ecdf(speedBefore$speed), verticals = TRUE,
las = 1, main = "", do.points = FALSE,
ylab = "Cumulative relative frequency",
xlab = "Running speed before amputation (cm/s)")
plot( ecdf(speedAfter$speed), verticals = TRUE,
las = 1, main = "", do.points = FALSE,
ylab = "Cumulative relative frequency",
xlab = "Running speed after amputation (cm/s)" )
```

`dev.off()`

```
## null device
## 1
```

`par(mfrow=c(1,1))`

You can pick any point on the vertical axis, connect a straight line to the plotted line, then drop a perpendicular to the x-axis … and you’ll be able to say “so-and-so proportion of the spider had speeds below the x-value”. Out of sheer curiosity, pick \(0.5\) on the vertical axis of both groups and draw the straight line to the plot. Now read the approximate running speed that would result if you dropped a perpendicular to the x-axis. What this says is that 50% of spiders in the before group had speeds of at most 2.9 while 50% of spiders in the after group had speeds of at most 3.51.

```
par(mfrow=c(1,2))
plot( ecdf(speedBefore$speed), verticals = TRUE,
las = 1, main = "", do.points = FALSE,
ylab = "Cumulative relative frequency",
xlab = "Running speed before amputation (cm/s)")
abline(h=0.5)
plot( ecdf(speedAfter$speed), verticals = TRUE,
las = 1, main = "", do.points = FALSE,
ylab = "Cumulative relative frequency",
xlab = "Running speed after amputation (cm/s)" )
abline(h=0.5)
```

`dev.off()`

```
## null device
## 1
```

`par(mfrow=c(1,1))`

You don’t have to pick the \(0.5\) point; you could have picked \(0.9\) and said “90% of spiders in group so-and-so had speeds of so much or less”, so on and so forth.

The `psych`

package (install it before you use it) is very handy here to calculate a number of statistics. We will see it in use with the data-set on sticklebacks:

```
sticklebackData = read.csv(url("http://whitlockschluter.zoology.ubc.ca/wp-content/data/chapter03/chap03e3SticklebackPlates.csv"))
head(sticklebackData)
```

```
## id plates genotype
## 1 4-1 11 mm
## 2 4-2 63 Mm
## 3 4-4 22 Mm
## 4 4-5 10 Mm
## 5 4-10 14 mm
## 6 4-12 11 mm
```

```
library(psych)
describeBy(sticklebackData$plates, group=sticklebackData$genotype, na.rm=TRUE)
```

```
## group: mm
## vars n mean sd median trimmed mad min max range skew kurtosis se
## 1 1 88 11.67 3.57 11 11.42 2.97 6 37 31 4.2 26.92 0.38
## --------------------------------------------------------
## group: Mm
## vars n mean sd median trimmed mad min max range skew kurtosis
## 1 1 174 50.38 15.15 59 52.52 7.41 10 68 58 -1.06 -0.18
## se
## 1 1.15
## --------------------------------------------------------
## group: MM
## vars n mean sd median trimmed mad min max range skew kurtosis se
## 1 1 82 62.78 3.41 63 63.21 1.48 42 69 27 -3.55 17.62 0.38
```

Histograms are easily drawn as well:

`histogram(~ plates | genotype, data=sticklebackData, layout=c(1,3), type="count", xlab="Number of Plates", ylab="Frequency")`

Again, in a very marked fashion the \(MM\) group has the most number of plates on average.

Constructing a table of frequencies and proportions of the stickleback genotypes is easy as well.

```
tab.s1 = table(sticklebackData$genotype, dnn = "genotype")
tab.s1
```

```
## genotype
## mm Mm MM
## 88 174 82
```

```
tab.s2 = addmargins(tab.s1)
tab.s2
```

```
## genotype
## mm Mm MM Sum
## 88 174 82 344
```

```
sticklebackFreq = data.frame(tab.s2)
sticklebackFreq
```

```
## genotype Freq
## 1 mm 88
## 2 Mm 174
## 3 MM 82
## 4 Sum 344
```

```
tab.sprop1 = prop.table(tab.s1)
tab.sprop1
```

```
## genotype
## mm Mm MM
## 0.2558140 0.5058140 0.2383721
```

```
tab.sprop2 = addmargins(tab.sprop1)
sticklebackRelFreq = data.frame(tab.sprop2)
sticklebackRelFreq
```

```
## genotype Freq
## 1 mm 0.2558140
## 2 Mm 0.5058140
## 3 MM 0.2383721
## 4 Sum 1.0000000
```