install.packages(c("colorspace", "lvplot", "patchwork", "janitor", "lubridate", "vcd", "ggbeeswarm", "kableExtra"))ETC5521 Tutorial 7
Making comparisons between groups and strata
🎯 Objectives
These are exercises so that you can make some numerical and graphical comparisons for various datasets and help to think about the comparisons being made.
🔧 Preparation
- The reading for this week is Wilke (2019) Chapters 9, 10.2-4, 11.2.
- Complete the weekly quiz, before the deadline!
- Install the following R-packages if you do not have them already:
- Open your RStudio Project for this unit, (the one you created in week 1,
ETC5521). Create a.qmddocument for this weeks activities.
📥 Exercises
Exercise 1: Melbourne daily maximum temperature
The csv file melb_temp_2023-09-08.csv contains data on the daily maximum temperature from 1970 to 2023 collected from the weather station at Melbourne Airport. Use this to answer the following questions, with the additional information that in Australia:
- Summer is from the beginning of December to the end of February,
- Autumn is from the beginning of March to the end of May,
- Winter is from the beginning of June to the end of August, and
- Spring is from the beginning of September to the end of November.
- There are four plots below. Write the code to make them yourself. Then think about the three questions (i), (ii) or (iii) below.
- Are there any winters where the daily maximum temperature is different to winter in other years?
- What is the general pattern of maximum daily temperatures in winter?
- Is there evidence that winters in Melbourne are getting warmer?
Which plot best matches each question? If none of them work, for any particular question, make an alternative plot. Also, if any of the plots don’t help answer any of the questions, think about a question that they might answer.
- Make a transformation of the data and a new plot with this variable, that will allow a more direct comparison to answer question (iii).
The data can be read and processed using this code:
melb_df <- read_csv("https://raw.githubusercontent.com/numbats/ddde/main/data/melb_temp_2023-09-08.csv") |>
clean_names() |>
rename(temp = maximum_temperature_degree_c) |>
dplyr::filter(!is.na(temp)) |>
dplyr::select(year, month, day, temp) |>
mutate(
date = as.Date(paste(year, month, day, sep = "-")))
Solution
- Code for the plots is:
# Plot (A)
winter_df <- melb_df |>
dplyr::filter(month %in% c("06", "07", "08")) |>
mutate(winter_day = as.numeric(
difftime(date,
ymd(paste0(year, "-06", "-01")),
units="days")))
p_a <- ggplot(winter_df, aes(x=winter_day, y=temp)) +
geom_point(
data = dplyr::select(winter_df, -year),
color = "gray", size = 0.1
) +
geom_line() +
facet_wrap(~year) +
labs(
tag = "(A)", x = "",
y = "Temperature (°C)"
) +
scale_x_continuous("", breaks = c(0, 31, 62, 92), labels=c("J", "J", "A", "S"))
# Plot (B)
p_b <- ggplot(winter_df, aes(x=year, y=temp)) +
#geom_boxplot() +
geom_jitter(width=0.1, height=0, alpha=0.2) +
stat_summary(fun="mean", geom="point", colour="orangered", size=3) +
geom_smooth(colour="red") +
labs(
tag = "(B)", x = "Year",
y = "Temperature (°C)"
)
# Plot (C)
p_c <- winter_df |>
mutate(year = fct_reorder(as.factor(year), temp)) |>
ggplot(aes(temp, year)) +
# geom_boxplot(aes(color = as.numeric(as.character(year)))) +
geom_quasirandom(aes(x=as.factor(year), y=temp,
colour = as.numeric(as.character(year)))) +
labs(
tag = "(C)", x = "Year",
y = "Temperature (°C)",
color = "Year"
) +
scale_color_continuous_divergingx(mid = 1995) +
coord_flip()
# Plot (D)
p_d <- winter_df |>
mutate(pre1995 = ifelse(year < 1995, "< 1995", "> 1995")) |>
ggplot(aes(x=pre1995, y=temp)) +
geom_lv(aes(fill=after_stat(LV))) +
labs(
tag = "(D)", y = "Temperature (°C)",
x = "Time Period"
) +
scale_fill_discrete_divergingx(palette = "Fall")
p_a + p_b + p_c + p_d +
plot_layout(ncol=1, heights=c(5,3,5,3))Plot A is the only one that allows examining the pattern for winter each year, which helps to address questions (i) and (ii).
In answer to (i) I don’t see any year that is especially different from any other year.
There is a lot of variability in the patterns for winter. For most years it appears to be fairly flat with an increase in August. However, some years temperature remain cool in August.
To better study the winter pattern it might be better to use a smoother model instead of connecting the day-to-day measurements with a line.
ggplot(winter_df, aes(x=winter_day, y=temp)) +
geom_point(
data = dplyr::select(winter_df, -year),
color = "gray", size = 0.1
) +
geom_smooth(se=F, colour="black") +
facet_wrap(~year) +
labs(
tag = "(E)", x = "",
y = "Temperature (°C)"
) +
scale_x_continuous("", breaks = c(0, 31, 62, 92), labels=c("J", "J", "A", "S"))
From plot E, it’s a bit easier to see that overall the years (grey points) there is the expected pattern of a dip to the coolest part of winter and then a warmup. However, this pattern is very rare in any year. Variability is more typical, some years the temperature is flat across June, July, August, and in some years there is a dramatic warming in August. Some years it is even warmer in July than August.
Both plots C, D address question (iii). Plot C is a bit awkward! Reading the colour trend is easy. The top of the diagram is more red, and the bottom more green. To understand what this means, and how it relates to the question requires thinking about the colour mapping, in relation to the sorting of the y axis. If there was no relationship between median temperature and year, the colouring would be very mixed. That’s not what is seen - the later years (more orange) occur with higher medians. Thus this plot would help us answer question (iii), with there is evidence that winters have been warmer in recent years.
Plot D is much more compact and simple to interpret. The colder winter temperatures seen before 1995 have not been observed after 1995. Also the median is higher in the after 1995 years. It gives evidence for “winters are warmer in recent years” too, but it is too coarse a view with only these two subsets of years. Interestingly, the maximums of the maximum temperatures don’t appear to have changed.
- Using the median or a median of an early period of measurements as a baseline, compute the relative change in temperature. Plot these against year, using yearly median as points with a smoother, or possibly using bars (above and below zero) to display the yearly median.
winter_df <- winter_df |>
mutate(rel_change = (temp - median(temp))/median(temp))
winter_df |>
ggplot(aes(x=year, y=rel_change)) +
geom_hline(yintercept=0, linewidth=3, colour="grey90") +
stat_summary(fun="mean", geom="point") +
geom_smooth(colour="red", se=F) +
labs(
tag = "(F)", x = "Year",
y = "Change in temperature",
color = "Year"
) 
Exercise 2: Evidence of Simpson’s paradox?
Check the following data set for evidence of Simpsons Paradox, in the sense that if group2 == "X" the pass rate is higher.
df <- tribble(
~group1, ~group2, ~result, ~count,
"A", "X", "pass", 100,
"B", "X", "pass", 50,
"C", "X", "pass", 25,
"A", "X", "fail", 10,
"B", "X", "fail", 20,
"C", "X", "fail", 20,
"A", "Y", "pass", 10,
"B", "Y", "pass", 70,
"C", "Y", "pass", 15,
"A", "Y", "fail", 20,
"B", "Y", "fail", 40,
"C", "Y", "fail", 30)
Solution
Overall, if group2 is “X” the pass rate is higher. This is true if the data is also examined separately for each of A, B, C of group1. The rates are slightly different though, with A having a big difference in pass rate, and B having a small difference. So there is evidence of Simpsons Paradox because the proportions differ.
df |> group_by(group1) |> summarise(sum(count))# A tibble: 3 × 2
group1 `sum(count)`
<chr> <dbl>
1 A 140
2 B 180
3 C 90df |> group_by(group2) |> summarise(sum(count))# A tibble: 2 × 2
group2 `sum(count)`
<chr> <dbl>
1 X 225
2 Y 185ggplot(df, aes(x=group2, y=count, fill=result)) +
geom_col(position="fill") +
scale_fill_discrete_divergingx()
doubledecker(xtabs(count ~ group1 + group2 + result, data = df),
gp = gpar(fill = c("grey90", "orangered"))
)
👌 Finishing up
Make sure you say thanks and good-bye to your tutor. This is a time to also report what you enjoyed and what you found difficult.