Worksheet 3

Questions are below. My solutions will be available after the tutorials are all finished. The whole point of these worksheets is for you to use your lecture notes to figure out what to do. In tutorial, the TAs are available to guide you if you get stuck. Once you have figured out how to do this worksheet, you will be prepared to tackle the assignment that depends on it.

If you are not able to finish in an hour, I encourage you to continue later with what you were unable to finish in tutorial.

Questions are respectively on choosing rows and columns, recoding values, and joins. There is a lot here, so you should plan to work at a good pace. You don’t have to do these questions in order. If you want practice on a certain thing, start with that problem.

Packages you’ll need:

library(tidyverse)

Student stats

A survey of students at the University of California Davis asked respondents:

• the gender they identify as (called Sex in the data set)
• hours of TV they watched in the last week
• hours they spent doing non-academic things on a computer in the last week (playing games, on social media etc. This survey was taken some years ago, and so pre-dates smartphones.)
• how many hours of sleep they got last night
• where they prefer to sit in their classes (front, middle, back of classroom)
• how many alcoholic drinks they consumed in the last week (defined as US standard drinks)
• their height (inches)
• their mother’s height (inches)
• their father’s height (inches)
• the number of hours of exercise they got in the last week
• their GPA
• their area of study, classified as Liberal Arts or NonLib (“not arts” meaning science, math, engineering etc).

Students were allowed to not respond to any of the items, so there are some missing values, labelled NA in the dataframe you are about to read in.

The data are in http://datafiles.ritsokiguess.site/UCDavis1.csv.

1. Read in and display (some of) the data.

Solution

As usual:

my_url <- "http://datafiles.ritsokiguess.site/UCDavis1.csv"
davis <- read_csv(my_url)

Rows: 173 Columns: 12
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (3): Sex, Seat, class
dbl (9): TV, computer, Sleep, alcohol, Height, momheight, dadheight, exercis...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

davis

Give the dataframe a suitable name. I called it davis, because it is a bunch of information about students at University of California-Davis.

\(\blacksquare\)

2. Find the mean and standard deviation of sleep times (for all the students taken together).

Solution

A summarize. Use your name for the dataframe, if it is different from mine:

davis %>% summarize(mean_sleep = mean(Sleep), sd_sleep = sd(Sleep))

Use names for the summaries that say something about what they are: what statistic (mean, sd here) and what variable. Using a name like mean by itself will make your reader wonder “mean of what?”, and using an even shorter name like m will make your reader wonder “what is m again?”, if not now, then in a few minutes. This is true even if the only reader is going to be you: when you come back to this work, future-you is going to thank past-you for being clear.

The one possible exception to this is for something you are going to use once and then never again. I often use d as a name for a “temporary” dataframe, for example one that I create to make a graph out of, and then never use again. But for numerical summaries, it is pretty much always good to be clear about what they are.

\(\blacksquare\)

3. For the students that sit in each part of the classroom (separately), find the number of students and their mean sleep time. (Hint: the column is called Seat with a capital S; the S being uppercase matters.)

Solution

This is a group-by and summarize, grouping by the categorical variable Seat first, and then computing the required summaries. We are finding something else in addition to the counts (the mean sleep time), so we have to use the n() idea here:

# summary(davis)
davis %>% 
  group_by(Seat) %>% 
  summarize(n_student = n(), mean_sleep = mean(Sleep))

There were two students who didn’t say where they liked to sit, so their Seat is missing.

Extra:

There seems to be a tendency for students that sit nearer the front of the class to sleep less. This is a bit hard to read from this table, because the values of Seat come out in alphabetical order, rather than a more easily interpretable order like “front, middle, back”.¹

Having seen this, a graph might give us a hint about whether this is a real trend or more likely just chance:

ggplot(davis, aes(x = Seat, y = Sleep)) + geom_boxplot()

The median sleep times (for the students with non-missing Seat) seem to be exactly equal, and there is a lot of variability, so the trend we saw is almost certainly just chance.

\(\blacksquare\)

4. For each of the students this data set, display how much time they have spent watching TV and on the computer.

Solution

These two things are columns in the dataframe, and choosing columns is select:

davis %>% select(TV, computer)

It is rather annoying that the column computer has a lowercase C, and this matters as far as R is concerned.²

\(\blacksquare\)

5. Several of the columns are heights. Display all of these columns, without naming any columns.

Solution

The first thing is to find out exactly what it is that these height column names have in common. They are Height, momheight, and dadheight. (They are the student’s height and the heights of their parents.) What they have in common is that they end with height:

davis %>% select(ends_with("height"))

Given what we were just saying about uppercase and lowercase, you might be surprised that this worked: after all, the one that is the student’s height starts with an uppercase H, which does not match what I put in the ends_width. It turns out that select-helpers are not case-sensitive, unless you make them so:

davis %>% select(ends_with("height", ignore.case = FALSE))

The option is the rather cumbersome double negative “don’t ignore case”, ie. pay attention to the case. This gets only the two columns whose names literally end in height, not the student’s height whose name is (ends with) Height with a capital H, and this no longer matches because we are paying attention to uppercase/lowercase.

Another way to get all three columns, if you are worried about upper and lower case, is to match on the ends of the names that are the same including case:

davis %>% select(ends_with("eight"))

Fortunately, there are no columns ending in weight here, so this gets only the columns we want.

Or you could say that these three columns are the only ones that contain height, which is also true:

davis %>% select(contains("height"))

or if you really want to grapple with regular expressions, even this:

davis %>% select(matches("height$"))

Inside matches, my regular expression searches for columns that match with the letters height at the end of the column name, so that eg height_in_cm wouldn’t match.³

\(\blacksquare\)

6. Display only the students who sit at the back of the classroom.

Solution

This is choosing rows, specifically the rows for which Seat is Back, so it needs filter:

davis %>% filter(Seat == "Back")

Don’t forget the two equals signs, which is how R tests that something is equal to something else: “give me all the rows for which it is true that Seat is equal to the text Back”.

\(\blacksquare\)

7. Display all the students that either slept 5 hours or less or watched over 30 hours of TV (or both).

Solution

This is either/or, so needs when_any:

davis %>% filter(when_any(
  Sleep <= 5,
  TV > 30))

To be precise, “5 hours or less” is less than or equal, and “over” is strictly greater than.

To convince yourself that you have the right ones, scroll through your results and say why each of the students is there. My first ten students are all there because they slept 5 hours or less (the 8th one also watched over 30 hours of TV). In the next ten, the first student got more than 5 hours of sleep but is there because they (claimed to have) watched 100 hours of TV. Students 24 and 25 got lots of sleep but also watched lots of TV. And so on.

If you don’t have when_any, you’ll need to look in the workarounds for how you say “or”:

davis %>% 
  filter(Sleep <= 5 | TV > 30)

\(\blacksquare\)

8. For the students who are 70 or more inches tall, what are the mean heights of their mother and their father? Hint: the first time you do this, your answers will probably be missing. Why is that? Hint: to get the answers you were expecting, search for the na.rm option. What does that do?

Solution

In this kind of problem, where we want summaries for only some of the observations, the first step is to select the observations we want, and the second step is to summarize the observations we selected:

davis %>% 
  filter(Height >= 70) %>% 
  summarize(mom_mean = mean(momheight), dad_mean = mean(dadheight))

OK, so those were both missing. To find out why, step back to the filter (run it without the summarize):

davis %>% filter(Height >= 70) 

Page down almost to the end to see that the 40th student didn’t give an answer for their mother’s and father’s heights. As far as R is concerned, the mean of a set of numbers including a missing value is itself missing, on the basis that the missing value could be anything.

Here, it seems reasonable to calculate the means of the non-missing parental heights.⁴ There are a couple of ways to do that, which you probably haven’t seen before. The first includes the missing values, but removes them when calculating the mean:⁵

davis %>% 
  filter(Height >= 70) %>% 
  summarize(mom_mean = mean(momheight, na.rm = TRUE), dad_mean = mean(dadheight, na.rm = TRUE))

When searching for an R function or option, it usually helps to search with r first. When I searched for “r na.rm”, I got several hits, any of which would be good. For example, https://www.statology.org/na-rm/ says clearly how to calculate the mean without including any missing values, which is exactly what we want to do.

Aside: I hope you know better than to use an LLM like chatgpt for this kind of job. When you search with an actual search engine, you can see where the results came from and can assess the sources for reliability.⁶ If you are using Google, this means ignoring anything that says “AI” and looking at the actual search results. As you gain experience searching for things in R, you will come to know which websites are accurate and helpful, and which ones read like AI slop.

A second way removes the missing momheights and dadheights first before trying to calculate the means:

davis %>% filter(Height >= 70) %>% 
  drop_na(momheight, dadheight) %>% 
  summarize(mom_mean = mean(momheight), dad_mean = mean(dadheight))

To search for this kind of approach, you can try something like “r omit rows with missing values”. The third answer here gives the Tidyverse answer (the one you want). You’ll see a lot of discussion of the old-fashioned na.omit, but the fourth example here also shows the “tidy” way to do it.

The answer, both ways, is the same. If you run the first two lines of this (and not the third one), you’ll see that there are now only 41 students chosen, and the one with the missing momheight and dadheight is no longer there. So when you run mean, you will not run into any missing values.

Extras:

Extra 1: when you feed drop_na more than one column, it removes any rows that contain missing values in any of the columns. For example, consider this dataframe:⁷

d

What does this do?

d %>% drop_na(y, z)

There is only the first row left, because the second row has a missing value for y and the third has a missing value for z. The values 5 for z and 7 for y, even though they are good data values, are not included because the rows they were in had missing values for something else. The moral of this story is that if you use drop_na on several columns, you stand to lose quite a lot of your data. The default for drop_na is to drop rows with any missing values at all. The original davis dataframe had 173 rows, but if we drop all the rows with any missing values:

davis %>% drop_na()

we are down to 148 rows, having lost almost 15% of our data. Some people prefer to “impute” missing data, which is to replace missing values with an estimate based on the values of other variables. For example, a student with a missing Height whose parents are both tall is probably tall themselves.

Extra 2: on that note, you would expect taller students to have taller parents. How do the values we calculated above compare with the mean heights of all the parents?

davis %>% 
  summarize(mom_mean = mean(momheight, na.rm = TRUE), dad_mean = mean(dadheight, na.rm = TRUE))

An inch or two less.

Extra 3: You’ll remember group-by and summarize from numerical summaries. Is there a way to use that here? Sort of:

davis %>% 
  group_by(Height >= 70) %>% 
  summarize(n = n(), mom_mean = mean(momheight, na.rm = TRUE), dad_mean = mean(dadheight, na.rm = TRUE))

The answers in the TRUE line are the ones we got before. We are used to using group_by with a categorical column, which we don’t have for this problem, but in fact you can use group_by with anything that makes a categorical, including something that evaluates to TRUE or FALSE. The expression Height >= 70 evaluates to true or false or missing (there were two students who didn’t answer how tall they were). The first line of the result gives the mean parental heights for the students who were less than 70 inches tall (most of them), and the last line gives the mean parental heights for the students who didn’t give their own heights.

Note that some of the parental heights were missing, even for students who gave their own heights, so I had to use the na.rm thing again (or I could have used drop_na as before, before calculating the means).

\(\blacksquare\)

Diet Wars 2

In an earlier worksheet, we investigated dropout rates for subjects on three different diets. The data were in http://datafiles.ritsokiguess.site/diet-wars-1.txt. To repeat the data description from earlier:

312 subjects took part in a comparative study of three diets: a low-fat diet (labelled 1 in the dataset), a Mediterranean diet (labelled 2), and a low-carbohydrate diet (labelled 3). Each subject was randomly allocated to one of the diets. One measure of success of a diet program is (of course) how much weight a subject lost, but another is whether they stayed with their diet until the end of the study (1 in the Status column) or whether they dropped out of the study (0 in the Status column).

The various different diets and status values are represented by numbers. In the questions that follow, we will be reorganizing the data so that the categories are represented by something more meaningful.

9. Read in and display (some of) the data. For this, re-use your code from the earlier worksheet.

You’ll recall that these were aligned columns separated by more than one space:

my_url <- "http://datafiles.ritsokiguess.site/diet-wars-1.txt"
diets <- read_table(my_url)

── Column specification ────────────────────────────────────────────────────────
cols(
  Diet = col_double(),
  Status = col_double()
)

diets

10. Change the values in both the Diet and Status columns to be text that describes what they are. Specifically, use LowFat, Medit and LowCarb for the diets, and Completed and DroppedOut for the Status, making sure to match the right description to the right number. Save your updated dataframe.

What we are doing is “recoding” the numbers to the text descriptions, and since we are changing the whole columns, recode_values is what we want. The diets are in the correct order to match the numbers, but the status descriptions correspond to 1 and 0 in that order (just to see whether you are paying attention):

diets %>% 
  mutate(Diet = recode_values(
    Diet,
    1 ~ "LowFat",
    2 ~ "Medit",
    3 ~ "LowCarb"
  )) %>% 
  mutate(Status = recode_values(
    Status,
    0 ~ "DroppedOut",
    1 ~ "Completed"
  )) -> diets1
diets1

I saved this into a new dataframe, because I am going to do this again in a different way shortly.

A couple of things to remember:

• the first input to recode_values is the name of the column where the values are coming from.⁸ This is usually the same as the name of the column where the values are going to, but it doesn’t have to be. For example, I might want to keep the original Diet column and create a column called Diet_name.
• any values, old or new, that are text have to be in quotes, otherwise recode_values will go looking for variables of those names (and probably fail to find them).

If you don’t have recode_values, then you probably need to use case_when and deal with the extra typing:

diets %>% 
  mutate(Diet = case_when(
    Diet == 1 ~ "LowFat",
    Diet == 2 ~ "Medit",
    Diet == 3 ~ "LowCarb"
  )) %>% 
  mutate(Status = case_when(
    Status == 0 ~ "DroppedOut",
    Status == 1 ~ "Completed"
  )) -> diets
diets

Extra: now we can see where the data you read in came from. This was ex1921 from the Sleuth3 package:

Sleuth3::ex1921

This is laid out as a contingency table, with diets in rows and what I called Status in columns. The numbers inside the table are frequencies: how many subjects ended up in that combination of diet and status. I want all the frequencies in one column, which is done with pivot_longer (that you will see later in the course, so don’t worry about the details now):

Sleuth3::ex1921 %>% 
  pivot_longer(-Diet, names_to = "Status", values_to = "frequency")

What this has done is to turn the columns (whether or not a subject completed the study) into one column called Status whose values are the old column names, and the numbers that were in the table are carried along into the column frequency, matched up with the correct Diet and Status. For example, the number 16 is the number of subjects on the Mediterranean diet that dropped out of the study, and you can see that after the pivot-longer it has ended up in the right place.

I wanted to give you a recode_values question, so I thought I would mess these up for you so that you would have to put them back together. I decided to turn the Diet and Status columns into numerical codes, like this:

Sleuth3::ex1921 %>% 
  pivot_longer(-Diet, names_to = "Status", values_to = "frequency") %>% 
  mutate(
    Diet = recode_values(
      Diet,
      "LowFat" ~ 1,
      "Medit" ~ 2,
      "LowCarb" ~ 3,
    ),
    Status = recode_values(
      Status,
      "DroppedOut" ~ 0,
      "Completed" ~ 1
    )
  )

You can see that my recode_values is the exact opposite of yours, where you turned the numbers back into text.

When you saw these data before, you did a grouped bar chart, and this is not quite the right format for that yet: we still have the column frequency, showing the number of subjects in that diet-status combination, and geom_bar requires one subject per row so that it can count them up. An aside before we fix that up: when you already have frequencies, you can still make a bar chart, but you use geom_col instead of geom_bar, and supply the column with the frequencies as y:

Sleuth3::ex1921 %>% 
  pivot_longer(-Diet, names_to = "Status", values_to = "frequency") %>% 
  mutate(
    Diet = recode_values(
      Diet,
      "LowFat" ~ 1,
      "Medit" ~ 2,
      "LowCarb" ~ 3,
    ),
    Status = recode_values(
      Status,
      "DroppedOut" ~ 0,
      "Completed" ~ 1
    )
  ) %>% 
  ggplot(aes(x = Diet, y = frequency, fill = factor(Status))) + geom_col(position = "dodge")

which should be the same bar chart you got before. End of aside.

So the last thing to do is to get rid of the frequency column and produce one row per subject. You know about count, which starts from one row per observation and produces counts of observations per group. It turns out that there is a function to do the opposite (start from counts per group and finish with one row per observation), and it is amusingly called uncount:

Sleuth3::ex1921 %>% 
  pivot_longer(-Diet, names_to = "Status", values_to = "frequency") %>% 
  mutate(
    Diet = recode_values(
      Diet,
      "LowFat" ~ 1,
      "Medit" ~ 2,
      "LowCarb" ~ 3,
    ),
    Status = recode_values(
      Status,
      "DroppedOut" ~ 0,
      "Completed" ~ 1
    )
  ) %>% 
  uncount(frequency)

You can check that there are 10 rows with diet 1 and status 0, and so on. The 312 rows come from the fact that there were 312 people in the study altogether (which you can check by adding up all the numbers in the contingency table we started from.) This is the data set I saved for you.

11. The principal investigator of the study meets with you and says “this is great, but could you change it so that the name of the Mediterranean diet is Mediterranean and not Medit?” Do that, and save your new dataframe.

This time, we are changing only some of the values in a column, so it is best to use replace_values, thus:

diets %>% 
  mutate(Diet = replace_values(
    Diet,
    "Medit" ~ "Mediterranean"
  )) -> diets2

The value of using replace_values is that you only have to specify the thing that changes; everything else stays the same, and you don’t have to specify a default value or anything like that.

To check that it worked, either scroll through your dataframe, or do something like this:

diets2 %>% count(Diet)

and you see that Medit has disappeared.

If you don’t have replace_when, you will need to use case_when again, and deal with the default, which is “whatever the Diet was before”:

diets %>% 
  mutate(Diet = case_when(
    Diet == "Medit" ~ "Mediterranean",
    .default        = Diet
  )) -> diets
diets

and to check:

diets %>% count(Diet)

That worked.

12. In a previous worksheet, you drew a grouped bar chart of these data. Redraw that graph with your new data, and comment on whether it is easier to interpret.

The two columns are very obviously categorical now, so there is no issue with drawing the grouped bar chart:

ggplot(diets, aes(x = Diet, fill = Status)) + geom_bar(position = "dodge")

The diets have been reorganized into alphabetical order, but the greatest fraction of dropouts is still with the low carb diet, and the lowest is still with the low fat diet. I think the interpretation is a lot easier because the diets and statuses have names that we can interpret directly, rather than having to go back to the data description to see what the numbers mean.

Extra: if you like the properly stacked chart better, you can re-do that as well:

ggplot(diets, aes(x = Diet, fill = Status)) + geom_bar(position = "fill")

and you see that the proportion of dropouts on the low fat diet is actually noticeably less than on the other two diets.

Counting seabirds

Each year, bird experts associated with the Kodiak National Wildlife Refuge in Alaska count the number of seabirds of different types on the water in each of four different bays in the area. This is done by drawing (on a map) a number of straight-line “transects”, then driving a boat along each transect and counting the number and type of birds visible within a certain distance of the boat.

Variables of interest are:

• the Year of observation
• the Transect number
• Temp: the temperature
• ObservCond: visibility, from Average up to Ideal
• Bay the name of the bay
• bird: an abbreviation for the species of bird observed
• count: how many of that type of bird were observed (in that year, bay, transect).

The data are in http://datafiles.ritsokiguess.site/seabird_long.csv.

13. Read in and display some of the data.

my_url <- "http://datafiles.ritsokiguess.site/seabird_long.csv"
seabird <- read_csv(my_url)

Rows: 56895 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (4): ObservCond, Bay, ObservCondFactor3, bird
dbl (5): Year, Site, Transect, Temp, count

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

seabird

There are over 56,000 observations, and the variables named (plus some others which we ignore). The pager displays a maximum of 1000 pages of 10 rows each, so you might think there are exactly 10,000 observations, but there are only this many displayed. The actual number of rows is at the top.

Extra: some more details:

• Transect reference. This talks about measuring plant cover on land, but the transects in our data are on water (so you cannot use a measuring tape!). The idea is that by always using the same transects, measurements are consistent from year to year, and by having a set of transects that cover the study area, you should be observing most or all of what you are trying to observe.
• The original paper that the data came from, with a lot more detail about how the observations were made. The aim of this study was to compare seabird numbers from one year to the next, so using the same transects every year was key. It was not important to estimate the number of seabirds in the entire wildlife refuge, but instead to see how the counts in the places you looked changed over time.

14. A data set containing the full bird names and their principal diet is in http://datafiles.ritsokiguess.site/bird_names.csv. Read in and display some of this data set.

The same again, really:

my_url <- "http://datafiles.ritsokiguess.site/bird_names.csv"
bird_names <- read_csv(my_url)

Rows: 15 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (3): bird, bird_name, diet

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

bird_names

There are 15 different bird species observed. Some of them, like the Loon, actually have subspecies, but they are difficult to tell apart (at least to an observer on a boat looking through binoculars⁹).

15. It is awkward to read the bird abbreviations in the first dataframe. Create and save a new dataframe that has the full bird names as well as the number that were observed and all the other information.

This is using the small dataframe of bird names as a lookup table (as we did with the product codes and product names in lecture). The strategy is to use the lookup table second in a left-join. Note that both dataframes have a column called bird with the bird abbreviations, and that is the only column they have in common. Hence there is no need for a join_by:

seabird %>% left_join(bird_names) -> birds_all

Joining with `by = join_by(bird)`

birds_all

Scroll across to see that you have a column containing the bird’s full name (in bird_name) and also its diet. The other way when you have a lot of columns is to instead use glimpse:

glimpse(birds_all)

Rows: 56,895
Columns: 11
$ Year              <dbl> 1986, 1986, 1986, 1986, 1986, 1986, 1986, 1986, 1986…
$ Site              <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, …
$ Transect          <dbl> 101, 101, 101, 101, 101, 101, 101, 101, 101, 101, 10…
$ Temp              <dbl> 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.…
$ ObservCond        <chr> "Good", "Good", "Good", "Good", "Good", "Good", "Goo…
$ Bay               <chr> "Uyak", "Uyak", "Uyak", "Uyak", "Uyak", "Uyak", "Uya…
$ ObservCondFactor3 <chr> "Good", "Good", "Good", "Good", "Good", "Good", "Goo…
$ bird              <chr> "BAGO", "BLSC", "COME", "COMU", "CORM", "HADU", "HOG…
$ count             <dbl> 0, 0, 0, 11, 0, 0, 0, 1, 4, 5, 0, 0, 0, 0, 0, 0, 0, …
$ bird_name         <chr> "Barrow's Goldeneye", "Black Scoter", "Common Mergan…
$ diet              <chr> "shellfish", "shellfish", "fish", "fish", "fish", "s…

and you see the full bird names and their diets at the end.

If you want to be clear about how the join is happening, you can specify the column to join by explicitly:¹⁰

seabird %>% left_join(bird_names, join_by("bird")) 

I would have accepted either version of this left_join if it were on an assignment, but in both cases you should explain why you did it the way you did:

• without the join_by: say that the only column that seabird and bird_names have in common is bird, and so left_join will match that only.
• with the join_by: say (for example) that you wanted to make sure that left_join only matched on the bird column (and not on any others).

Your reader might appreciate your including the join_by so that it is clear what you are doing.

16. Which three species of bird were seen the most often altogether?

This is a group-by and summarize. Since we’re summing up over years, bays, and transects, we don’t need to specify any of them in the group-by: just group by species, and summarize by summing the counts. Use your dataframe with the full bird names in it. The last line shows the three rows of the result with the highest total count:

birds_all %>% 
  group_by(bird_name) %>% 
  summarize(total_count = sum(count)) %>% 
  slice_max(total_count, n = 3)

Common Murre, by a long way, then long-tailed duck and black scoter.

Another way to do this is to arrange the total counts in descending order, and pick off the top three.

17. Make a graph that shows the trend in total counts of each bird species over time. Think about what would make the most appealing graph to understand the time trends.

The first step is to work out the totals over species and year, which will make the raw material for our graph:¹¹

birds_all %>% group_by(bird_name, Year) %>% 
  summarize(total = sum(count)) -> d

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by bird_name and Year.
ℹ Output is grouped by bird_name.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(bird_name, Year))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

d

So now we have total and year (quantitative) and species (categorical). This suggests a scatterplot with the species indicated by colour. Joining the points by lines often works well on a time plot:

ggplot(d, aes(x = Year, y = total, colour = bird_name)) + geom_point() + geom_line()

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_line()`).

The problem here is that some of the birds are much more commonly seen than others. That is the Common Murre at the top, and pretty much everything else mixed up at the bottom, which makes it hard to see what is going on. Also, it’s hard to tell that many colours apart.

Two ideas, then:

• instead of using colours for the species, how about using facets?
• because some birds are much more commonly seen than others, give each facet its own scale:

ggplot(d, aes(x = Year, y = total)) + geom_point() + geom_line() +
  facet_wrap(~ bird_name, scales = "free")

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_line()`).

and now you can see what is going on. There seem to be a few birds reaching a peak around 2000 or just before, and then dropping off, or maybe a bit of up-and-down cyclic behaviour, but there aren’t very many common themes. It would probably be best to talk about each bird species separately, if you were doing that (say in a report about the bird observation counts).

Footnotes

1. R has no way to know what a logical order is for these data, unless we tell it. Sometimes, the order that the categories first appear in the data is a sensible one. Here, Back appears first, Middle second, and Front third, which is a sensible order, so, as we will see later, you could use fct_inorder to make the categories come out in that order, which makes more sense than alphabetical order.

2. You could have another column called Computer in here as well, and R would treat them as different. But that means that in R, case matters.

3. Though height_in_cm, if we had such a column, would match contains.

4. You might want not to do that if being missing is itself informative about heights. The sort of thing I mean is if students with short parents are more likely not to report their parents’ heights. In that case, a mean of parental heights would be biased because you have too few small ones.

5. rm is shorthand for “remove”, and NA is R’s notation for “missing”. Thus the function option na.rm = TRUE means “it is true that you remove the missing values” before calculating the mean. rm comes from the Unix operating system, which lives on as Linux, and is a command for removing files.

6. This is one of the standard recommendations in web literacy: you might not be able to assess an individual claim for accuracy, but you can assess whether a whole website is saying something sensible.

7. tribble is a handy way to construct small dataframes. For this, you lay out things as if they were in a data file, variables in columns and observations in rows.

8. I usually seem to forget this and have to come back and fix it up.

9. If you were able to trap them, examine them, and then release them, it might be a different story.t

10. Note that if you don’t specify what you are joining by, left_join will tell you what it joined by, and you should make sure that this is correct.

11. I am using my “temporary name” d for this dataframe, since I am using it to make a graph and then never again.