Reading: R Cookbook, Chapter 6
Agenda for today:
Homework 1 solutions
Avoiding iteration in base R
Split/apply/combine
We've already seen a lot of examples of functions in R that are written so as to avoid iteration.
Built-in functions for making sequences or repeating elements, for example, we can say:
x = seq(1, 11, by = 2)instead of
x = numeric(6)
for(i in 1:length(x)) {
x[i] = 1 + 2 * (i - 1)
}All the vectorized functions, so we can use
log_x = log(x)instead of
log_x = numeric(length(x))
for(i in 1:length(x)) {
log_x[i] = log(x[i])
}Plus a whole bunch we haven't seen yet, like rowSums, colMeans, sweep, scale
However, not everything is vectorized, and you will often want to go beyond the built-in functions, and for that you use the *apply family of functions.
lapply, and sapply will apply a function to a one-dimensional structure.
The difference between them is how the output is formatted.
They have the same syntax: lapply(x, fun, ...), sapply(x, fun, ...)
x is a vector or list that the function should be applied to
fun is the function that should be applied to each element of the vector or the list
You can specify other arguments in the ... portion. fun will be called with first argument equal to each element of x, and subsequent arguments anything you specify in the ... portion.
Return values:
lapply returns a list with length equal to the length of x containing the result of the function evaluation on the corresponding element of x.
sapply tries to simplify the output to a vector or matrix. If each element of the output is of length 1 it returns a vector. If each element has the same length of two or more, it returns a matrix. Otherwise it returns a list like lapply.
Examples from the homework:
library(readr)
book_dir = "~/GitHub/stat-comp-spring-19/homeworks/hw1/books"
book_file_names = list.files(book_dir, full.names = TRUE)
books = lapply(book_file_names, read_file)We need to do this because read_file doesn't vectorize. Check what happens if you do read_file(book_file_names) (you'll have to change book_dir to wherever you saved the books to).
Example using extra function arguments:
get_num_characters = function(file_name, honorific_regex) {
text = read_file(file_name)
name_matches = regmatches(text, gregexpr(honorific_regex, text))[[1]]
return(length(unique(name_matches)))
}
first_five_books = list.files(book_dir, full.names = TRUE)[1:5]
lapply(first_five_books, get_num_characters, "(((Mr|Ms|Mrs|Dr)\\.?)|Miss) ([A-Z]([a-z]+|\\.)[ ]+)+[A-Z][a-z]+")## [[1]]
## [1] 0
##
## [[2]]
## [1] 7
##
## [[3]]
## [1] 8
##
## [[4]]
## [1] 40
##
## [[5]]
## [1] 2sapply(first_five_books, get_num_characters, "(((Mr|Ms|Mrs|Dr)\\.?)|Miss) ([A-Z]([a-z]+|\\.)[ ]+)+[A-Z][a-z]+")## /Users/julia/GitHub/stat-comp-spring-19/homeworks/hw1/books/1289-0.txt
## 0
## /Users/julia/GitHub/stat-comp-spring-19/homeworks/hw1/books/1400-0.txt
## 7
## /Users/julia/GitHub/stat-comp-spring-19/homeworks/hw1/books/564-0.txt
## 8
## /Users/julia/GitHub/stat-comp-spring-19/homeworks/hw1/books/580-0.txt
## 40
## /Users/julia/GitHub/stat-comp-spring-19/homeworks/hw1/books/653-0.txt
## 2Syntax: apply(mat, dim, fun, ...)
mat: The matrix or data frame you want to apply the function to.
dim: Either 1 or 2: 1 means apply the function to each row in mat, 2 means apply the function to each column in mat. (1 and 2 are related to the way the dimensions are laid out: the first dimension in the matrix is the number of rows, the second is the number of columns.)
fun: The function you want to apply to each row or to each column. It should take a vector as its first argument.
...: Extra arguments to be passed to fun.
As with lapply and sapply, fun will be called with each row or each column of mat as its first argument, along with any extra arguments specified after the function.
Some simple examples:
X = matrix(sample(-5:5, size = 20, replace = TRUE), nrow = 4, ncol = 5)
X## [,1] [,2] [,3] [,4] [,5]
## [1,] 4 4 2 -4 0
## [2,] -3 -3 1 2 2
## [3,] -1 4 -5 -1 5
## [4,] 1 5 -3 3 -1apply(X, 2, min)## [1] -3 -3 -5 -4 -1log_X = log(X)## Warning in log(X): NaNs producedlog_X## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.386294 1.386294 0.6931472 NaN -Inf
## [2,] NaN NaN 0.0000000 0.6931472 0.6931472
## [3,] NaN 1.386294 NaN NaN 1.6094379
## [4,] 0.000000 1.609438 NaN 1.0986123 NaNapply(log_X, 1, min, na.rm = TRUE)## [1] -Inf 0.000000 1.386294 0.000000Be careful when using apply on rows of a data frame: apply assumes the argument you pass to the function is a vector, and will coerce the rows of the data frame to vectors to make it so.
steak_combinations = data.frame(
temp = c(116, 120, 135, 105),
type = c("rare", "med_rare", "med_rare", "rare"))
steak_directions_for_apply = function(temp_and_time) {
temp = temp_and_time[1]
steak_type = temp_and_time[2]
if(steak_type == "rare" & temp > 115) {
return("take your steak off!")
} else if(steak_type == "med_rare" & temp > 125) {
return("take your steak off!")
}
"you can keep cooking"
}
## this works, but it really shouldn't
apply(steak_combinations, 1, steak_directions_for_apply)## [1] "take your steak off!" "you can keep cooking" "take your steak off!"
## [4] "you can keep cooking"If you want to apply a function that uses columns of a data frame that are of different types, better/safer to use mapply:
mapply(f, vec1, vec2, ..., vecN)
f is the function you want to apply
Output is a list, ith element will be f(vec1[i], vec2[i], ..., vecN[i])
steak_directions = function(steak_type, temp) {
if(steak_type == "rare" & temp > 115) {
return("take your steak off!")
} else if(steak_type == "med_rare" & temp > 125) {
return("take your steak off!")
}
"you can keep cooking"
}
mapply(FUN = steak_directions, steak_combinations$type, steak_combinations$temp)## [1] "take your steak off!" "you can keep cooking" "take your steak off!"
## [4] "you can keep cooking"The *apply functions we've seen so far have applied a function element-by-element to a vector or matrix.
What if we need to apply a function not to individual elements/rows/columns, but to subsets of different sizes?
Turns out this is common enough that it has a name and a bunch of built-in functions.
The programming pattern is pretty clear from the name, what we want to do is:
split the data into groups
apply a function to each of the groups
combine the results of applying the function to each group into an output dataset
Imagine how you would do this if you only had for-loops.
You would have to first
Compute what size the output dataset should be
Create a data structure to hold the output.
Then, for each level of the grouping factor:
Compute how many data elements correspond to that grouping level
Create a data structure to hold the data subset
Loop through the data to extract the correct subset
Apply your function to that subset and put the output in the data structure your created at the beginning.
Compare with split/apply/combine abstraction:
Specify data
Specify grouping variable
Specify function to apply
The abstraction is useful both for writing and reading the resulting code: you think about what you want to do instead of how to do it.
Simple case: we have a vector, and a grouping factor, we want to apply a function to groups of elements of the vector defined by the grouping factor
Solution is tapply
tapply(x, f, fun)
Syntax:
x: A vector with the data.
f: A factor variable describing how the data should be split up.
fun: The function you want to apply to each subset of the data.
Example: data on police stops in New York City in 1998--1999, information on the number of stops as categorized by certain variables.
In particular, we have counts of police stops for all combinations of
precinct: 75 total
eth: Ethnicity of the person stopped, three possibilities (1 = black, 2 = Hispanic, 3 = white), and
crime: The type of crime, four possibilities (1 = violent, 2 = weapons, 3 = property, and 4 = drug)
frisk =
read.table("http://www.stat.columbia.edu/~gelman/arm/examples/police/frisk_with_noise.dat",
skip = 6, header = TRUE)
head(frisk)## stops pop past.arrests precinct eth crime
## 1 75 1720 191 1 1 1
## 2 36 1720 57 1 1 2
## 3 74 1720 599 1 1 3
## 4 17 1720 133 1 1 4
## 5 37 1368 62 1 2 1
## 6 39 1368 27 1 2 2Suppose we want to get the total number of stops in each precinct.
We can use tapply with the sum function to do so:
tapply(frisk$stops, frisk$precinct, sum)## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
## 385 347 1603 1411 1600 1297 649 1572 285 1156 838 884 2215 1165 2373
## 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
## 965 946 2359 1311 1566 1613 4030 2334 2622 3685 1541 2005 902 2888 2157
## 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
## 1805 1041 2966 2810 724 1189 997 1808 2464 859 1979 2679 1645 1669 2953
## 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
## 3443 1001 1574 1041 3054 1208 2116 1216 995 1140 572 856 2671 1147 2941
## 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
## 1625 1635 1248 1351 2290 2941 2430 1109 1267 2476 2922 3560 3650 1327 322Just as a check, make sure that the numbers match up:
precinct_1 = subset(frisk, precinct == 1)
sum(precinct_1$stops)## [1] 385And remember how much better this is than iterating using a for loop would have been:
## find the number of precincts and make a vector of the right length to store the output
precincts = unique(frisk$precinct)
stops_per_precinct = numeric(length(precincts))
## iterate over precinct, get the subset for each one, compute the number of stops, and store
for(i in 1:length(precincts)) {
precinct_subset = subset(frisk, precinct == precincts[i])
precinct_stops = sum(precinct_subset$stops)
stops_per_precinct[i] = precinct_stops
}
stops_per_precinct## [1] 385 347 1603 1411 1600 1297 649 1572 285 1156 838 884 2215 1165
## [15] 2373 965 946 2359 1311 1566 1613 4030 2334 2622 3685 1541 2005 902
## [29] 2888 2157 1805 1041 2966 2810 724 1189 997 1808 2464 859 1979 2679
## [43] 1645 1669 2953 3443 1001 1574 1041 3054 1208 2116 1216 995 1140 572
## [57] 856 2671 1147 2941 1625 1635 1248 1351 2290 2941 2430 1109 1267 2476
## [71] 2922 3560 3650 1327 322Splitting a data frame into groups of rows and applying a function to the groups:
by(df, fact, fun)
df: A data frame containing the data you want to split up and apply a function to.
fact: A factor variable describing the groups you want to split df into.
fun: The function you want to apply to each group. It should take as its first argument a data frame.
Example:
data(diamonds)
get_diamond_coefficients = function(data_subset) {
diamond_lm = lm(log(price) ~ carat, data = data_subset)
diamond_coefficients = coef(diamond_lm)
return(diamond_coefficients)
}
## by does the split and apply steps, and an ugly version of a combine step
out = by(diamonds, diamonds$color, get_diamond_coefficients)The output is as a list with some extra attributes (the class is by, as you can check).
The form is not ideal on its own, usually you'll want to simplify it. You have a couple of options, including the examples below:
simplify2array(out)## D E F G H I J
## (Intercept) 6.048811 6.034513 6.088442 6.109554 6.180284 6.175315 6.254074
## carat 2.383864 2.348335 2.272790 2.178489 1.906300 1.799199 1.627947## changing the output from a list to a matrix
do.call(rbind, out)## (Intercept) carat
## D 6.048811 2.383864
## E 6.034513 2.348335
## F 6.088442 2.272790
## G 6.109554 2.178489
## H 6.180284 1.906300
## I 6.175315 1.799199
## J 6.254074 1.627947## it's kind of bad that the splitting variable is only available as row names
## we can fix that but it's a bit of a pain
diamond_coef_matrix = do.call(rbind, out)
diamond_coef_df = data.frame(diamond_coef_matrix,
color = attributes(out)$dimnames$`diamonds$color`)
diamond_coef_df## X.Intercept. carat color
## D 6.048811 2.383864 D
## E 6.034513 2.348335 E
## F 6.088442 2.272790 F
## G 6.109554 2.178489 G
## H 6.180284 1.906300 H
## I 6.175315 1.799199 I
## J 6.254074 1.627947 Jby can also be used with multiple grouping factors:
by(diamonds, diamonds[, c("color", "cut")], get_diamond_coefficients)## color: D
## cut: Fair
## (Intercept) carat
## 6.723804 1.514938
## --------------------------------------------------------
## color: E
## cut: Fair
## (Intercept) carat
## 6.246003 1.955816
## --------------------------------------------------------
## color: F
## cut: Fair
## (Intercept) carat
## 6.466073 1.637304
## --------------------------------------------------------
## color: G
## cut: Fair
## (Intercept) carat
## 6.659407 1.358722
## --------------------------------------------------------
## color: H
## cut: Fair
## (Intercept) carat
## 6.921231 1.122101
## --------------------------------------------------------
## color: I
## cut: Fair
## (Intercept) carat
## 6.759325 1.202279
## --------------------------------------------------------
## color: J
## cut: Fair
## (Intercept) carat
## 7.0504140 0.8761416
## --------------------------------------------------------
## color: D
## cut: Good
## (Intercept) carat
## 5.925962 2.415309
## --------------------------------------------------------
## color: E
## cut: Good
## (Intercept) carat
## 5.976449 2.322896
## --------------------------------------------------------
## color: F
## cut: Good
## (Intercept) carat
## 6.016284 2.273226
## --------------------------------------------------------
## color: G
## cut: Good
## (Intercept) carat
## 6.180670 2.037179
## --------------------------------------------------------
## color: H
## cut: Good
## (Intercept) carat
## 6.106822 1.952841
## --------------------------------------------------------
## color: I
## cut: Good
## (Intercept) carat
## 6.182043 1.762664
## --------------------------------------------------------
## color: J
## cut: Good
## (Intercept) carat
## 6.133176 1.728532
## --------------------------------------------------------
## color: D
## cut: Very Good
## (Intercept) carat
## 5.953898 2.487391
## --------------------------------------------------------
## color: E
## cut: Very Good
## (Intercept) carat
## 5.909810 2.480529
## --------------------------------------------------------
## color: F
## cut: Very Good
## (Intercept) carat
## 5.947350 2.446917
## --------------------------------------------------------
## color: G
## cut: Very Good
## (Intercept) carat
## 5.965330 2.343517
## --------------------------------------------------------
## color: H
## cut: Very Good
## (Intercept) carat
## 6.110954 1.993643
## --------------------------------------------------------
## color: I
## cut: Very Good
## (Intercept) carat
## 6.207955 1.807934
## --------------------------------------------------------
## color: J
## cut: Very Good
## (Intercept) carat
## 6.211086 1.696800
## --------------------------------------------------------
## color: D
## cut: Premium
## (Intercept) carat
## 6.109571 2.258907
## --------------------------------------------------------
## color: E
## cut: Premium
## (Intercept) carat
## 6.133275 2.185096
## --------------------------------------------------------
## color: F
## cut: Premium
## (Intercept) carat
## 6.189666 2.124037
## --------------------------------------------------------
## color: G
## cut: Premium
## (Intercept) carat
## 6.181514 2.052932
## --------------------------------------------------------
## color: H
## cut: Premium
## (Intercept) carat
## 6.270586 1.809367
## --------------------------------------------------------
## color: I
## cut: Premium
## (Intercept) carat
## 6.229239 1.723257
## --------------------------------------------------------
## color: J
## cut: Premium
## (Intercept) carat
## 6.305114 1.578440
## --------------------------------------------------------
## color: D
## cut: Ideal
## (Intercept) carat
## 5.970971 2.626088
## --------------------------------------------------------
## color: E
## cut: Ideal
## (Intercept) carat
## 5.989833 2.526368
## --------------------------------------------------------
## color: F
## cut: Ideal
## (Intercept) carat
## 6.052536 2.408753
## --------------------------------------------------------
## color: G
## cut: Ideal
## (Intercept) carat
## 6.032829 2.367710
## --------------------------------------------------------
## color: H
## cut: Ideal
## (Intercept) carat
## 6.075676 2.071994
## --------------------------------------------------------
## color: I
## cut: Ideal
## (Intercept) carat
## 6.071246 1.932255
## --------------------------------------------------------
## color: J
## cut: Ideal
## (Intercept) carat
## 6.104631 1.784315But how do you get at which coefficients correspond to which factor groups?
split does just the splitting part of split/apply/combine
Syntax: split(x, f)
x: Your data, either a vector or a data frame, that you want to divide into groups.
f: A factor variable used to split x. Length is equal to the length of x if it is a vector or equal to the number of rows of x if it is a data frame.
Output:
A list of length equal to the number of levels in f
Each element of the list is either a vector or a data frame, depending on what x was, containing the elements of x corresponding to one of the levels of f.
We can re-do the stop-and-frisk example that we used tapply for with split followed by lapply or sapply
frisk_split_by_precinct = split(frisk$stops, frisk$precinct)
## check that this has length 75, the number of precincts
length(frisk_split_by_precinct)## [1] 75sapply(frisk_split_by_precinct, sum)## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
## 385 347 1603 1411 1600 1297 649 1572 285 1156 838 884 2215 1165 2373
## 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
## 965 946 2359 1311 1566 1613 4030 2334 2622 3685 1541 2005 902 2888 2157
## 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
## 1805 1041 2966 2810 724 1189 997 1808 2464 859 1979 2679 1645 1669 2953
## 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
## 3443 1001 1574 1041 3054 1208 2116 1216 995 1140 572 856 2671 1147 2941
## 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
## 1625 1635 1248 1351 2290 2941 2430 1109 1267 2476 2922 3560 3650 1327 322If you have a function you only want to use once, you don't have to assign it to anything.
When you use a function this way it is called an anonymous function because it doesn't have a name in your code.
These are often used in the context of the apply family.
diamonds_split = split(diamonds, diamonds$color)
## here the function for computing the coefficients is an anonymous function
diamond_coef_list = lapply(diamonds_split, FUN = function(data_subset) {
diamond_lm = lm(log(price) ~ carat, data = data_subset)
diamond_coefficients = coef(diamond_lm)
return(diamond_coefficients)
})
## changing the output from a list to a matrix
do.call(rbind, diamond_coef_list)## (Intercept) carat
## D 6.048811 2.383864
## E 6.034513 2.348335
## F 6.088442 2.272790
## G 6.109554 2.178489
## H 6.180284 1.906300
## I 6.175315 1.799199
## J 6.254074 1.627947In general, base R apply-family functions don't produce nice output from the combine step.
We'll see implementations of split/apply/combine next time that do a nicer job.