Vectorization + Functions

Week 5 · Thursday, Oct 22 · Lecture

Dr. Marcela Alfaro Córdoba

Today’s plan

Time What
1:30 DS 4 discussion check
1:45 Vectorization: one operation, every element
2:15 Functions: name a piece of work
2:40 Functions that fail loudly, and functions on data frames
2:55 Tuesday: how the day will run

Learning goals: replace a loop with one vectorized line, predict what recycling does, and write a small function that takes arguments, returns a value, and stops on bad input.

Laptops closed for the first 15 minutes. Still no AI.

DS 4 discussion check

How this works

  1. Your DS 4 PDF (Lesson 8) in front of you: printed, or open on your laptop. Nothing else open
  2. I ask for volunteers. Raise your hand when you have a plot you’re ready to talk through
  3. I point at a line of your own code and ask what it does, or what changes if you move it
  4. “I’m not sure, but I think…” followed by reasoning is a good answer

Note

From now on, printing is optional. A printout is still the easiest way to scribble on your code, so keep doing it if it helps you. The PDF on Canvas is what’s graded.

The kind of question you’ll get

1. In Exercise 3, which line made the lines draw on top of the points? What happens if you swap the two layers?

The layer added later is drawn on top. Swap them and the points cover the lines where they overlap.

2. In Exercise 5, how many lm lines did you get? Point to the line of code that decides that number.

Five, one per continent, because colour = continent is in ggplot() and geom_smooth() inherits it. Move it into geom_point() and you get one line.

3. You wrote colour = "red" inside aes(). What does the legend say, and what colour are the points?

A legend with one entry, red, and salmon points: inside aes(), "red" is treated as data.

Vectorization

Today’s script

Create week05/vectorization-functions.R:

# Week 5 · Vectorization + Functions
source("scripts/load-gapminder.R")
source("scripts/check-gapminder.R")
source("scripts/clean-survey.R")

Three source() lines and you have both datasets, loaded and checked. That’s what Weeks 3 and 4 built.

One operation, every element

x <- c(2, 4, 6)
x * 10          #> [1] 20 40 60
x + c(1, 2, 3)  #> [1]  3  6  9
sqrt(x)         #> [1] 1.414214 2.000000 2.449490

You’ve been doing this since Week 2:

gapminder$pop_millions <- gapminder$pop / 1e6   # in load-gapminder.R

1704 divisions, one line, no loop.

Two columns: row 1 with row 1

gapminder$gdp_billions <- gapminder$pop * gapminder$gdpPercap / 1e9
head(gapminder[, c("country", "year", "pop", "gdpPercap", "gdp_billions")], 3)
#>       country year      pop gdpPercap gdp_billions
#> 1 Afghanistan 1952  8425333  779.4453     6.567086
#> 2 Afghanistan 1957  9240934  820.8530     7.585449
#> 3 Afghanistan 1962 10267083  853.1007     8.758856
  • pop[1] * gdpPercap[1], then pop[2] * gdpPercap[2], and so on, for every row
  • Both columns have 1704 values, so they pair up exactly

The loop you already know, and the line that replaces it

Week 3: a loop

squares <- numeric(5)
for (i in 1:5) {
  squares[i] <- i ^ 2
}
squares
#> [1]  1  4  9 16 25

Today: one line

squares <- (1:5) ^ 2
squares
#> [1]  1  4  9 16 25

Tip

Vectorized code is shorter, easier to read, and much faster. For a million values, the loop above is many times slower. Loops still have their place: when each step depends on the previous one, or when you repeat a whole block (a plot, a check) per group.

Recycling: when lengths differ

R reuses the shorter vector from the start until it covers the longer one:

1:6 * c(1, -1)          #> [1]  1 -2  3 -4  5 -6
1:6 + c(100, 200, 300)  #> [1] 101 202 303 104 205 306

If the longer length isn’t a multiple of the shorter one, you get an answer and a warning:

c(2, 4, 6) * c(1, 10)
#> [1]  2 40  6
#> Warning message:
#> longer object length is not a multiple of shorter object length

Warning

A length-1 vector (x * 10) is recycling too, and that’s the one you want. Any other mismatch is almost always a bug. Read that warning.

Your turn: predict, then run

⏱ 5 min

On paper first. Write the result of each, then run them:

m <- matrix(1:6, nrow = 2)
m

m * c(1, -1)

m + c(10, 20, 30)

Hint: a matrix is a vector stored column by column.

m * c(1, -1)            m + c(10, 20, 30)
#>      [,1] [,2] [,3]   #>      [,1] [,2] [,3]
#> [1,]    1    3    5   #> [1,]   11   33   25
#> [2,]   -2   -4   -6   #> [2,]   22   14   36

The first one is clean: row 2 is always multiplied by −1. The second one runs 10, 20, 30 down the columns, not across the rows. No warning, because 6 is a multiple of 3. If you expected 11 23 35 on the first row, that’s the trap.

Comparisons are vectorized too

gm07 <- gapminder[gapminder$year == 2007, ]

gm07$lifeExp > 70          # 142 TRUE/FALSE values, one per country
sum(gm07$lifeExp > 70)     #> [1] 83        TRUE counts as 1
mean(gm07$lifeExp > 70)    #> [1] 0.584507  the proportion

ifelse(): a vectorized if, one answer per element:

gm07$le_band <- ifelse(gm07$lifeExp > 70, "long", "short")
table(gm07$le_band)
#>  long short
#>    83    59

sum(condition) counts, mean(condition) gives a proportion. You’ll use these two in every project.

NA spreads, again

x <- c(3, NA, 5)
x > 4                    #> [1] FALSE    NA  TRUE
sum(x > 4)               #> [1] NA
sum(x > 4, na.rm = TRUE) #> [1] 1

Note

Same rule as Week 3: one NA in, NA out. With na.rm = TRUE you drop the NAs, so first know how many you’re dropping: sum(is.na(x)).

& and && are not the same

Works on Use it in
&, \| whole vectors, element by element subsetting, ifelse(), new columns
&&, \|\| one TRUE/FALSE each if () conditions
gapminder[gapminder$continent == "Asia" & gapminder$year == 2007, ]   # vector: &

if (nrow(survey) == 45 && !any(duplicated(survey))) message("ok")    # one value: &&

any() and all() turn a vector of TRUE/FALSE into one value, so they can go inside if (). That’s how your check scripts work.

Your turn: the survey, vectorized

⏱ 8 min

In your script, one line each, no loops:

  1. A new column commute_h: the commute in hours
  2. How many people commute more than 20 minutes? What proportion?
  3. A new column long_commute: "yes" or "no", for more than 20 minutes
  4. Check: does table(survey$long_commute, useNA = "ifany") agree with step 2?
survey$commute_h <- survey$commute_min / 60
sum(survey$commute_min > 20, na.rm = TRUE)     #> [1] 17
mean(survey$commute_min > 20, na.rm = TRUE)    #> [1] 0.3863636
survey$long_commute <- ifelse(survey$commute_min > 20, "yes", "no")

17 yes, 27 no, 1 NA (id 35, the missing commute). ifelse() keeps the NA: it doesn’t guess. The proportion is out of 44, not 45.

Functions

Why write a function?

You just wrote this:

mean(survey$commute_min > 20, na.rm = TRUE)
mean(survey$commute_min > 30, na.rm = TRUE)
mean(survey$visits > 3, na.rm = TRUE)
mean(gapminder$lifeExp > 70)

Copied the same line three times? Time for a function.

  • Name the idea: “proportion above a cutoff”
  • Fix a mistake in one place
  • Test it once, trust it everywhere

Anatomy of a function

prop_above <- function(x, cutoff) {
  above <- x > cutoff
  mean(above, na.rm = TRUE)
}
Part Here
name prop_above
the keyword function
arguments x, cutoff
body inside { }
returned value the last line, or whatever you give to return()

Calling it

prop_above(survey$commute_min, 20)    #> [1] 0.3863636
prop_above(survey$commute_min, 30)    #> [1] 0.1590909
prop_above(survey$visits, 3)          #> [1] 0.5555556
prop_above(gapminder$lifeExp, 70)     #> [1] 0.2893192

By name, in any order:

prop_above(cutoff = 20, x = survey$commute_min)   # same answer

And it’s vectorized for free: inside, x > cutoff works on the whole column, because > does.

What happens inside stays inside

prop_above(survey$commute_min, 20)
above
#> Error: object 'above' not found
  • x, cutoff and above exist only while the function runs
  • When it finishes, only the returned value comes back
  • That’s a feature: your function can’t accidentally overwrite a variable in your script

Warning

The other direction does leak. A function can read variables from your script. Don’t let it:

bad_mean <- function(cont) mean(gapminder$lifeExp[gapminder$continent == cont])

It only works on gapminder. Pass the data in as an argument, and it works on any data frame.

Small functions, combined

min_to_h <- function(min) {
  min / 60
}

h_to_days <- function(h) {
  h / 24
}

min_to_days <- function(min) {
  h_to_days(min_to_h(min))
}

min_to_days(1440)    #> [1] 1

Each piece is easy to check on its own. When min_to_days() gives a weird answer, you test the two small ones first. Same idea as Week 3’s check script: small steps you can trust.

Your turn: write two functions

⏱ 7 min

  1. count_na(x): how many NAs are in x. Test it on survey$age (should be 3) and survey$site (1)
  2. pct(x): turns a proportion into a percentage, rounded to 1 decimal: pct(0.3863636) gives 38.6
  3. Use both: what percentage of the survey commutes more than 20 minutes? One line, using prop_above() and pct()
count_na <- function(x) {
  sum(is.na(x))
}

pct <- function(p) {
  round(100 * p, 1)
}

pct(prop_above(survey$commute_min, 20))    #> [1] 38.6

count_na() on every column at once? sapply(survey, count_na). A function you wrote, applied to each column: we’ll see more of this later.

Fail loudly

Wrong input, wrong answer, no error

prop_above(survey$site, 20)
#> [1] 1

site is text. R compared "east" > 20 as text (alphabetically), and every site name sorts after "20".

“100% of sites are above 20” is nonsense, and R said nothing.

This is the worst kind of bug: a number that looks plausible. A function should stop when it gets input it can’t handle.

stopifnot() and stop()

prop_above <- function(x, cutoff) {
  stopifnot(is.numeric(x), is.numeric(cutoff))
  if (all(is.na(x))) {
    stop("x has no non-missing values")
  }
  mean(x > cutoff, na.rm = TRUE)
}

prop_above(survey$site, 20)
#> Error in prop_above(survey$site, 20) : is.numeric(x) is not TRUE
Use it for
stopifnot(cond1, cond2) quick checks; R writes the message
if (...) stop("message") your own message, for a reader who isn’t you

Same idea as your check scripts from Week 3, now built into the tool.

Functions on data frames, with a default

continent_mean <- function(dat, cont, col = "lifeExp") {
  rows <- dat[dat$continent == cont, ]
  mean(rows[[col]])
}

continent_mean(gapminder, "Africa")                      #> [1] 48.86533
continent_mean(gapminder, "Europe", col = "gdpPercap")   #> [1] 14469.48
  • col = "lifeExp" is a default: leave it out and that’s what you get
  • dat is an argument, so it works on gapminder, gm07, or any gapminder-style data frame
  • rows[[col]] takes a column whose name is stored in a variable

Your turn: make it fail loudly

⏱ 6 min

continent_mean(gapminder, "africa")
#> [1] NaN
  1. Why NaN? (Hint: how many rows does rows have?)
  2. Add a check: if cont is not one of the continents in dat, stop with a message that says which continent wasn’t found
  3. Test it with "africa", then with "Africa"
continent_mean <- function(dat, cont, col = "lifeExp") {
  if (!(cont %in% dat$continent)) {
    stop("continent not found: ", cont)
  }
  rows <- dat[dat$continent == cont, ]
  mean(rows[[col]])
}

Lowercase a, zero matching rows, mean() of nothing is NaN. A silent wrong answer turns into a clear error.

Where functions live in a project

stat266a/
├── data/
├── scripts/
│   ├── load-gapminder.R
│   ├── check-gapminder.R
│   ├── clean-survey.R
│   └── functions.R        ← new
└── week05/
# scripts/functions.R: helper functions for this project
prop_above <- function(x, cutoff) { ... }
count_na <- function(x) { ... }

Then, at the top of any analysis script: source("scripts/functions.R").

Tip

For your project: if you’ve cleaned two columns the same way, that’s your first function.

Key points

  • Arithmetic, comparisons and most functions work on whole vectors, element by element
  • Two columns pair up row by row; a shorter vector is recycled, and a length mismatch warns: read it
  • Matrices recycle down the columns
  • sum(condition) counts, mean(condition) gives a proportion; ifelse() is a vectorized if
  • & for vectors, && inside if ()
  • name <- function(args) { body }: the last line is the result
  • Variables inside a function stay inside. Pass data in as an argument, never use a global
  • Small functions, combined. Fail early with stopifnot() and stop()

Tuesday, Oct 27

Tuesday is a full day

9:50-10:55 am · DS 5

Lesson 9: Vectorization and Lesson 10: Functions, alone, no AI. Two PDFs.

1:30-3:05 pm · Regression + AI audit #1

Starts with the DS 5 check: you swap laptops with a classmate and check each other’s DS. Then linear regression, and your first AI audit.

Important

Bring your laptop, charged, to the afternoon. The DS check happens on your classmate’s screen, not on paper. No printing needed.

Before you go

Before Tuesday: read the Lesson 9 and Lesson 10 pages. The exercises will feel familiar after today.

Next: Tuesday, Oct 27 · DS 5 in the morning, regression in the afternoon.