From 2b744d9d6e0602c9552105427d16eaa7a7a0dadc Mon Sep 17 00:00:00 2001 From: Paul Bastide Date: Wed, 6 Jul 2022 18:13:01 +0200 Subject: [PATCH 1/4] Error when the network has negative or missing branches for phylolm. --- src/traits.jl | 14 ++++++++++++++ test/test_lm.jl | 20 ++++++++++++++++++++ 2 files changed, 34 insertions(+) diff --git a/src/traits.jl b/src/traits.jl index bc357d6c..b8a8e3e2 100644 --- a/src/traits.jl +++ b/src/traits.jl @@ -337,6 +337,7 @@ Returns an object of type [`MatrixTopologicalOrder`](@ref). """ function sharedPathMatrix(net::HybridNetwork; checkPreorder=true::Bool) + checkBranchLengths(net::HybridNetwork) recursionPreOrder(net, checkPreorder, initsharedPathMatrix, @@ -382,6 +383,19 @@ function initsharedPathMatrix(nodes::Vector{Node}, params) return(zeros(Float64,n,n)) end +function checkBranchLengths(net::HybridNetwork) + branches = [e.number for e in net.edge] + undefined = branches[[e.length == -1.0 for e in net.edge]] + negatives = branches[[e.length <= 0 for e in net.edge]] + str = "The variance-covariance matrix of the network is not defined, and the phylogenetic regression cannot be done." + if (length(undefined) > 0) + error(string("Branches ", undefined, " have no length in the network. ", str)) + end + if (length(negatives) > 0) + error(string("Branches ", negatives, " have a negative or zero length in the network. ", str)) + end +end + ############################################################################### """ descendenceMatrix(net::HybridNetwork; checkPreorder=true::Bool) diff --git a/test/test_lm.jl b/test/test_lm.jl index 492e5839..89703eea 100644 --- a/test/test_lm.jl +++ b/test/test_lm.jl @@ -624,6 +624,26 @@ lmSH = phylolm(@formula(trait ~ pred), dfr, net, model="scalingHybrid") # λ so large?? largest γ = 0.056, so λγ = 1.34 is > 1... end +############################################################################### +### Undefined branch lengths +############################################################################### +@testset "Undefined branch length" begin + ## No branch length + net = readTopology("(A,((B:1,#H1:1::0.1):1,(C:1,(D:1)#H1:1::0.9):1):0.5);"); + dfr = DataFrame(trait = [11.6,8.1,10.3,9.1], tipNames = ["A","B","C","D"]); + @test_throws ErrorException("Branches [1] have no length in the network. The variance-covariance matrix of the network is not defined, and the phylogenetic regression cannot be done.") phylolm(@formula(trait ~ 1), dfr, net); + ## Negative branch length + net.edge[1].length = -0.5; + @test_throws ErrorException("Branches [1] have a negative or zero length in the network. The variance-covariance matrix of the network is not defined, and the phylogenetic regression cannot be done.") phylolm(@formula(trait ~ 1), dfr, net); + ## Zero branch length + net.edge[1].length = 0.0; + @test_throws ErrorException("Branches [1] have a negative or zero length in the network. The variance-covariance matrix of the network is not defined, and the phylogenetic regression cannot be done.") phylolm(@formula(trait ~ 1), dfr, net); + ## Non zero branch length + net.edge[1].length = 0.1; + fit = phylolm(@formula(trait ~ 1), dfr, net); + @test loglikelihood(fit) ≈ -6.522868090996417 +end + ############################ ## Against no regressor ########################### From 0100906fda3e6eeb6a13fece09a504827f7d7246 Mon Sep 17 00:00:00 2001 From: Paul Bastide Date: Wed, 6 Jul 2022 19:56:52 +0200 Subject: [PATCH 2/4] zero branch lengths actually allowed, most of the time --- src/traits.jl | 2 +- test/test_lm.jl | 20 ++++++++++++-------- 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/src/traits.jl b/src/traits.jl index b8a8e3e2..ff85b16c 100644 --- a/src/traits.jl +++ b/src/traits.jl @@ -386,7 +386,7 @@ end function checkBranchLengths(net::HybridNetwork) branches = [e.number for e in net.edge] undefined = branches[[e.length == -1.0 for e in net.edge]] - negatives = branches[[e.length <= 0 for e in net.edge]] + negatives = branches[[e.length < 0 for e in net.edge]] str = "The variance-covariance matrix of the network is not defined, and the phylogenetic regression cannot be done." if (length(undefined) > 0) error(string("Branches ", undefined, " have no length in the network. ", str)) diff --git a/test/test_lm.jl b/test/test_lm.jl index 89703eea..f0b67da8 100644 --- a/test/test_lm.jl +++ b/test/test_lm.jl @@ -629,19 +629,23 @@ end ############################################################################### @testset "Undefined branch length" begin ## No branch length - net = readTopology("(A,((B:1,#H1:1::0.1):1,(C:1,(D:1)#H1:1::0.9):1):0.5);"); + net = readTopology("(A:2.5,((B,#H1:1::0.1):1,(C:1,(D:1)#H1:1::0.9):1):0.5);"); dfr = DataFrame(trait = [11.6,8.1,10.3,9.1], tipNames = ["A","B","C","D"]); - @test_throws ErrorException("Branches [1] have no length in the network. The variance-covariance matrix of the network is not defined, and the phylogenetic regression cannot be done.") phylolm(@formula(trait ~ 1), dfr, net); + @test_throws ErrorException("Branches [2] have no length in the network. The variance-covariance matrix of the network is not defined, and the phylogenetic regression cannot be done.") phylolm(@formula(trait ~ 1), dfr, net); ## Negative branch length - net.edge[1].length = -0.5; - @test_throws ErrorException("Branches [1] have a negative or zero length in the network. The variance-covariance matrix of the network is not defined, and the phylogenetic regression cannot be done.") phylolm(@formula(trait ~ 1), dfr, net); + net.edge[2].length = -0.5; + @test_throws ErrorException("Branches [2] have a negative or zero length in the network. The variance-covariance matrix of the network is not defined, and the phylogenetic regression cannot be done.") phylolm(@formula(trait ~ 1), dfr, net); ## Zero branch length - net.edge[1].length = 0.0; - @test_throws ErrorException("Branches [1] have a negative or zero length in the network. The variance-covariance matrix of the network is not defined, and the phylogenetic regression cannot be done.") phylolm(@formula(trait ~ 1), dfr, net); + net.edge[2].length = 0.0; + fit = phylolm(@formula(trait ~ 1), dfr, net); + @test loglikelihood(fit) ≈ -6.245746681512051 ## Non zero branch length - net.edge[1].length = 0.1; + net.edge[2].length = 0.1; fit = phylolm(@formula(trait ~ 1), dfr, net); - @test loglikelihood(fit) ≈ -6.522868090996417 + @test loglikelihood(fit) ≈ -6.2197697662066815 + ## Illicit zero branch length + net.edge[1].length = 0.0; + @test_throws PosDefException phylolm(@formula(trait ~ 1), dfr, net); end ############################ From faf1aa57a21b4c56c6da02e6b1a36e34560a52da Mon Sep 17 00:00:00 2001 From: Paul Bastide Date: Wed, 6 Jul 2022 20:57:23 +0200 Subject: [PATCH 3/4] Bring PosDefException in test scope --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 622be528..153ddc5d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,7 +10,7 @@ using CSV # for reading files using DataFrames using Distributed # parallel in test_correctLik.jl and test_bootstrap.jl using GLM # for coef, nobs, residuals etc. -using LinearAlgebra: norm, diag, logdet # LinearAlgebra.rotate! not brought into scope +using LinearAlgebra: norm, diag, logdet, PosDefException # LinearAlgebra.rotate! not brought into scope using Random using StaticArrays # for rate substitution matrices using Statistics From 79ecc713c20dd5c7741056affcdf1b0a757990c4 Mon Sep 17 00:00:00 2001 From: Cecile Ane Date: Thu, 7 Jul 2022 22:49:14 +0200 Subject: [PATCH 4/4] more specific name, docstring, cosmetic changes --- src/traits.jl | 27 +++++++++++++++++---------- test/test_lm.jl | 17 +++++++++-------- 2 files changed, 26 insertions(+), 18 deletions(-) diff --git a/src/traits.jl b/src/traits.jl index ff85b16c..98a68dd9 100644 --- a/src/traits.jl +++ b/src/traits.jl @@ -337,7 +337,9 @@ Returns an object of type [`MatrixTopologicalOrder`](@ref). """ function sharedPathMatrix(net::HybridNetwork; checkPreorder=true::Bool) - checkBranchLengths(net::HybridNetwork) + check_nonmissing_nonnegative_edgelengths(net, + """The variance-covariance matrix of the network is not defined. + A phylogenetic regression cannot be done.""") recursionPreOrder(net, checkPreorder, initsharedPathMatrix, @@ -383,16 +385,21 @@ function initsharedPathMatrix(nodes::Vector{Node}, params) return(zeros(Float64,n,n)) end -function checkBranchLengths(net::HybridNetwork) - branches = [e.number for e in net.edge] - undefined = branches[[e.length == -1.0 for e in net.edge]] - negatives = branches[[e.length < 0 for e in net.edge]] - str = "The variance-covariance matrix of the network is not defined, and the phylogenetic regression cannot be done." - if (length(undefined) > 0) - error(string("Branches ", undefined, " have no length in the network. ", str)) +""" + check_nonmissing_nonnegative_edgelengths(net, str="") + +Throw an Exception if `net` has undefined edge lengths (coded as -1.0) or +negative edge lengths. The error message indicates the number of the offending +edge(s), followed by `str`. +""" +function check_nonmissing_nonnegative_edgelengths(net::HybridNetwork, str="") + if any(e.length == -1.0 for e in net.edge) + undefined = [e.number for e in net.edge if e.length == -1.0] + error(string("Branch(es) number ", join(undefined,","), " have no length.\n", str)) end - if (length(negatives) > 0) - error(string("Branches ", negatives, " have a negative or zero length in the network. ", str)) + if any(e.length < 0 for e in net.edge) + negatives = [e.number for e in net.edge if e.length < 0.0] + error(string("Branch(es) number ", join(negatives,","), " have negative length.\n", str)) end end diff --git a/test/test_lm.jl b/test/test_lm.jl index f0b67da8..504b5d99 100644 --- a/test/test_lm.jl +++ b/test/test_lm.jl @@ -631,19 +631,20 @@ end ## No branch length net = readTopology("(A:2.5,((B,#H1:1::0.1):1,(C:1,(D:1)#H1:1::0.9):1):0.5);"); dfr = DataFrame(trait = [11.6,8.1,10.3,9.1], tipNames = ["A","B","C","D"]); - @test_throws ErrorException("Branches [2] have no length in the network. The variance-covariance matrix of the network is not defined, and the phylogenetic regression cannot be done.") phylolm(@formula(trait ~ 1), dfr, net); + @test_throws ErrorException("""Branch(es) number 2 have no length. + The variance-covariance matrix of the network is not defined. + A phylogenetic regression cannot be done.""") phylolm(@formula(trait ~ 1), dfr, net); ## Negative branch length net.edge[2].length = -0.5; - @test_throws ErrorException("Branches [2] have a negative or zero length in the network. The variance-covariance matrix of the network is not defined, and the phylogenetic regression cannot be done.") phylolm(@formula(trait ~ 1), dfr, net); - ## Zero branch length + @test_throws ErrorException("""Branch(es) number 2 have negative length. + The variance-covariance matrix of the network is not defined. + A phylogenetic regression cannot be done.""") phylolm(@formula(trait ~ 1), dfr, net); + ## Zero branch length: allowed net.edge[2].length = 0.0; fit = phylolm(@formula(trait ~ 1), dfr, net); @test loglikelihood(fit) ≈ -6.245746681512051 - ## Non zero branch length - net.edge[2].length = 0.1; - fit = phylolm(@formula(trait ~ 1), dfr, net); - @test loglikelihood(fit) ≈ -6.2197697662066815 - ## Illicit zero branch length + net.edge[2].length = 0.1; # back to non-zero + ## Illicit zero length for 1st edge: from root to single "outgroup" taxon net.edge[1].length = 0.0; @test_throws PosDefException phylolm(@formula(trait ~ 1), dfr, net); end