Network with undefined branches makes phylolm fail

[pbastide]

This issue was raised by Edward A Myers.

If the entry network has some undefined branch lengths, then phylolm fails with an un-informative message:

# Terminal branche A is not defined
net_str= "(A,((B:1,#H1:1::0.1):1,(C:1,(D:1)#H1:1::0.9):1):0.5);"
net = readTopology(net_str);

# phylolm fails
dfr = DataFrame(trait = [11.6,8.1,10.3,9.1], tipNames = ["A","B","C","D"]);
phylolm(@formula(trait ~ 1), dfr, net);

ERROR: PosDefException: matrix is not positive definite; Cholesky factorization failed.
Stacktrace:
  [1] checkpositivedefinite
    @ /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/factorization.jl:18 [inlined]
  [2] cholesky!(A::Hermitian{Float64, Matrix{Float64}}, ::Val{false}; check::Bool)
    @ LinearAlgebra /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:266
  [3] cholesky!(A::Matrix{Float64}, ::Val{false}; check::Bool)
    @ LinearAlgebra /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:298
  [4] #cholesky#143
    @ /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:394 [inlined]
  [5] cholesky (repeats 2 times)
    @ /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:394 [inlined]
  [6] pgls(X::Matrix{Float64}, Y::Vector{Float64}, V::PhyloNetworks.MatrixTopologicalOrder; nonmissing::BitVector, ind::Vector{Int64})
    @ PhyloNetworks ~/.julia/packages/PhyloNetworks/Wtq8B/src/traits.jl:1741
  [7] phylolm(::PhyloNetworks.BM, X::Matrix{Float64}, Y::Vector{Float64}, net::HybridNetwork, reml::Bool; nonmissing::BitVector, ind::Vector{Int64}, kwargs::Base.Pairs{Symbol, Union{Missing, Float64}, Tuple{Symbol, Symbol}, NamedTuple{(:startingValue, :fixedValue), Tuple{Float64, Missing}}})
    @ PhyloNetworks ~/.julia/packages/PhyloNetworks/Wtq8B/src/traits.jl:1689
  [8] #phylolm#632
    @ ~/.julia/packages/PhyloNetworks/Wtq8B/src/traits.jl:1677 [inlined]
  [9] phylolm(f::FormulaTerm{Term, ConstantTerm{Int64}}, fr::DataFrame, net::HybridNetwork; model::String, tipnames::Symbol, no_names::Bool, reml::Bool, ftolRel::Float64, xtolRel::Float64, ftolAbs::Float64, xtolAbs::Float64, startingValue::Float64, fixedValue::Missing, withinspecies_var::Bool, y_mean_std::Bool)
    @ PhyloNetworks ~/.julia/packages/PhyloNetworks/Wtq8B/src/traits.jl:2425
 [10] phylolm(f::FormulaTerm{Term, ConstantTerm{Int64}}, fr::DataFrame, net::HybridNetwork)
    @ PhyloNetworks ~/.julia/packages/PhyloNetworks/Wtq8B/src/traits.jl:2348
 [11] top-level scope
    @ REPL[55]:1

This is because the missing branch lengths are silently replaced by a -1, which causes the network matrix to be ill-defined:

# Branch length is set to -1 by default
net.edge[1]

PhyloNetworks.Edge:
 number:1
 length:-1.0
 attached to 2 node(s) (parent first): -2 1

# This causes a wrong vcv matrix
V = vcv(net)

4×4 DataFrame
 Row │ A        B        C        D       
     │ Float64  Float64  Float64  Float64 
─────┼────────────────────────────────────
   1 │    -1.0      0.0      0.0     0.0
   2 │     0.0      2.5      0.5     0.6
   3 │     0.0      0.5      2.5     1.4
   4 │     0.0      0.6      1.4     3.14

# Which indeed has negative eigen values
eigen(sharedPathMatrix(net)[:Tips])

Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
4-element Vector{Float64}:
 -1.0
  1.3835766421973135
  2.2033313913326182
  4.553091966470069
vectors:
4×4 Matrix{Float64}:
 1.0   0.0         0.0        0.0
 0.0   0.0187359   0.93449   -0.355496
 0.0  -0.784956   -0.206476  -0.584133
 0.0   0.619268   -0.289993  -0.729665

# And for which the cholesky fails
cholesky(sharedPathMatrix(net)[:Tips])

ERROR: PosDefException: matrix is not positive definite; Cholesky factorization failed.
Stacktrace:
 [1] checkpositivedefinite
   @ /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/factorization.jl:18 [inlined]
 [2] cholesky!(A::Hermitian{Float64, Matrix{Float64}}, ::Val{false}; check::Bool)
   @ LinearAlgebra /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:266
 [3] cholesky!(A::Matrix{Float64}, ::Val{false}; check::Bool)
   @ LinearAlgebra /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:298
 [4] #cholesky#143
   @ /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:394 [inlined]
 [5] cholesky (repeats 2 times)
   @ /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/cholesky.jl:394 [inlined]
 [6] top-level scope
   @ REPL[63]:2

Posible Fixes:

• Set missing branch length to NA instead of -1 (same as ape in R) ? But that might break some other things elsewhere ?
• Check that all the branch lengths are positive before doing any lm, and throw an explicit error otherwise ?
• Any other ideas ?

[cecileane]

Yes, we should check that all edges are non-negative / non-missing (as currently coded). Thanks for the PR!

[cecileane]

should be fixed in v0.15.1 (see #180)