Ecological Analyses in Julia

Eduardo Santos

doi:10.5281/zenodo.23126524

4 Introduction to Julia

This chapter introduces the language for readers who already program in R. The focus is on the differences that matter in practice, not an exhaustive sweep of the syntax.

using Statistics
using Random
Random.seed!(42)

TaskLocalRNG()

How the exercises work

Most sections end with a Your turn box. Work on them in VS Code with the book folder open (Section 3.8), and load the checker once per session:

include("exercises/check.jl")

Then pass your answer — a value, or a function you wrote — with the exercise number:

check("4.3", simpson)

check says whether the answer is right and, if not, what it did and what was expected. Each box has a hint and a solution folded away: try first, then peek.

4.1 Vectors and indexing

heights = [12.4, 8.1, 15.7, 9.3, 11.0]

length(heights), first(heights), last(heights)

(5, 12.4, 11.0)

Separating values with commas makes a tuple: a fixed group of values, shown in parentheses. This is how Julia returns several results at once.

Two things work as in R, and one does not:

heights[1]        # indexing starts at 1, as in R
heights[2:4]      # ranges are inclusive at both ends

3-element Vector{Float64}:
  8.1
 15.7
  9.3

# Julia does NOT silently recycle vectors.
# try ... catch captures the error and shows it, instead of stopping the book:
try
    [1, 2, 3] + [1, 2]
catch e
    e
end

DimensionMismatch("a has size (3,), b has size (2,), mismatch at dim 1")

In R

c(1,2,3) + c(1,2) recycles the shorter vector and emits only a warning; Julia raises DimensionMismatch, as the block above shows. The difference is one of policy: R prioritises carrying on, Julia prioritises failing early. Which serves you better depends on your workflow — what matters is knowing which one you are using.

Your turn 4.1 — Drop a bad value

These are the areas (m²) of five ponds. The last one is a data-entry error.

areas = [174.1, 161.6, 233.9, 174.5, 1980.0]

5-element Vector{Float64}:
  174.1
  161.6
  233.9
  174.5
 1980.0

Create clean, holding every pond except the last. Do not type the position 4: your code should still work if more ponds are added.

Hint

Inside square brackets, end means “the last position”, so end - 1 is the one before it.

Solution

clean = areas[1:end-1]
check("4.1", clean)

✓ Correct!

4.2 Broadcasting: the dot

In R, nearly every function is vectorized by default. In Julia you declare vectorization with a dot. This is explicit and works for any function, including the ones you write.

log.(heights)              # elementwise log
heights .^ 2               # elementwise power
heights .> 10              # elementwise comparison

5-element BitVector:
 1
 0
 1
 0
 1

# works identically with your own function
standardize(x, μ, σ) = (x - μ) / σ

standardize.(heights, mean(heights), std(heights))

5-element Vector{Float64}:
  0.3724002529860428
 -1.0833461905048525
  1.4896010119441712
 -0.6770913690655326
 -0.10156370535983013

The dot rule

If you mean “apply this to each element”, add a dot. f.(x), x .+ y, x .== y. An entire chain can be dotted at once with @., a macro (Section 4.8):

@. log(heights) + 1

5-element Vector{Float64}:
 3.517696472610991
 3.0918640616783932
 3.7536607123542622
 3.2300144001592104
 3.3978952727983707

Your turn 4.2 — Log-transform the areas

The area plot in Section 3.8 used a log scale, because pond areas are strongly right-skewed. Compute the base-10 logarithm of every value in areas and store it in log_areas.

Hint

log10(areas) fails. Read the error, then reread the dot rule just above.

Solution

log_areas = log10.(areas)
check("4.2", log_areas)

✓ Correct!

log10 is defined for a single number; the dot applies it to each element.

4.3 Functions

# short form, for one-line functions
richness(abundances) = count(>(0), abundances)

# long form
function shannon(abundances)
    total = sum(abundances)
    p = abundances ./ total
    p = filter(>(0), p)          # log(0) is undefined
    return -sum(p .* log.(p))
end

community = [45, 12, 0, 8, 23, 0, 3]
richness(community), shannon(community)

(5, 1.289250407583511)

>(0) is a function that asks “is this greater than 0?”. count(>(0), x) counts the elements of x for which the answer is yes, and filter(>(0), x) keeps only those.

Your turn 4.3 — Simpson’s index

Write simpson(abundances), returning Simpson’s index $1 - \sum p_i^2$, where $p_i$ is the proportion of individuals belonging to species $i$. The first lines of shannon already compute the proportions.

Before calling check, try simpson([3, 1]) yourself: it should give 0.375.

Hint

p .^ 2 squares each proportion and sum adds them up. Zeros need no special care here: $0^2 = 0$.

Solution

function simpson(abundances)
    p = abundances ./ sum(abundances)
    return 1 - sum(p .^ 2)
end

check("4.3", simpson)

✓ Correct!

4.4 Multiple dispatch

This is the central idea of the language, and the reason Julia packages compose without glue code. A function can have several methods, and Julia picks which one to run based on the types of all the arguments.

# one method for integer counts
diversity(x::Vector{<:Integer}) = "index for count data"

# another for proportional cover
diversity(x::Vector{<:AbstractFloat}) = "index for cover data"

diversity([45, 12, 8]), diversity([0.45, 0.12, 0.08])

("index for count data", "index for cover data")

:: declares the type an argument must have, and Vector{<:Integer} means “a vector of some kind of integer” (<: reads “is a subtype of”). A vector of decimal numbers does not fit the first method, so Julia picks the second.

In practice this means that when someone defines a new ecological data type, mean, plot and fit start working with it as soon as the appropriate methods are defined — without the authors of Statistics or Makie knowing the new type exists.

Your turn 4.4 — One index per site

simpson handles one community. Field data usually come as a site × species matrix, one row per site. Add a second method, simpson(M::Matrix), that returns one index per row. The vector method must keep working.

Hint

eachrow(M) goes through the rows one at a time. Combine it with the dot, or with a comprehension: [f(row) for row in eachrow(M)]. If check complains about the argument types, read its message carefully — it is a common trap.

Solution

simpson(M::Matrix) = simpson.(eachrow(M))

sites = [45 12  0 8;
          5  5  5 5;
         30  0  0 0]
simpson(sites)

3-element Vector{Float64}:
 0.47147928994082844
 0.75
 0.0

check("4.4", simpson)
length(methods(simpson))

✓ Correct!

The matrix method reuses the vector method for each row. One function, two methods, and Julia picks by the argument’s type; run methods(simpson) in the REPL to see both.

4.5 Missing data

measurements = [12.4, missing, 15.7, 9.3]

mean(measurements)              # propagates, as in R

missing

mean(skipmissing(measurements))   # equivalent to na.rm = TRUE

12.466666666666669

In R

NA corresponds to missing. The difference is that missing is a type of its own, and the compiler knows whether a column can contain it — a Vector{Float64} never contains missing, a Vector{Union{Missing,Float64}} may. That eliminates a whole class of silent bugs.

Your turn 4.5 — How much is missing?

Canopy cover was not recorded at every pond:

canopy = [70.0, 42.6, missing, 49.4, missing, 81.2]

6-element Vector{Union{Missing, Float64}}:
 70.0
 42.6
   missing
 49.4
   missing
 81.2

Write missing_summary(x) returning a tuple of two numbers: how many values are missing, and the mean of the rest. For canopy that is (2, 60.8).

Hint

count(f, x) counts the elements for which f is true, and ismissing is such a function. The mean is in the block just above. Return both with return (n, m).

Solution

missing_summary(x) = (count(ismissing, x), mean(skipmissing(x)))

missing_summary(canopy)

(2, 60.8)

check("4.5", missing_summary)

✓ Correct!

4.6 Chaining

# the native pipe |> chains single-argument functions
community |> shannon |> exp    # Hill number of order 1

3.63006446685595

For pipelines with multiple arguments, @chain (from Chain.jl, re-exported by TidierData.jl) plays the role of R’s |>. It is the backbone of Chapter 5.

4.7 Control flow

result = String[]

for (i, a) in enumerate(community)
    if a == 0
        push!(result, "species $i absent")
    elseif a < 10
        push!(result, "species $i rare")
    else
        push!(result, "species $i common")
    end
end

result

7-element Vector{String}:
 "species 1 common"
 "species 2 common"
 "species 3 absent"
 "species 4 rare"
 "species 5 common"
 "species 6 absent"
 "species 7 rare"

String[] makes an empty vector of strings, which push! fills. enumerate gives the position i and the value a of each element. Inside a string, $i inserts the value of i (this is interpolation, like paste0("species ", i) in R).

Loops are fast

Unlike R, you do not need to vectorize for performance. Write the loop that expresses the idea. This changes how null models and simulations are written — a topic that returns in Chapter 11 (see Section 11.6).

Your turn 4.6 — Shannon, one species at a time

Rewrite shannon as shannon_loop, using a for loop instead of broadcasting. It must give the same result as shannon, including for communities with zero counts such as community.

Hint

Start with H = 0.0, go through the counts one by one, turn each into a proportion and subtract p * log(p) from H. If you get NaN, type 0 * log(0) in the REPL.

Solution

function shannon_loop(abundances)
    total = sum(abundances)
    H = 0.0
    for a in abundances
        a == 0 && continue       # absent species contribute nothing
        p = a / total
        H -= p * log(p)
    end
    return H
end

shannon_loop(community), shannon(community)

(1.289250407583511, 1.289250407583511)

check("4.6", shannon_loop)

✓ Correct!

a == 0 && continue is Julia shorthand for if a == 0 continue end.

4.8 Naming conventions

Convention	Meaning
`f!(x)`	modifies `x` in place (`sort!`, `push!`)
`f.(x)`	applies `f` elementwise
`@f x`	a macro, not a function — rewrites code before running
`Type`	types in `CamelCase`
`function_name`	functions in lowercase; underscores are fine when they help readability (`missing_summary`)

Challenges

These have no checker and no single right answer. They are about what your code does — and what it does not tell you.

simpson([-3, 5]) returns a number. Which number, and should it? Change the vector method of simpson so that negative counts raise an informative error with throw(ArgumentError("...")). Does the matrix method inherit the fix?
Communities [10, 4, 3, 3] and [8, 7, 4, 1] have the same richness and the same number of individuals. Compute Shannon and Simpson for both. Do the two indices agree on which community is more diverse? What does this imply for a paper that reports “diversity” as a single number?
Load data/ponds.csv as in Section 3.8 and run missing_summary on ponds.canopy. Does the count agree with what describe(ponds) reports?