Delta Statsitics running forever

[luisevonderwiese]

Hi!
I have the problem that delta(trait,tree,lambda0,se,sim,thin,burn) is running forever with the attached tree and this trait:

    "0"     "0"     "0"     "0"     "0"     "0"     "1"     "0"     "0"     "1" 
     OE       E     Ger      Da     Ice     Swe     Far     Dut     Blg     Mac 
    "0"     "1"     "0"     "0"     "0"     "0"     "0"     "0"     "0"     "0" 
     SC     Slo     OCS     Cze    Slvk      Po     Rus     Ukr    BelR    Latv 
    "0"     "0"     "?"     "0"     "0"     "0"     "0"     "0"     "0"     "0" 
   Lith      Ir     Wel   ScGae     Bre  Sanskr     Mar     Hin     Guj   Kashm 
    "0"     "0"     "1"     "1"     "0"     "?"     "1"     "1"     "1"     "1" 
Bengali  Nepali   Farsi  Pashto     Alb     Arm 
    "1"     "1"     "1"     "1"     "0"     "0"  

What is going on here? Does this indicate high or low phylogenetic signal?
Thank you very much for your support!
Luise

cognate_ie_compatible_rd.tree.txt

[bbitarello]

I have the same problem. I was able to find where the problem is, but not the solution.

Here is a reproducible example:

library("ape")
source("scripts/code_delta.R")

#SOME PARAMETERS...
lambda0 <- 0.1   #rate parameter of the proposal
se      <- 0.5   #standard deviation of the proposal
sim     <- 10000 #number of iterations
thin    <- 10    #we kept only each 10th iterate
burn    <- 100   #100 iterates are burned-in

#RANDOM TREE:
#same for both examples, only the trait vector varies.
set.seed(25)

ns=20 # nr species
k=8 # nr discrete character states possible

tree <- rtree(ns)
trait<-sample(1:k, ns, replace=T)
deltaA <- delta(trait,tree,lambda0,se,sim,thin,burn)

This will run forever. Why?

Looking at the delta() code, and running line by line:

ar <- ace(trait,tree,type="discret",method="ML",model="ARD")$lik.anc
head(ar)
                  1                2                 3                   4
21 0.57331449861-0i 0.39619261451-0i 0.000092604563+0i 4.566863e-06+0e+00i
22 0.15353848462+0i 0.04345525530-0i 0.222362656675+0i 5.776166e-01-0e+00i
23 0.00008723593+0i 0.00001643572+0i 0.771910095217-0i 1.476295e-01-0e+00i
24 0.00034768848-0i 0.00006462666-0i 0.164177230233+0i 1.313340e-03+0e+00i
25 0.00884309886+0i 0.01086807711+0i 0.000036860471+0i 4.923686e-09+0e+00i
26 0.00923148924+0i 0.00861618361+0i 0.000000222462+0i 4.216449e-10+0e+00i
                  5                   7                  8
21 0.00001922104-0i  0.0009210151467-0i  0.029455479266-0i
22 0.00061247739+0i -0.0000006813249-0i  0.002415181924+0i
23 0.08173394087+0i -0.0013756929057-0i -0.000001500631+0i
24 1.15028765676-0i -0.3161989467451-0i  0.000008404311-0i
25 0.01127825178-0i  0.2166637112836-0i  0.752309995569+0i
26 0.00171083152-0i  0.1680875157688-0i  0.812353756983+0i

 # deletes the complex part whenever it pops up
  if (class(ar[1,1]) == "complex"){
    ar <- Re(ar)
  }

              1            2            3            4            5
21 5.733145e-01 3.961926e-01 9.260456e-05 4.566863e-06 1.922104e-05
22 1.535385e-01 4.345526e-02 2.223627e-01 5.776166e-01 6.124774e-04
23 8.723593e-05 1.643572e-05 7.719101e-01 1.476295e-01 8.173394e-02
24 3.476885e-04 6.462666e-05 1.641772e-01 1.313340e-03 1.150288e+00
25 8.843099e-03 1.086808e-02 3.686047e-05 4.923686e-09 1.127825e-02
26 9.231489e-03 8.616184e-03 2.224620e-07 4.216449e-10 1.710832e-03
               7             8
21  9.210151e-04  2.945548e-02
22 -6.813249e-07  2.415182e-03
23 -1.375693e-03 -1.500631e-06
24 -3.161989e-01  8.404311e-06
25  2.166637e-01  7.523100e-01
26  1.680875e-01  8.123538e-01

x  <- nentropy(ar)
x
         21          22          23          24          25          26
 0.20224170  0.38756261  0.26053716 -0.20020903  0.20286401  0.18948518
         27          28          29          30          31          32
 0.17568399  0.19269650  0.32998488  0.33077020  0.18066524  0.33325460
         33          34          35          36          37          38
 0.34783233  0.33302131  0.33498690  0.16748766  0.32887182  0.19344868
         39
 0.06753229

Everything look ok so far. The next step in the delta() function is the emcmc() function:

mc1    <- emcmc(rexp(1),rexp(1),x,lambda0,se,sim,thin,burn)

This keeps running on and on. Looking into that function line by line:

 usim <- seq(burn,sim,thin)
  gibbs <- matrix(NA,ncol=2,nrow=length(usim))
  p <- 1
usim

> head(usim)
[1] 100 110 120 130 140 150

Next in the function there is the following loop:

 for (i in 1:sim) {
    alpha <- mhalpha(alpha,beta,x,l0,se)
    beta  <- mhbeta(alpha,beta,x,l0,se)

    if (i == usim[p]) {
      gibbs[p,] <- c(alpha,beta)
      p <- p+1
    }
  }

Let's run one iteration:

 alpha <- mhalpha(alpha,beta,x,l0,se)

Runs forever. Looking into mhalpha:

mhalpha <- function(a,b,x,l0,se) {       #metropolis hastings alpha
  a0 <- a
  a1 <- exp(rnorm(1,log(a0),se))

  r  <- min(1, exp(lpalpha(a1,b,x,l0) - lpalpha(a0,b,x,l0) ) )
  i=0 # I added this to understand what's going on
  while (is.na(r) == T) {
        cat(i, "th iteration of while loop","\r") #I added this
    a1 <- exp(rnorm(1,log(a0),se))
 r  <- min(1, exp(lpalpha(a1,b,x,l0) - lpalpha(a0,b,x,l0) ) )

    i=i+1 #I added this
  }

So the reason the function runs forever is because that "while" loop keeps trying to solve the issue of "r" being NA, but it never gets solves because the part of the code that makes is.na(r)==T is lpalpha(a0,b,x,l0), and that part remains unchanged in each iteration. Why is lpalpha(a0,b,x,l0) == NA?

lpalpha <- function(a,b,x,l0) {          #log posterior alpha
  N  <- length(x)
  lp <- N*(lgamma(a+b)-lgamma(a)) - a*(l0-sum(log(x)))
  return(lp)
}

The culprit: sum(log(x). If x<0, log(x) is undefined. Then we get is.na(r)==T and the infinite loop.

What I don't know: why do we get negative x values?

As a recap, x is the output from the nentropy function:

>x
         21          22          23          24          25          26
 0.20224170  0.38756261  0.26053716 -0.20020903  0.20286401  0.18948518
         27          28          29          30          31          32
 0.17568399  0.19269650  0.32998488  0.33077020  0.18066524  0.33325460
         33          34          35          36          37          38
 0.34783233  0.33302131  0.33498690  0.16748766  0.32887182  0.19344868
         39
 0.06753229
> log(x)
        21         22         23         24         25         26         27
-1.5982918 -0.9478779 -1.3450098        NaN -1.5952194 -1.6634445 -1.7390684
        28         29         30         31         32         33         34
-1.6466389 -1.1087085 -1.1063314 -1.7111095 -1.0988485 -1.0560347 -1.0995488
        35         36         37         38         39
-1.0936638 -1.7868456 -1.1120872 -1.6427430 -2.6951495
Warning message:
In log(x) : NaNs produced
> sum(log(x))
[1] NaN
Warning message:
In log(x) : NaNs produced

To make the code keep going in spite of this, a simple chang to sum(log(x), na.rm=T) would do the trick, but because I don't know why nentropy is returning some negative values, I am unsure about this trick.

Hope this helps!

[bbitarello]

By the way, I changed one line of code and the problem is resolved, although I don't quite understand why the ancestral reconstruction is leading to negative probabilities in the first place.

Here it is, in the function delta(). Basically, adding marginal=T to ace():

ar <- ace(trait,tree,type="discrete",method="ML",model="ARD", marginal=T)$lik.anc