Ecological Analyses in Julia

Eduardo Santos

doi:10.5281/zenodo.23126524

5 Data wrangling with Tidier

In the original book this is the tidyverse chapter. Julia has a remarkably faithful translation: the Tidier family, a pure-Julia reimplementation of the tidyverse mini-language on top of DataFrames.jl.

using TidierData
using DataFrames
using Statistics

5.1 The data

Throughout the book we use a simulated survey of anurans in 180 ponds across three regions. The data are simulated on purpose: because the generating process is known, every model fitted in Chapter 8 and Chapter 9 can be checked against the truth. The generator lives in scripts/make_data.jl.

using CSV

ponds = CSV.read(joinpath("..", "data", "ponds.csv"), DataFrame)
first(ponds, 5)

5×9 DataFrame

Row	pond_id	region	hydro	area_m2	canopy	temp_c	richness	focal	body_size
	String7	String15	String15	Float64	Float64?	Float64	Int64	Int64	Float64
1	P001	Atlantic	permanent	174.1	70.0	22.4	8	0	34.55
2	P002	Amazon	permanent	161.6	42.6	23.9	12	1	30.36
3	P003	Cerrado	temporary	233.9	57.5	25.2	6	1	34.44
4	P004	Amazon	temporary	174.5	49.4	25.6	13	1	39.3
5	P005	Cerrado	permanent	314.4	30.0	27.6	6	0	43.65

The book is rendered from inside the en/ folder, hence the ".." (the folder above) in this and the following chapters. In your own session, with the book folder open, use joinpath("data", "ponds.csv"), as in the exercise box below.

describe(ponds, :eltype, :nmissing, :mean, :min, :max)

9×6 DataFrame

Row	variable	eltype	nmissing	mean	min	max
	Symbol	Type	Int64	Union…	Any	Any
1	pond_id	String7	0		P001	P180
2	region	String15	0		Amazon	Cerrado
3	hydro	String15	0		permanent	temporary
4	area_m2	Float64	0	240.215	10.0	1981.8
5	canopy	Union{Missing, Float64}	7	49.1422	6.4	95.8
6	temp_c	Float64	0	24.4406	17.0	33.9
7	richness	Int64	0	7.13333	0	31
8	focal	Int64	0	0.45	0	1
9	body_size	Float64	0	35.6777	26.91	44.26

Column	Meaning
`pond_id`	pond identifier
`region`	region (Atlantic, Cerrado, Amazon)
`hydro`	hydroperiod (permanent, temporary)
`area_m2`	pond area (m²)
`canopy`	canopy cover (%)
`temp_c`	water temperature (°C)
`richness`	species richness (count)
`focal`	focal species occurrence (0/1)
`body_size`	mean body size (mm)

Exercises in this chapter

The Your turn boxes work as in Chapter 4. Here the answers are tables and numbers computed from the pond data, so start each session by loading the checker, the packages and the data, with the book folder open in VS Code:

include("exercises/check.jl")
using TidierData, DataFrames, Statistics, CSV
ponds = CSV.read(joinpath("data", "ponds.csv"), DataFrame)

Note

Column names are in English in both versions of the book, so that the code is literally identical in English and Portuguese.

5.2 `@chain`: the pipe

@chain is the backbone of everything that follows. It passes each line’s result as the first argument of the next — exactly the role of |> (or %>%) in R.

@chain ponds begin
    @filter(region == "Amazon")
    @select(pond_id, area_m2, canopy, richness)
    @arrange(desc(richness))
    first(5)
end

5×4 DataFrame

Row	pond_id	area_m2	canopy	richness
	String7	Float64	Float64?	Int64
1	P065	1048.6	18.0	31
2	P044	568.1	31.7	24
3	P067	411.4	30.9	21
4	P085	204.6	58.2	19
5	P073	140.5	17.7	18

In R

ponds %>%
  filter(region == "Amazon") %>%
  select(pond_id, area_m2, canopy, richness) %>%
  arrange(desc(richness)) %>%
  head(5)

The correspondence is line for line. The only syntactic difference is the @ and the begin ... end.

Your turn 5.1 — The largest ponds in the Cerrado

Write a @chain that returns the three largest ponds (by area_m2) in the Cerrado, largest first. Store the result in top3 and check it with check("5.1", top3).

Hint

Three steps, one per line: keep the Cerrado rows, sort by area from largest to smallest, keep the first three. The block above does something very similar.

Solution

top3 = @chain ponds begin
    @filter(region == "Cerrado")
    @arrange(desc(area_m2))
    first(3)
end

3×9 DataFrame

Row	pond_id	region	hydro	area_m2	canopy	temp_c	richness	focal	body_size
	String7	String15	String15	Float64	Float64?	Float64	Int64	Int64	Float64
1	P130	Cerrado	temporary	1981.8	84.9	26.0	9	0	38.53
2	P034	Cerrado	temporary	1126.3	missing	22.9	6	0	33.53
3	P143	Cerrado	temporary	932.7	39.8	25.1	11	0	34.15

check("5.1", top3)

✓ Correct!

5.3 The five verbs

5.3.1 `@filter` — pick rows

@chain ponds begin
    @filter(richness > 8 && hydro == "permanent")
    nrow
end

Note: inside @filter, use && and ||, not & and | as in dplyr.

Your turn 5.2 — Several conditions at once

How many temporary ponds have canopy cover above 60% and richness below the overall median richness? Check the number with check("5.2", n).

Hint

Compute the median first, outside the chain: med = median(ponds.richness). Inside a Tidier macro a bare name means a column, so to use the variable med write !!med.

If you get TypeError: non-boolean (Missing) used in boolean context, some ponds have no canopy value, and missing > 60 is neither true nor false. Remove those rows first with @drop_missing(canopy) (Section 5.4).

Solution

med = median(ponds.richness)

n = @chain ponds begin
    @drop_missing(canopy)
    @filter(hydro == "temporary" && canopy > 60 && richness < !!med)
    nrow
end

check("5.2", n)

✓ Correct!

Dropping ponds without canopy is the right call here: we cannot know whether their canopy is above 60%. The median, though, is computed on all 180 ponds, before dropping anything.

5.3.2 `@select` — pick columns

@chain ponds begin
    @select(pond_id, region, richness)
    first(3)
end

3×3 DataFrame

Row	pond_id	region	richness
	String7	String15	Int64
1	P001	Atlantic	8
2	P002	Amazon	12
3	P003	Cerrado	6

# exclusion and ranges work as in dplyr
@chain ponds begin
    @select(-pond_id, -focal)
    names
end

9-element Vector{String}:
 "region"
 "hydro"
 "area_m2"
 "canopy"
 "temp_c"
 "richness"
 "focal"
 "body_size"
 "pond_id"

5.3.3 `@mutate` — create columns

ponds_t = @chain ponds begin
    @mutate(
        log_area   = log(area_m2),
        canopy_std = (canopy - mean(skipmissing(canopy))) / std(skipmissing(canopy)),
        open       = canopy < 40,
    )
end

@chain ponds_t begin
    @select(pond_id, area_m2, log_area, canopy, canopy_std, open)
    first(5)
end

5×6 DataFrame

Row	pond_id	area_m2	log_area	canopy	canopy_std	open
	String7	Float64	Float64	Float64?	Float64?	Bool?
1	P001	174.1	5.15963	70.0	0.946503	false
2	P002	161.6	5.08512	42.6	-0.296877	false
3	P003	233.9	5.45489	57.5	0.379268	false
4	P004	174.5	5.16192	49.4	0.0116988	false
5	P005	314.4	5.75067	30.0	-0.868651	true

Your turn 5.3 — Richness per 100 m²

Bigger ponds tend to hold more species. Add a column rich_per_100m2 with the number of species per 100 m² of pond, keeping every row and column of ponds. Store the result in ponds_d and check it with check("5.3", ponds_d).

Hint

richness / area_m2 is species per m²; multiply by 100.

Solution

ponds_d = @chain ponds begin
    @mutate(rich_per_100m2 = richness / area_m2 * 100)
end

check("5.3", ponds_d)

✓ Correct!

5.3.4 `@group_by` + `@summarize` — aggregate

@chain ponds begin
    @group_by(region, hydro)
    @summarize(
        n            = n(),
        mean_rich    = mean(richness),
        sd_rich      = std(richness),
        prop_focal   = mean(focal),
    )
    @arrange(region, hydro)
end

GroupedDataFrame with 3 groups based on key: region

First Group (2 rows): region = "Amazon"

Row	region	hydro	n	mean_rich	sd_rich	prop_focal
	String15	String15	Int64	Float64	Float64	Float64
1	Amazon	permanent	26	11.2692	5.99038	0.538462
2	Amazon	temporary	22	10.9091	5.95946	0.318182

⋮

Last Group (2 rows): region = "Cerrado"

Row	region	hydro	n	mean_rich	sd_rich	prop_focal
	String15	String15	Int64	Float64	Float64	Float64
1	Cerrado	permanent	31	4.77419	2.90643	0.580645
2	Cerrado	temporary	35	4.8	2.95854	0.4

Your turn 5.4 — Where is the focal species?

Compute, for each region, the proportion of ponds where the focal species occurs, in a column called prop_focal, sorted from highest to lowest. Check it with check("5.4", by_region).

Hint

focal is 0 or 1, so its mean is the proportion of ponds where it is 1. Sorting from highest to lowest is @arrange(desc(...)).

Solution

by_region = @chain ponds begin
    @group_by(region)
    @summarize(prop_focal = mean(focal))
    @arrange(desc(prop_focal))
end

check("5.4", by_region)

✓ Correct!

5.3.5 `@arrange` — sort

@chain ponds begin
    @arrange(region, desc(richness))
    @select(region, pond_id, richness)
    first(6)
end

6×3 DataFrame

Row	region	pond_id	richness
	String15	String7	Int64
1	Amazon	P065	31
2	Amazon	P044	24
3	Amazon	P067	21
4	Amazon	P085	19
5	Amazon	P073	18
6	Amazon	P164	18

5.4 Missing data

The canopy column has missing values on purpose.

sum(ismissing.(ponds.canopy))

@chain ponds begin
    @drop_missing(canopy)
    nrow
end

# aggregating while ignoring missings
@chain ponds begin
    @group_by(region)
    @summarize(mean_canopy = mean(skipmissing(canopy)))
end

3×2 DataFrame

Row	region	mean_canopy
	String15	Float64
1	Atlantic	51.7125
2	Amazon	45.313
3	Cerrado	49.327

Warning

mean(canopy) returns missing if any value is missing — the same behaviour as R’s mean() without na.rm = TRUE. Be explicit with skipmissing.

Your turn 5.5 — Where are the gaps?

Before dropping missing values, it is worth knowing where they are: if they concentrate in one region, dropping them changes that region more than the others. Count, per region, how many ponds have no canopy value, in a column called n_missing. Check it with check("5.5", gaps).

Hint

ismissing.(canopy) gives true for each missing value, and sum counts the trues. Both work inside @summarize.

Solution

gaps = @chain ponds begin
    @group_by(region)
    @summarize(n_missing = sum(ismissing.(canopy)))
end

check("5.5", gaps)

✓ Correct!

Seven missing values spread over three regions: no region is hit much harder than the others, so dropping them is unlikely to bias comparisons between regions.

5.5 Long and wide formats

long = @chain ponds begin
    @select(pond_id, region, area_m2, canopy, temp_c)
    @pivot_longer(area_m2:temp_c, names_to = "variable", values_to = "value")
end

first(long, 6)

6×4 DataFrame

Row	pond_id	region	variable	value
	String7	String15	String	Float64?
1	P001	Atlantic	area_m2	174.1
2	P002	Amazon	area_m2	161.6
3	P003	Cerrado	area_m2	233.9
4	P004	Amazon	area_m2	174.5
5	P005	Cerrado	area_m2	314.4
6	P006	Cerrado	area_m2	620.2

@chain long begin
    @pivot_wider(names_from = variable, values_from = value)
    first(3)
end

3×5 DataFrame

Row	pond_id	region	area_m2	canopy	temp_c
	String7	String15	Float64?	Float64?	Float64?
1	P001	Atlantic	174.1	70.0	22.4
2	P002	Amazon	161.6	42.6	23.9
3	P003	Cerrado	233.9	57.5	25.2

5.6 Joins

# an auxiliary table with region-level information
regions = DataFrame(
    region     = ["Atlantic", "Cerrado", "Amazon"],
    biome_type = ["forest", "savanna", "forest"],
    mean_rain  = [1400, 1500, 2300],
)

@chain ponds begin
    @left_join(regions, region)
    @select(pond_id, region, biome_type, mean_rain, richness)
    first(5)
end

5×5 DataFrame

Row	pond_id	region	biome_type	mean_rain	richness
	String7	String15	String?	Int64?	Int64
1	P001	Atlantic	forest	1400	8
2	P002	Amazon	forest	2300	12
3	P003	Cerrado	savanna	1500	6
4	P004	Amazon	forest	2300	13
5	P005	Cerrado	savanna	1500	6

Your turn 5.6 — Forest versus savanna

Using the regions table above, compute the mean richness of forest and savanna ponds, in a column called mean_rich. Check it with check("5.6", by_biome).

Hint

Join first, so every pond gets its biome_type; then group by biome_type and summarize.

Solution

by_biome = @chain ponds begin
    @left_join(regions, region)
    @group_by(biome_type)
    @summarize(mean_rich = mean(richness))
end

check("5.6", by_biome)

✓ Correct!

5.7 `DataFrames.jl` underneath

Tidier is a convenience layer. It is worth knowing the native form, which is what you will meet in most of the ecosystem’s documentation:

combine(
    groupby(ponds, :region),
    :richness => mean => :mean_rich,
    :richness => length => :n,
)

3×3 DataFrame

Row	region	mean_rich	n
	String15	Float64	Int64
1	Atlantic	6.59091	66
2	Amazon	11.1042	48
3	Cerrado	4.78788	66

In Tidier, this would be @group_by(region) followed by @summarize(mean_rich = mean(richness), n = n()): the same idea as the @group_by/@summarize above, but grouped by region only. :region is a Symbol, Julia’s way of naming a column, and :richness => mean => :mean_rich reads “take richness, apply mean and call the result mean_rich”. Use Tidier to write, and recognize the native form to read.

Your turn 5.7 — The same, natively

Rewrite your answer to 5.4 with groupby, combine and sort from DataFrames.jl, without Tidier macros. Check it with check("5.7", by_region_native).

Hint

sort(df, :col, rev = true) sorts from highest to lowest.

Solution

by_region_native = sort(
    combine(groupby(ponds, :region), :focal => mean => :prop_focal),
    :prop_focal, rev = true)

check("5.7", by_region_native)

✓ Correct!

Challenges

These have no checker. They are about what the numbers do — and do not — tell you.

In 5.4 the Cerrado comes first, but the three proportions only range from 0.42 to 0.48, each from 48 to 66 ponds. Open scripts/make_data.jl: does region affect the focal species in the simulation? What does that tell you about reading a ranking off a summary table?
In 5.6 forest ponds are richer than savanna ponds. But “forest” here means Atlantic plus Amazon. Compute mean richness per region and decide whether “forest” is a useful grouping for this question.
Use @pivot_longer to put area_m2, canopy and temp_c in one column, then compute the mean of each variable per region in a single pipeline. How do you deal with the missing canopy values?

5.1 The data

5.2 @chain: the pipe

5.3 The five verbs

5.3.1 @filter — pick rows

5.3.2 @select — pick columns

5.3.3 @mutate — create columns

5.3.4 @group_by + @summarize — aggregate

5.3.5 @arrange — sort

5.4 Missing data

5.5 Long and wide formats

5.6 Joins

5.7 DataFrames.jl underneath

Challenges

5.2 `@chain`: the pipe

5.3.1 `@filter` — pick rows

5.3.2 `@select` — pick columns

5.3.3 `@mutate` — create columns

5.3.4 `@group_by` + `@summarize` — aggregate

5.3.5 `@arrange` — sort

5.7 `DataFrames.jl` underneath