Merge pull request #206 from EcoJulia/fix/trophiclevel

Fix trophic level

Project.toml

@@ -2,7 +2,7 @@ name = "EcologicalNetworks"
 uuid = "f03a62fe-f8ab-5b77-a061-bb599b765229"
 authors = ["Timothée Poisot <[email protected]>"]
 repo = "https://github.com/EcoJulia/EcologicalNetworks.jl"
-version = "0.5.1"
+version = "0.5.2"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"

docs/make.jl

@@ -30,7 +30,8 @@ makedocs(
             "Nestedness" => "properties/nestedness.md",
             "Motifs" => "properties/motifs.md",
             "Centrality and paths" => "properties/paths.md",
-            "Overlap and similarity" => "properties/overlap.md"
+            "Overlap and similarity" => "properties/overlap.md",
+            "Food web measures" => "properties/foodwebs.md"
         ],
         "Advanced information" => [
             "Beta-diversity" => "properties/betadiversity.md",

docs/src/properties/foodwebs.md

@@ -0,0 +1,7 @@
+# Food web measures
+
+```@docs
+distance_to_producer
+trophic_level
+omnivory
+```

src/EcologicalNetworks.jl

@@ -162,10 +162,7 @@ export KGL01, KGL02, KGL03, KGL04, KGL05, KGL06, KGL07, KGL08, KGL09, KGL10,
 
 # Food webs
 include(joinpath(".", "foodwebs/trophiclevels.jl"))
-export fractional_trophic_level, trophic_level
-
-include(joinpath(".", "foodwebs/omnivory.jl"))
-export omnivory
+export distance_to_producer, trophic_level, omnivory
 
 include(joinpath(".", "information/entropy.jl"))
 export entropy, make_joint_distribution, mutual_information, conditional_entropy,

src/foodwebs/trophiclevels.jl

@@ -1,24 +1,107 @@
-function fractional_trophic_level(N::T) where {T<:UnipartiteNetwork}
-  Y = nodiagonal(N)
-  producers = keys(filter(spedeg -> spedeg.second == 0, degree(Y; dims=1)))
-  sp = shortest_path(Y)
-  prod_id = findall(isequal(0), vec(sum(sp; dims=2)))
-  return Dict(zip(species(Y; dims=1), maximum(sp[:,prod_id]; dims=2).+1))
+"""
+    distance_to_producer(N::UnipartiteNetwork; f=maximum)
+
+Returns the distance from a primary producert of every species in the food
+web. Primary producers have a trophic level of 1, after that, it's complicated
+(see *e.g.* Williams & Martinez 2004, Thompson et al. 2007).
+
+Post (2002) notes that the trophic level can be obtained from the maximum,
+mean, or minimum distance to a producer. Given that consumers may be connected
+to more than one producer, one might argue that the *mode* of this distribution
+is also relevant.
+
+In favor of the minimum, one can argue that most energy transfer should happen
+along short chains; but imagining a consumer atop a chain of length 5, also
+connected directly to a producer, the minimum would give it a trophic level
+of 2, hiding its position at the top of the food web.
+
+In favor of the maximum, one can argue that the higher chains give a better
+idea of how far energy coming from the bottom of the food web can go. This
+is a strong indication of how *vertically* diverse it is.
+
+In this package, the keyword `f` defaults to `maximum`. Using any other
+function, as long as it accepts an array (of chain distances) and return a
+scalar, will work; `minimum` and `Statistics.mean` are natural alternatives,
+as is `StatsBase.mode` or `Statistics.median`.
+
+The chain length from any species to any producer is taken as the shortest
+possible path between these two species, as identified by the Dijkstra
+algorithm.
+"""
+function distance_to_producer(N::T; f=maximum) where {T<:UnipartiteNetwork}
+    Y = nodiagonal(N)
+    paths = dijkstra(Y)
+    consumers = collect(keys(filter(p -> iszero(p.second), degree(Y; dims=1))))
+    filter!(int -> int.to ∈ consumers, paths)
+    tl = Dict{eltype(species(Y)),Float64}()
+    for sp in species(Y)
+        sp_paths = filter(int -> isequal(sp)(int.from), paths)
+        chain_lengths = [int.weight for int in sp_paths]
+        tl[sp] = isempty(chain_lengths) ? 1.0 : f(chain_lengths) + 1.0
+    end
+    return tl
 end
 
-function trophic_level(N::T) where {T<:UnipartiteNetwork}
-  TL = fractional_trophic_level(N)
-  Y = nodiagonal(N)
-  D = zeros(Float64, size(Y))
-  ko = degree_out(Y)
+"""
+    trophic_level(N::UnipartiteNetwork; kwargs...)
+
+Returns the *fractional* trophic level (after Odum & Heald 1975)
+of species in a binary unipartite network. The trophic level is
+calculated after Pauly & Christensen (1995), specifically as TLᵢ = 1 +
+∑ⱼ(TLⱼ×DCᵢⱼ)/∑ⱼ(DCᵢⱼ), wherein TLᵢ is the trophic level
+of species i, and DCᵢⱼ is the proportion of species j in the diet of
+species i. Note that this function is calculated on a network where the
+diagonal (i.e. self loops) are removed.
+
+This function uses a *pseudo*-inverse to solve the linear system
+described above *iff* the determinant of the diet matrix is 0 (it is
+non-invertible). Otherwise, it will use an exact inverse.
+"""
+function trophic_level(N::UnipartiteNetwork)
+    # Get the degree as a vector, ordered the same way as species
+    kout = EcologicalNetworks.degree_out(N)
+    𝐤 = [kout[s] for s in species(N)]
+    # Adjacency matrix to solve the TL
+    𝐀 = adjacency(N)
+    # Diet matrix
+    𝐃 = .-(𝐀 ./ 𝐤)
+    replace!(𝐃, NaN => -0.0)
+    𝐃[diagind(𝐃)] .= 1.0 .- 𝐃[diagind(𝐃)]
+    # Solve with the inverse matrix
+    𝐛 = ones(eltype(𝐃), richness(N))
+    invfunc = iszero(det(𝐃)) ? pinv : inv
+    return Dict(zip(species(N), invfunc(𝐃) * 𝐛))
+end
+
+"""
+    omnivory(N::T) where {T <: UnipartiteNetwork}
+
+Returns a vector of length richness(N) with an index of consumption across
+trophic levels. Note we explicitly assume that diet contribution is spread
+evenly across prey species, in proportion to their degree. Omnivory is
+measured based on the output of `trophic_level`.
+
+#### References
+- Christensen, Villy, and Daniel Pauly. "ECOPATH II—a software for balancing
+  steady-state ecosystem models and calculating network characteristics."
+  Ecological modelling 61, no. 3-4 (1992): 169-185.
+
+"""
+function omnivory(N::T) where {T<:UnipartiteNetwork}
+    # Get the trophic level as an array
+    TL = trophic_level(N)
+
+    # Prepare the omnivory index dictionary
+    OI = Dict(zip(species(N), fill(0.0, richness(N))))
 
-  # inner loop to avoid dealing with primary producers
-  for i in 1:richness(Y)
-    if ko[species(Y)[i]] > 0.0
-      D[i,:] = Y[i,:]./ko[species(Y)[i]]
+    # Loop
+    for sp in species(N)
+        if !isempty(N[sp,:])
+            TLj = [TL[prey] for prey in N[sp,:]]
+            # Note that we don't need the normalisation here, because ∑Dᵢⱼ=1
+            OI[sp] = sum((TLj .- mean(TLj)).^2.0)
+        end
     end
-  end
 
-  # return
-  return Dict(zip(species(N), 1 .+ D * [TL[s] for s in species(Y)]))
+    return OI
 end

test/community/foodwebs.jl

@@ -1,27 +1,26 @@
 module TestFoodWebs
-  using Test
-  using EcologicalNetworks
+using Test
+using EcologicalNetworks
 
-  # test trophic levels on motifs
-  s1 = unipartitemotifs()[:S1]
-  s1_ftl = fractional_trophic_level(s1)
-  s1_tl = trophic_level(s1)
-  @test maximum(values(s1_ftl)) == 3
-  @test maximum(values(s1_tl)) == 3.0
-  @test s1_ftl[:a] == 3
-  @test s1_ftl[:b] == 2
-  @test s1_ftl[:c] == 1
+# test trophic levels on motifs
+s1 = unipartitemotifs()[:S1]
+s1_ftl = trophic_level(s1)
+s1_tl = distance_to_producer(s1)
+@test maximum(values(s1_ftl)) ≈ 3
+@test maximum(values(s1_tl)) ≈ 3.0
+@test s1_ftl[:a] ≈ 3
+@test s1_ftl[:b] ≈ 2
+@test s1_ftl[:c] ≈ 1
 
-  s2 = unipartitemotifs()[:S2]
-  @test minimum(values(fractional_trophic_level(s1))) == 1.0
-  @test minimum(values(trophic_level(s1))) == 1.0
+s2 = unipartitemotifs()[:S2]
+@test minimum(values(trophic_level(s1))) ≈ 1.0
+@test minimum(values(distance_to_producer(s1))) == 1.0
 
-  # test on a quantitative bipartite network
-  @test web_of_life("A_HP_002") |>
-    x -> convert(BinaryNetwork, x) |> x -> convert(UnipartiteNetwork, x) |>
-    trophic_level |> values |> minimum |> x -> x == 1.0
-  @test web_of_life("A_HP_002") |>
-    x -> convert(BinaryNetwork, x) |> x -> convert(UnipartiteNetwork, x) |>
-    trophic_level |> values |> maximum |> x -> x == 2.0
+# Test for omnivory
+for s in species(s1)
+    if iszero(EcologicalNetworks.degree_out(s1)[s])
+        @test iszero(omnivory(s1)[s])
+    end
+end
 
 end