Problem with profiles.plot

[nhcooper123]

This won't work for alpha or theta...

yule25.trees<-structure(list(edge = structure(c(26L, 27L, 28L, 29L, 29L, 30L, 
30L, 31L, 31L, 28L, 32L, 33L, 34L, 35L, 35L, 34L, 33L, 36L, 37L, 
37L, 36L, 38L, 39L, 40L, 40L, 39L, 41L, 41L, 42L, 42L, 38L, 32L, 
43L, 43L, 44L, 44L, 27L, 45L, 46L, 46L, 45L, 26L, 47L, 48L, 49L, 
49L, 48L, 47L, 27L, 28L, 29L, 1L, 30L, 2L, 31L, 3L, 4L, 32L, 
33L, 34L, 35L, 5L, 6L, 7L, 36L, 37L, 8L, 9L, 38L, 39L, 40L, 10L, 
11L, 41L, 12L, 42L, 13L, 14L, 15L, 43L, 16L, 44L, 17L, 18L, 45L, 
46L, 19L, 20L, 21L, 47L, 48L, 49L, 22L, 23L, 24L, 25L), .Dim = c(48L, 
2L)), Nnode = 24L, tip.label = c("t1", "t2", "t3", "t4", "t5", 
"t6", "t7", "t8", "t9", "t10", "t11", "t12", "t13", "t14", "t15", 
"t16", "t17", "t18", "t19", "t20", "t21", "t22", "t23", "t24", 
"t25"), edge.length = c(0.5904258688, 0.8249415727, 2.470189261, 
0.932987671, 0.7898975338, 0.1430901373, 0.01988442392, 0.1232057133, 
0.1232057133, 0.4996568724, 0.138501434, 1.739476059, 0.3999097194, 
0.6256328474, 0.6256328474, 1.025542567, 1.180119854, 1.198822986, 
0.386075786, 0.386075786, 0.3431831631, 0.2316611904, 0.3650179037, 
0.6450365148, 0.6450365148, 0.1565813483, 0.8534730703, 0.7328424265, 
0.1206306438, 0.1206306438, 1.241715609, 2.049187285, 0.8543327749, 
0.2658454834, 0.5884872915, 0.5884872915, 0.9436889333, 2.704861689, 
0.5795678831, 0.5795678831, 3.284429572, 3.050576216, 0.1083979538, 
0.1808128489, 1.478757356, 1.478757356, 1.659570205, 1.767968158
)), .Names = c("edge", "Nnode", "tip.label", "edge.length"), class = "phylo", order = "cladewise")

yule25.data<-data.frame(t1 = 2.24, t10= 2.86, t11= 0.576, t12=0.244, t13= 1.82, t14= 1.33,t15= 2.78, t16= 0.985, t17=1.77, t18= 1.65, t19 = 3.8, t2 =-0.531, t20 = 2.7, t21 = 2.06, t22= -3.73, t23= -1.19, t24= -2.29, t25= -0.507, t3 = -0.065, t4 = -0.568, t5 =1.59, t6 =1.48,  t7 =1.48, t8 = 2.45, t9 =0.688)

#Make likelihood function for OU
lik.ou <- make.ou(yule25.trees, yule25.data, control = list(method="pruning", backend = "C"))

#Get MLE
fit.ou.mle<-find.mle(lik.ou, c(0.1,0.1,0.1))

#Set start parameters (s2, alpha, theta) for MCMC from MLE
pars<-fit.ou.mle$par

#Set prior (uniformative exponential)
prior.exp<-make.prior.exponential(0.1)

#Find values for w (tuning parameter) using short MCMC run
# = 95% confidence interval of parameters from short MCMC run
#Needs abs or theta < 0

tmp <- mcmc(lik.ou, abs(pars), nsteps = 100, w = 0.1, prior = prior.exp)

    w <- diff(sapply(tmp[2:4], quantile, c(.05, .95)))

#Run MCMC
out <- mcmc(lik.ou, abs(pars), nsteps=10000, w=w, prior = prior.exp)

# Look at the marginal distribution of parameters
profiles.plot(out["alpha"], col.line="red") #Error in xy.coords(x, y) : 'x' and 'y' lengths differ
profiles.plot(out["theta"], col.line="blue") #Error in xy.coords(x, y) : 'x' and 'y' lengths differ
profiles.plot(out["s2"], col.line="green") #Runs

[richfitz]

This looks like it was fixed in 5bf5eb4db3597d023c04df15b1c39940fb093cb5 - can you install the current source version and confirm? (I can't replicate this problem with the current github version of diversitree).

[nhcooper123]

Thanks Rich. This is fixed now. Though compiling was a pain!

[richfitz]

:thumbsup: glad it worked.