Is 'x' properly defined in the code?

[JCGiron]

Hi, I mentioned issues with your code through ResearchGate.

One thing I noticed (and it might totally be that I'm misinterpreting the code), is that in the second function (#FUNCTION FOR BAYESIAN INFERENCES) there is a call for 'x' (#x: node entropies), but the first function does not define 'x' itself.

Should 'tent'(as defined/calculated from the first function) equal 'x' at some point?

---

The other thing is that in the last function you calculate delta, in which the last defined item is deltaA:

--
#DELTA
#calculate delta statistic
#trait: trait vector
delta <- function(trait, tree,lambda0,se,sim,thin,burn) {

ar <- ace(trait,tree,type="discret",method="ML",model="ARD")$lik.anc
x <- nentropy(ar)
mc1 <- emcmc(rexp(1),rexp(1),x,lambda0,se,sim,thin,burn)
mc2 <- emcmc(rexp(1),rexp(1),x,lambda0,se,sim,thin,burn)
mchain <- rbind(mc1,mc2)
deltaA <- mean(mchain[,2]/mchain[,1])

return(deltaA)
}

---

When you run the example, you start by defining deltaA:

deltaA <- delta(trait,tree,lambda0,se,sim,thin,burn)
print(deltaA)

And I was wondering if this 'redefinition' of deltaA would mess up the calculation of the last function entered.

---

Also, I just re-ran the examples and the same error/warning appears twice:
"Warning in sqrt(diag(solve(h))) : NaNs produced"
Is that normal?

--

#CALCULATE DELTA A
deltaA <- delta(trait,tree,lambda0,se,sim,thin,burn)
Warning in sqrt(diag(solve(h))) : NaNs produced
print(deltaA)
[1] 4.783886

#DRAW THE TREE...
par(mfrow=c(1,2))
tree$tip.label <- rep("",ns)
plot(tree,main="SCENARIO A")
ar <- ace(trait,tree,type="discret",method="ML",model="ARD")$lik.anc
Warning in sqrt(diag(solve(h))) : NaNs produced
nodelabels(pie = ar, cex = 1,frame="n")
mtrait <- matrix(0,ncol=3,nrow=ns)
for ( i in 1:ns) {

• mtrait[i,trait[i]] <- 1
  

• }

tiplabels(pie=mtrait,cex=0.5)

trait <- c(2,3,1,3,3,3,3,2,2,3,1,2,1,2,3,1,2,3,1,2)

#CALCULATE DELTA B
deltaB <- delta(trait,tree,lambda0,se,sim,thin,burn)
Warning in sqrt(diag(solve(h))) : NaNs produced
print(deltaB)
[1] 0.5067339

#DRAW THE TREE...
ar <- ace(trait,tree,type="discret",method="ML",model="ARD")$lik.anc
Warning in sqrt(diag(solve(h))) : NaNs produced
plot(tree,main="SCENARIO B")
nodelabels(pie = Re(ar), cex = 1)
mtrait <- matrix(0,ncol=3,nrow=ns)
for ( i in 1:ns) { mtrait[i,trait[i]] <- 1 }
tiplabels(pie=mtrait,cex=0.5)

---

Thanks for taking a look at this.

Cheers,

Jennifer