#prior lambda = 6
#for the geometric mean (body size)
x<-data$GM
names(x)<-rownames(data)
prior <- make.prior(tree, model = "OU", param=list("dalpha"=list(),"dsig2"=list(),"dtheta"=list(),
                                         "dk"=list(lambda=6,kmax=length(tree$tip.label)*2-2),

                                          "dloc"=list(min=0,max=1),
                                         "dsb"=list(bmax=Inf,prob=tree$edge.length)))


fit.1<-bayou.mcmc(tree, x, model = "OU",prior= prior, ngen = 500000)
fit.2<-bayou.mcmc(tree, x,model = "OU", prior= prior, ngen = 500000)
chain.1 <- load.bayou(fit.1, save.Rdata=FALSE, cleanup=TRUE)
chain.2 <- load.bayou(fit.2, save.Rdata=FALSE, cleanup=TRUE)

combine<-combine.chains(chain.1, chain.2, burnin.prop = 0.2)

plotSimmap.mcmc(tree,combine, burnin = .1, fsize = 0.5,circle.lwd = 0.5 )

> effectiveSize(combine$lnL)
    var1
47696.89
> effectiveSize(combine$alpha)
    var1
17668.65
> effectiveSize(combine$k)
    var1
200.3831

http://andrewgelman.com/2010/10/29/could_someone_p/

#GM
#lambda=15
prior <- make.prior(tree, dists=list(dalpha="dlnorm", dsig2="dlnorm",dsb="dsb", dk="cdpois",
                                     dtheta="dnorm"), param=list(dalpha=list(meanlog=-5, sdlog=2),
                                     dsig2=list(meanlog=-1, sdlog=5), dk=list(lambda=15, kmax=200),
                                     dsb=list(bmax=Inf,prob=tree$edge.length),
                                     dtheta=list(mean=mean(x), sd=2)))

fit.1<-bayou.mcmc(tree, x, model = "OU",prior= prior, ngen = 500000)
fit.2<-bayou.mcmc(tree, x,model = "OU", prior= prior, ngen = 500000)
chain.1 <- load.bayou(fit.1, save.Rdata=FALSE, cleanup=TRUE)
chain.2 <- load.bayou(fit.2, save.Rdata=FALSE, cleanup=TRUE)

combine<-combine.chains(chain.1, chain.2, burnin.prop = 0.2)

> effectiveSize(combine$lnL)
    var1
54733.55
> effectiveSize(combine$alpha)
    var1
21377.18
> effectiveSize(combine$k)
    var1
460.0945
> median(combine$alpha)
[1] 0.005785393
> median(combine$k)
[1] 1
> median(combine$lnL)
[1] -226.6656
> median(combine$sig2)
[1] 0.1995975
> median(combine$ntheta)
[1] 2
> median(combine$theta)

plotSimmap.mcmc(tree,combine, burnin = 0, fsize = 0.5,circle.lwd = 0.5 )

#pPC1, lambda = 6
x<-data$PC1
names(x)<-rownames(data)
prior <- make.prior(tree, model = "OU", param=list("dalpha"=list(),"dsig2"=list(),"dtheta"=list(),
                                         "dk"=list(lambda=6,kmax=length(tree$tip.label)*2-2),
                                          "dloc"=list(min=0,max=1),
                                         "dsb"=list(bmax=Inf,prob=tree$edge.length)))
> effectiveSize(combine$lnL)
    var1
77556.68
> effectiveSize(combine$alpha)
    var1
750.1128
> effectiveSize(combine$k)
    var1
1098.867

> median(combine$alpha)

[1] 0.01054099

> median(combine$k)

[1] 5

> median(combine$lnL)

[1] -232.5772

> median(combine$sig2)

[1] 0.2309304

> median(combine$ntheta)

[1] 6

#pPC1

#lambda=15
prior <- make.prior(tree, dists=list(dalpha="dlnorm", dsig2="dlnorm",dsb="dsb", dk="cdpois",
                                     dtheta="dnorm"), param=list(dalpha=list(meanlog=-5, sdlog=2),
                                     dsig2=list(meanlog=-1, sdlog=5), dk=list(lambda=15, kmax=200),
                                     dsb=list(bmax=Inf,prob=tree$edge.length),
                                     dtheta=list(mean=mean(x), sd=2)))

fit.1<-bayou.mcmc(tree, x, model = "OU",prior= prior, ngen = 500000, plot.freq = NULL)
fit.2<-bayou.mcmc(tree, x,model = "OU", prior= prior, ngen = 500000, plot.freq = NULL)
chain.1 <- load.bayou(fit.1, save.Rdata=FALSE, cleanup=TRUE)
chain.2 <- load.bayou(fit.2, save.Rdata=FALSE, cleanup=TRUE)

combine<-combine.chains(chain.1, chain.2, burnin.prop = 0.2)

> effectiveSize(combine$lnL)
    var1
46687.59
> effectiveSize(combine$alpha)
    var1
1882.031
> effectiveSize(combine$k)
    var1
1916.088
> median(combine$alpha)
[1] 0.01097906
> median(combine$k)
[1] 6
> median(combine$lnL)
[1] -682.1394
> median(combine$sig2)
[1] 187.9879
> median(combine$ntheta)
[1] 7