Using R for data analysis

at Permafrost Research Group
November 6, 2013

Ken Butler
Dept of Computer and Mathematical Sciences
School of Graduate Studies, DPES

(http://www.utsc.utoronto.ca/%7ebutler/prg/)

(http://www.utsc.utoronto.ca/%7ebutler/prg/prg_pres.html)

Why not use a spreadsheet?

This: spreadsheet addiction
Excellent tool for what they are good at:
displaying data
simple calculations and plots
Values and formulas: sorting cells with formulas (especially =rand()) problematic.
Defining a formula and adding data to spreadsheet: did you update the formula?
Easy to get an answer. Can be hard to get the right answer.
Command history: what did you do, and in what order?
Excel-specific problems (eg. treatment of blank cells)
This: most dangerous software

So what to do? Use R!

Use software designed for doing statistics.
R free and open source.
Has huge user community that wants to help.
Makes superb graphics.
Well-tested: if a bug is found, it is fixed.
Downsides:
need to learn how to use it
not as obviously straightforward as a spreadsheet
"R way to do things", mindset.
command prompt

Using R

I recommend front end R Studio. Download and install R and R Studio. Run R Studio.
R
R Studio
Basic operation: read in a few values, calculate mean and standard deviation, draw a histogram and boxplot:

x = c(10, 11, 11, 13, 14, 15, 15, 16, 17, 19, 24, 33)
mean(x)

[1] 16.5

sd(x)

[1] 6.474

hist(x)

boxplot(x)

Reading in real data

Data in a spreadsheet: save as .csv
Slawek Kowal's Hudson Bay ice breakup data: 36 locations in Hudson Bay, date (in days since Jan 1) that ice cover first < 50%, for years starting 1971.
Locations across top, years down side.
I omitted 2001 and 2002 (missing values).
Read in data and look at it ("data frame"):

breakup = read.csv("breakup.csv", header = T)
breakup

      X  X1  X2  X3  X4  X5  X6  X7  X8  X9 X10 X11 X12 X13 X14 X15 X16
1  1971 187 187 201 215 201 201 208 194 194 201 187 208 194 173 187 208
2  1972 182 189 205 210 217 217 217 217 203 196 217 210 217 203 196 203
3  1973 205 191 191 212 198 205 226 198 191 198 205 212 198 177 219 212
4  1974 197 197 190 197 204 218 211 204 204 204 211 204 204 169 218 211
5  1975 189 182 168 210 217 210 196 196 196 182 203 203 189 168 210 203
6  1976 194 194 187 222 208 194 208 194 194 194 201 208 201 173 194 201
7  1977 165 179 186 200 207 200 207 207 193 165 186 200 200 193 193 200
8  1978 163 184 177 191 198 184 205 198 177 205 219 219 212 191 205 219
9  1979 176 162 204 204 190 204 204 204 211 155 197 204 204 176 197 197
10 1980 160 160 167 209 209 209 209 209 153 202 209 188 209 153 188 188
11 1981 201 166 194 187 208 215 208 208 208 187 208 208 208 215 187 208
12 1982 214 165 200 207 200 214 221 207 193 207 207 207 200 179 221 207
13 1983 201 180 194 201 215 201 222 215 187 180 222 215 215 159 201 215
14 1984 183 176 183 187 208 215 208 208 208 162 208 208 208 162 187 208
15 1985 193 186 186 207 214 221 207 209 186 214 228 200 179 136 214 207
16 1986 181 167 174 202 209 188 209 209 195 195 195 209 209 181 209 209
17 1987 166 173 166 208 194 180 208 208 215 180 208 208 208 208 208 208
18 1988 157 171 178 206 199 185 199 199 192 171 213 199 192 157 192 213
19 1989 184 163 170 198 219 198 212 191 205 191 205 198 219 219 184 198
20 1990 162 169 183 169 197 204 218 197 218 169 190 197 211 197 190 190
21 1991 168 161 168 175 196 203 175 203 189 161 203 203 196 175 189 196
22 1992 180 194 173 222 222 236 236 229 194 208 229 236 215 201 236 229
23 1993 193 172 165 193 186 200 200 193 193 221 186 200 193 179 186 200
24 1994 185 171 192 185 206 206 227 213 213 185 213 220 213 234 185 206
25 1995 163 177 177 212 191 205 212 191 198 170 212 212 198 184 191 191
26 1996 168 168 189 189 203 189 217 203 210 175 224 217 203 189 203 203
27 1997 178 178 185 199 199 199 227 206 199 164 220 206 199 171 192 199
28 1998 153 174 153 195 188 188 188 181 174 167 181 188 181 174 188 181
29 1999 166 166 166 180 173 180 173 173 173 122 166 166 173 153 166 166
30 2000 192 178 136 192 227 227 213 213 199 213 213 206 192 167 199 206
31 2003 182 153 182 182 196 189 182 203 210 167 189 189 203 210 182 189
32 2004 223 137 202 237 223 237 223 223 230 202 216 216 216 181 202 209
33 2005 172 179 153 153 179 186 186 172 172 179 200 193 200 172 186 193
34 2006 185 178 192 206 185 185 206 185 185 122 185 178 178 164 185 185
35 2007 177 163 177 191 191 191 184 191 191 121 198 191 198 149 177 191
36 2008 183 190 148 197 207 204 218 197 204 113 218 204 190 134 197 204
37 2009 188 146 188 174 209 209 216 209 188 174 223 209 202 188 216 209
38 2010 152 138 117 208 187 180 180 180 180 124 180 180 180 131 187 187
39 2011 172 165 172 179 186 179 193 193 172 165 200 193 186 172 165 179
   X17 X18 X19 X20 X21 X22 X23 X24 X25 X26 X27 X28 X29 X30 X31 X32 X33 X34
1  208 194 173 208 208 208 187 173 187 201 201 187 173 180 180 180 180 180
2  217 217 203 210 210 217 217 196 210 196 210 217 203 189 189 196 217 175
3  212 198 184 212 205 212 205 177 191 198 198 198 184 191 184 177 177 170
4  211 176 169 204 204 204 176 169 204 197 204 176 169 183 197 211 169 169
5  203 168 168 196 203 196 182 168 196 196 196 168 168 196 210 196 182 217
6  201 194 166 201 201 208 208 171 187 201 201 208 166 194 173 187 187 173
7  200 200 163 193 186 193 200 193 171 186 171 200 207 158 158 158 200 158
8  219 205 198 219 205 205 219 212 177 205 205 219 212 163 198 205 205 163
9  204 204 169 197 197 197 204 155 190 197 183 190 155 169 183 183 183 183
10 195 181 136 188 188 202 195 181 181 188 209 195 181 174 181 202 167 188
11 208 208 187 201 194 208 201 201 173 180 187 201 194 136 173 166 194 136
12 200 186 179 200 200 200 186 186 193 193 200 193 193 193 193 200 193 186
13 215 208 173 201 208 208 208 159 208 201 208 194 173 194 208 208 194 194
14 208 176 152 201 208 194 187 136 194 208 201 194 136 194 194 187 162 208
15 200 200 165 200 200 200 179 165 200 207 193 179 165 193 193 193 179 186
16 209 209 188 204 188 209 202 209 167 167 202 202 202 160 160 181 181 160
17 208 208 208 208 208 208 208 208 159 208 208 194 208 159 201 194 208 159
18 199 192 164 227 199 206 185 192 199 206 213 185 164 206 199 199 157 171
19 198 184 177 184 198 198 198 198 156 198 198 191 198 156 170 191 191 156
20 197 211 204 162 199 197 204 204 162 162 176 204 211 142 162 169 197 162
21 196 196 196 189 203 203 182 196 189 189 203 196 189 189 182 203 203 154
22 236 222 201 201 222 229 222 215 201 201 215 215 173 166 166 201 201 166
23 193 172 142 186 200 186 186 165 186 186 186 193 142 193 186 186 193 193
24 206 213 206 185 199 213 206 206 178 185 199 199 169 122 178 164 157 157
25 191 191 170 163 184 198 191 184 163 184 184 177 177 142 170 170 170 170
26 203 203 189 189 189 203 203 189 142 189 189 203 189 142 154 175 189 122
27 199 199 164 178 185 199 199 164 171 185 192 199 164 164 164 164 164 164
28 181 174 153 167 188 181 174 167 167 167 188 211 181 153 153 181 174 153
29 166 173 153 173 173 180 180 187 180 180 180 180 166 187 187 173 153 187
30 199 199 153 185 185 192 185 136 185 185 178 171 153 167 192 178 167 185
31 203 203 210 153 167 189 203 203 153 167 175 196 203 136 167 175 175 153
32 223 209 202 188 216 223 216 188 154 181 202 209 188 154 167 195 160 154
33 193 186 179 200 179 193 186 186 179 179 193 193 172 179 179 179 179 179
34 178 178 122 178 178 178 178 171 178 178 178 171 150 171 171 178 150 150
35 191 191 156 128 184 191 191 163 177 177 198 198 184 142 170 184 142 149
36 197 190 169 162 197 190 176 127 162 176 176 183 176 155 176 183 155 155
37 209 202 181 167 209 209 202 188 146 202 202 195 181 146 153 181 146 146
38 180 187 145 180 187 187 180 131 180 180 180 159 131 138 173 180 131 138
39 200 186 172 165 186 200 200 186 137 179 200 200 200 137 179 179 151 137
   X35 X36
1  180 173
2  182 224
3  170 191
4  204 169
5  182 168
6  173 173
7  158 186
8  184 171
9  183 190
10 167 181
11 152 187
12 186 193
13 187 173
14 162 194
15 193 193
16 181 181
17 166 194
18 164 199
19 170 191
20 169 204
21 189 175
22 166 166
23 179 193
24 122 178
25 163 184
26 142 161
27 164 171
28 153 174
29 166 173
30 153 185
31 136 175
32 154 202
33 122 193
34 136 192
35 170 156
36 155 183
37 146 167
38 131 145
39 144 186

Column X is year, X1-X36 locations 1-36.

Accessing rows and columns

By number or name (columns).
Bits of data frame accessed as breakup[x,y], x row(s), y column(s).
1980 breakup (row 10) for location 3 (column 4, since year first):

breakup[10, 4]

[1] 167

All the 1980 figures (row 10):

breakup[10, ]

      X  X1  X2  X3  X4  X5  X6  X7  X8  X9 X10 X11 X12 X13 X14 X15 X16
10 1980 160 160 167 209 209 209 209 209 153 202 209 188 209 153 188 188
   X17 X18 X19 X20 X21 X22 X23 X24 X25 X26 X27 X28 X29 X30 X31 X32 X33 X34
10 195 181 136 188 188 202 195 181 181 188 209 195 181 174 181 202 167 188
   X35 X36
10 167 181

All the location 3 figures (column 4):

breakup[, 4]

 [1] 201 205 191 190 168 187 186 177 204 167 194 200 194 183 186 174 166
[18] 178 170 183 168 173 165 192 177 189 185 153 166 136 182 202 153 192
[35] 177 148 188 117 172

Statistics for rows and columns

For just one row/col, feed appropriate column into mean, sd etc:

mean(breakup[, 4])

[1] 177.9

For all rows or cols, use apply:

apply(breakup, 2, mean)

     X     X1     X2     X3     X4     X5     X6     X7     X8     X9 
1990.5  180.5  172.5  177.9  197.7  201.7  201.4  206.6  200.8  194.8 
   X10    X11    X12    X13    X14    X15    X16    X17    X18    X19 
 177.2  204.5  202.9  199.8  177.4  195.7  200.7  201.4  194.7  174.1 
   X20    X21    X22    X23    X24    X25    X26    X27    X28    X29 
 188.5  195.9  200.4  195.2  179.6  177.8  188.8  194.4  193.4  178.2 
   X30    X31    X32    X33    X34    X35    X36 
 167.0  178.8  184.9  176.5  166.8  164.2  181.9

These are mean breakup dates for all columns (for each location, averaged over years).

Or, for each year, averaged over locations:

apply(breakup, 1, mean)

 [1] 239.9 253.1 244.5 242.8 239.4 241.0 235.1 246.3 237.4 234.6 240.1
[12] 245.5 246.7 237.8 242.0 240.2 244.7 238.8 239.0 236.5 237.1 254.5
[23] 235.1 240.3 232.6 235.4 234.4 224.1 220.5 235.3 231.2 247.7 229.9
[34] 223.6 225.1 228.6 236.1 214.7 227.2

Columns by name

Long way: all the breakup dates for location 25:

breakup$X25

 [1] 187 210 191 204 196 187 171 177 190 181 173 193 208 194 200 167 159
[18] 199 156 162 189 201 186 178 163 142 171 167 180 185 153 154 179 178
[35] 177 162 146 180 137

Short way: attach data frame, then don't need to use the $. (When done with data frame, detach it.):

attach(breakup)
X25

 [1] 187 210 191 204 196 187 171 177 190 181 173 193 208 194 200 167 159
[18] 199 156 162 189 201 186 178 163 142 171 167 180 185 153 154 179 178
[35] 177 162 146 180 137

mean(X25)

[1] 177.8

Time plots

To plot a collection of values ("vector") against order:

plot(X25)

Joining points by lines

plot(X25, type = "b")

Scatterplots: plotting by year (`X`)

plot(X25 ~ X, type = "b")

Titles

plot(X25 ~ X, type = "b", ylab = "Location 25", xlab = "Year", main = "Time plot for location 25")

Regression line

trend = lm(X25 ~ X)
summary(trend)


Call:
lm(formula = X25 ~ X)

Residuals:
    Min      1Q  Median      3Q     Max 
-30.858 -12.489  -0.217  12.802  24.595 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1942.775    404.227    4.81  2.6e-05 ***
X             -0.887      0.203   -4.37  9.8e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 15.1 on 37 degrees of freedom
Multiple R-squared:  0.34,  Adjusted R-squared:  0.322 
F-statistic: 19.1 on 1 and 37 DF,  p-value: 9.78e-05

Scatterplot with regression line

plot(X25 ~ X, type = "b")
abline(trend)

Scatterplot with lowess curve

Lowess curve is nonparametric trend curve, not affected by outliers
Gives sense of whether trend is straight or not.

plot(X25 ~ X, type = "b")
lines(lowess(X25 ~ X))

Lowess curve for location 33

plot(X33 ~ X, type = "b")
lines(lowess(X33 ~ X))

Colours, line types and such (example)

plot(X33 ~ X, type = "b", lty = "dotted", col = "blue")
lines(lowess(X33 ~ X), lty = "dashed", col = "red")

Drawing maps

Data often spatial, need to plot things on a map.
R offers 2 ways:
blank map
Google map
Plot points and lines on map
Annotate map (a bit)

Blank map

Idea: get world map, take relevant piece of it
Uses package rworldmap. If not used before, need to do this (once) first:

install.packages("rworldmap")

Then, load it (once per R session):

library(rworldmap)

Loading required package: sp
### Welcome to rworldmap ###
For a short introduction type :      vignette('rworldmap')

Get a world map. Low-resolution actually good enough.

newmap = getMap(resolution = "low")

Plot with x (longitude, east=positive) and y (latitude, north=positive) limits.

Hudson Bay

plot(newmap, xlim = c(-90, -80), ylim = c(50, 65))

New Zealand

plot(newmap, xlim = c(166, 180), ylim = c(-47, -33))

Put some cities on NZ map

Data file nz.txt

nz = read.table("nz.txt", header = T)
nz

             city  long    lat
1        Auckland 174.7 -36.83
2        Hamilton 175.3 -37.79
3      Wellington 174.8 -41.29
4          Napier 176.9 -39.52
5     NewPlymouth 174.1 -39.08
6 PalmerstonNorth 175.7 -40.40
7    Christchurch 172.7 -43.52
8         Dunedin 170.5 -45.87
9    Invercargill 168.4 -46.41

Adding cities to map

Add city locations and names to plot with text:

plot(newmap, xlim = c(166, 180), ylim = c(-47, -33))
points(nz$long, nz$lat, col = "red")
text(nz$long, nz$lat, nz$city, col = "blue", pos = 4)

Gagnon and Gough

36 locations within Hudson Bay, in file locations.txt:

hb = read.table("locations.txt", header = T)
head(hb)

   lat long
1 54.0   81
2 54.5   82
3 55.0   80
4 56.0   86
5 56.0   83
6 57.0   88

Plot them on Hudson Bay map, with numbers as labels:

plot(newmap, xlim = c(-90, -80), ylim = c(50, 65))
points(-hb$long, hb$lat)
text(-hb$long, hb$lat, 1:36, pos = 4)

Plotting stuff at points

Normally want to label the points (eg. by colour) according to values of some variable.
Take 2011 breakup dates. Omit 1st column (years), turn it from a row into a column (t()). Find five-number summary.

b = t(breakup[39, -1])
quantile(b)

   0%   25%   50%   75%  100% 
137.0 170.2 179.0 193.0 200.0

Start from 135 and go up in 15s:

bc = cut(b, c(135, 150, 165, 180, 195, 210))
bc

 [1] (165,180] (150,165] (165,180] (165,180] (180,195] (165,180] (180,195]
 [8] (180,195] (165,180] (150,165] (195,210] (180,195] (180,195] (165,180]
[15] (150,165] (165,180] (195,210] (180,195] (165,180] (150,165] (180,195]
[22] (195,210] (195,210] (180,195] (135,150] (165,180] (195,210] (195,210]
[29] (195,210] (135,150] (165,180] (165,180] (150,165] (135,150] (135,150]
[36] (180,195]
Levels: (135,150] (150,165] (165,180] (180,195] (195,210]

bc.col = as.numeric(bc)
bc.col

 [1] 3 2 3 3 4 3 4 4 3 2 5 4 4 3 2 3 5 4 3 2 4 5 5 4 1 3 5 5 5 1 3 3 2 1 1
[36] 4

Choosing colours

Want ordering of colours that makes sense. One way:

library(RColorBrewer)
my.palette = brewer.pal(5, "YlOrRd")

5 colours on a yellow-orange-red scale:

display.brewer.pal(5, "YlOrRd")

Use this as colour (turned into numbers 1, 2, 3, 4, 5).

Map with points labelled by colour for 2011 breakup date

pch=19 is "fillable point", col fills it with colour. Use colours in my.palette:

plot(newmap, xlim = c(-90, -80), ylim = c(50, 65))
points(-hb$long, hb$lat, pch = 19, col = my.palette[bc])
text(-hb$long, hb$lat, 1:36, pos = 4)

Something more sophisticated

Do same thing, but for mean 1971-1990, 1991-2011. Draw both on map.
Strategy:
find "early" and "late" mean for each location
turn them into colours on common scale
plot both colours at each location (move one of them over so that you can see both)

early.mean = apply(breakup[1:20, -1], 2, mean)
late.mean = apply(breakup[21:39, -1], 2, mean)
range(early.mean)

[1] 174.7 210.2

range(late.mean)

[1] 152.2 202.9

150 to 220 should catch everything.

Convert breakup dates to categories to numbers

intervals = seq(150, 220, 10)
intervals

[1] 150 160 170 180 190 200 210 220

early.cut = cut(early.mean, intervals)
early.cut

 [1] (180,190] (170,180] (180,190] (200,210] (200,210] (200,210] (210,220]
 [8] (200,210] (190,200] (180,190] (200,210] (200,210] (200,210] (170,180]
[15] (190,200] (200,210] (200,210] (190,200] (170,180] (200,210] (200,210]
[22] (200,210] (190,200] (180,190] (180,190] (190,200] (190,200] (190,200]
[29] (180,190] (170,180] (180,190] (180,190] (180,190] (170,180] (170,180]
[36] (180,190]
7 Levels: (150,160] (160,170] (170,180] (180,190] (190,200] ... (210,220]

late.cut = cut(late.mean, intervals)
early.val = as.numeric(early.cut)
early.val

 [1] 4 3 4 6 6 6 7 6 5 4 6 6 6 3 5 6 6 5 3 6 6 6 5 4 4 5 5 5 4 3 4 4 4 3 3
[36] 4

late.val = as.numeric(late.cut)
late.val

 [1] 3 2 3 5 5 5 6 5 5 2 6 6 5 3 5 5 5 5 3 3 5 5 5 3 2 4 5 5 3 1 3 4 2 1 1
[36] 3

Getting our colours

7 categories
So need 7 yellow-orange-red colours:

my.palette = brewer.pal(7, "YlOrRd")
display.brewer.pal(7, "YlOrRd")

Making the plot

2 points lines: one for 1991-2010, another (shifted left a little) for the 1971-1990 values:

plot(newmap, xlim = c(-90, -80), ylim = c(50, 65))
points(-hb$long, hb$lat, pch = 19, col = my.palette[late.val])
points(-hb$long - 0.5, hb$lat, pch = 19, col = my.palette[early.val])
text(-hb$long, hb$lat, 1:36, pos = 4)

Interpreting the plot

Yellow=early breakup, orange=medium, red=late, dark red=very late.
Overall colour of two spots for each location tells about overall earliness/lateness.
pattern similar at both time points.
Top left earliest
then right
then bottom left and centre.
Comparison of colours for same location tells about time trends.
Usually right-hand spot lighter than left one: trend towards earlier breakup over time.