This code performs two different versions of phylogenetic GLMMs to test the effect of environmental variables on reproductive mode. Reproductive mode, a discrete variable, is re-coded into a binary, terrestrial/aquatic variable and logistic link regressions can therefore be performed.
packages
library(ape)
library(phytools)
library(geiger)
library(gclus)
library(MCMCglmm)
library(coda)
data
tree<-read.nexus("~/r_scripts/data_files/beast_MCC.tre")
my.data<-read.csv("~/r_scripts/data_files/my_data.csv", header=T, row.names = 1)
match tree and data
tree<-drop.tip(tree, name.check(tree, my.data)$tree_not_data)
my.data<-my.data[tree$tip.label,]
all(rownames(my.data)==tree$tip)
## [1] TRUE
name.check(tree, my.data)
## [1] "OK"
plot(tree, cex=0.5)
trim the dataset down to only the variables we are interested in.
var.of.intr<-c("forest","slope","Q","twi","bio8","bio16", "PC1", "PC2")
dat<-my.data[,var.of.intr]
Re-code reproductive modes to be either terrestrial or not
# define terrestrial species
dat$terrestrial<-rep(0, nrow(my.data))
dat$terrestrial[my.data$rep_mode=="viviparity" | my.data$rep_mode=="terr_eggs"]<-1
We can make a function that produces pairwise plots with colour coding (yellow to pink) to depict the strength of correlation of variables
#make plotting function
corr.plot<-function(mod){
my.abs = abs(cor(mod))
my.colors = dmat.color(my.abs)
my.ordered = order.single(cor(mod))
cpairs(mod,my.ordered,panel.colors=my.colors,gap=0.3, pch=20)
}
#plot variables of interest
corr.plot(dat[,var.of.intr])
The GLMMs will be run on standardized predictors with means of 0 and variances of 1
dat.scaled<-data.frame(cbind(terrestrial=dat$terrestrial, scale(dat[,var.of.intr])))
round(colMeans(dat.scaled[,var.of.intr]),3) # means should be close to 0
## forest slope Q twi bio8 bio16 PC1 PC2
## 0 0 0 0 0 0 0 0
apply(dat.scaled[,var.of.intr], 2, var) # vriances should be close to 1
## forest slope Q twi bio8 bio16 PC1 PC2
## 1 1 1 1 1 1 1 1
first we will run each of the models detailed in the manuscript using the binaryPGLMM() function in the ape package.
m1.pca<-binaryPGLMM(data=dat.scaled, terrestrial~PC1+PC2, tree)
m1.pca$convergeflag
## [1] "converged"
m1.pca$iteration
## [1] 73
m1.pca
##
##
## Call:terrestrial ~ PC1 + PC2
##
## Random effect (phylogenetic signal s2):
## s2 Pr
## 1 5.823 1.299e-05
##
## Fixed effects:
## Value Std.Error Zscore Pvalue
## (Intercept) -3.90516 1.65957 -2.3531 0.01862 *
## PC1 2.66483 1.38066 1.9301 0.05359 .
## PC2 0.46536 0.72663 0.6404 0.52189
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1.slope<-binaryPGLMM(data=dat.scaled, terrestrial~slope, tree)
m1.slope$convergeflag
## [1] "converged"
m1.slope$iteration
## [1] 43
m1.slope
##
##
## Call:terrestrial ~ slope
##
## Random effect (phylogenetic signal s2):
## s2 Pr
## 1 4.606 5.316e-05
##
## Fixed effects:
## Value Std.Error Zscore Pvalue
## (Intercept) -2.65676 1.03062 -2.5778 0.009942 **
## slope 1.43853 0.72375 1.9876 0.046856 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1.twi<-binaryPGLMM(data=dat.scaled, terrestrial~twi, tree, tol.pql = 0.001)
m1.twi$convergeflag
## [1] "converged"
m1.twi$iteration
## [1] 13
m1.twi
##
##
## Call:terrestrial ~ twi
##
## Random effect (phylogenetic signal s2):
## s2 Pr
## 1 4.81 6.381e-05
##
## Fixed effects:
## Value Std.Error Zscore Pvalue
## (Intercept) -3.5910 1.4774 -2.4306 0.01507 *
## twi -2.4430 1.2669 -1.9282 0.05383 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1.forest<-binaryPGLMM(data=dat.scaled, terrestrial~forest, tree)
m1.forest$convergeflag
## [1] "converged"
m1.forest$iteration
## [1] 44
m1.forest
##
##
## Call:terrestrial ~ forest
##
## Random effect (phylogenetic signal s2):
## s2 Pr
## 1 4.872 4.847e-05
##
## Fixed effects:
## Value Std.Error Zscore Pvalue
## (Intercept) -2.35610 1.00993 -2.3329 0.01965 *
## forest 1.28756 0.75988 1.6944 0.09018 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1.Q<-binaryPGLMM(data=dat.scaled, terrestrial~Q, tree)
m1.Q$convergeflag
## [1] "converged"
m1.Q$iteration
## [1] 46
m1.Q
##
##
## Call:terrestrial ~ Q
##
## Random effect (phylogenetic signal s2):
## s2 Pr
## 1 5.437 1.346e-07
##
## Fixed effects:
## Value Std.Error Zscore Pvalue
## (Intercept) -2.41449 1.05471 -2.2892 0.02207 *
## Q 0.91772 0.74243 1.2361 0.21642
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1.bio8<-binaryPGLMM(data=dat.scaled, terrestrial~bio8, tree)
m1.bio8$convergeflag
## [1] "converged"
m1.bio8$iteration
## [1] 41
m1.bio8
##
##
## Call:terrestrial ~ bio8
##
## Random effect (phylogenetic signal s2):
## s2 Pr
## 1 5.152 4.136e-07
##
## Fixed effects:
## Value Std.Error Zscore Pvalue
## (Intercept) -2.28930 1.00023 -2.2888 0.02209 *
## bio8 -0.59522 0.61735 -0.9642 0.33496
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1.bio16<-binaryPGLMM(data=dat.scaled, terrestrial~bio16, tree)
m1.bio16$convergeflag
## [1] "converged"
m1.bio16$iteration
## [1] 38
m1.bio16
##
##
## Call:terrestrial ~ bio16
##
## Random effect (phylogenetic signal s2):
## s2 Pr
## 1 5.047 4.815e-07
##
## Fixed effects:
## Value Std.Error Zscore Pvalue
## (Intercept) -2.18864 0.97019 -2.2559 0.02408 *
## bio16 0.59145 0.66643 0.8875 0.37482
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
second we will run each of the above models again, but this time using the MCMCglmm() function in the MCMCglmm package.
for this, the binary response variable has to be converted to be a factor:
dat.scaled$terrestrial<-as.factor(dat.scaled$terrestrial)
and the priors set up including the phylogenetic error structure:
V <- vcv(tree)
V <- V/max(V)
detV <- exp(determinant(V)$modulus[1])
V <- V/detV^(1/length(tree$tip.label))
invV <- Matrix(solve(V),sparse=T)
dat.scaled$species <- rownames(dat.scaled)
rownames(invV) <- dat.scaled$species
prior <- list(R=list(V=1, fix=1), G=list(G1=list(V=1, nu=1000, alpha.mu=0, alpha.V=1)))
chain parameters are:
##the used parameters were as follows:
#nitt <- 5010000
#thin <- 500
#burnin <- 10000
#(nitt-burnin)/thin
## ideally (nitt-burnin)/thin = 20000
##but for demonstration purposes, these are drastically reduced to:
nitt <- 1001000
thin <- 100
burnin <- 1000
(nitt-burnin)/thin
## [1] 10000
and each model is run twice to check convergence of the chains
m2a.pca<-MCMCglmm(terrestrial~PC1+PC2, random=~species, ginvers=list(species=invV),data=dat.scaled, slice=TRUE, nitt=nitt, thin=thin, burnin=burnin,family="categorical", prior=prior, verbose=FALSE)
m2b.pca<-MCMCglmm(terrestrial~PC1+PC2, random=~species, ginvers=list(species=invV),data=dat.scaled, slice=TRUE, nitt=nitt, thin=thin, burnin=burnin,family="categorical", prior=prior, verbose=FALSE)
chain diagnostics
#check chains:
plot(m2a.pca$Sol)
plot(as.matrix(m2a.pca$Sol[,1]), type="l")
lines(as.matrix(m2b.pca$Sol[,1]), col="red")
## check ESS
as.matrix(effectiveSize(m2a.pca$Sol)) # ESS should be above 200
## [,1]
## (Intercept) 387.1374
## PC1 397.7251
## PC2 3323.8071
## check whether sampling interval (thin) is high enough to avoid autocorrelation
autocorr.plot(m2a.pca$Sol, lag.max=50, auto.layout=T) ### this should show a rappid tapering off to 0
model summary
#check chains:
summary(m2a.pca)
##
## Iterations = 1001:1000901
## Thinning interval = 100
## Sample size = 10000
##
## DIC: 17.41539
##
## G-structure: ~species
##
## post.mean l-95% CI u-95% CI eff.samp
## species 7.673 1.174 16.36 558.6
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 1 1 1 0
##
## Location effects: terrestrial ~ PC1 + PC2
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -9.2191 -16.4742 -2.8418 387.1 0.0002 ***
## PC1 6.2359 1.2528 11.7979 397.7 0.0010 ***
## PC2 0.9275 -1.8468 3.6310 3323.8 0.4658
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m2a.slope<-MCMCglmm(terrestrial~slope, random=~species, ginvers=list(species=invV),data=dat.scaled, slice=TRUE, nitt=nitt, thin=thin, burnin=burnin,family="categorical", prior=prior, verbose=FALSE)
m2b.slope<-MCMCglmm(terrestrial~slope, random=~species, ginvers=list(species=invV),data=dat.scaled, slice=TRUE, nitt=nitt, thin=thin, burnin=burnin,family="categorical", prior=prior, verbose=FALSE)
chain diagnostics
#check chains:
plot(m2a.slope$Sol)
plot(as.matrix(m2a.slope$Sol[,1]), type="l")
lines(as.matrix(m2b.slope$Sol[,1]), col="red")
## check ESS
as.matrix(effectiveSize(m2a.slope$Sol)) # ESS should be above 200
## [,1]
## (Intercept) 1262.202
## slope 1417.944
## check whether sampling interval (thin) is high enough to avoid autocorrelation
autocorr.plot(m2a.slope$Sol, lag.max=50, auto.layout=T) ### this should show a rappid tapering off to 0
model summary
#check chains:
summary(m2a.slope)
##
## Iterations = 1001:1000901
## Thinning interval = 100
## Sample size = 10000
##
## DIC: 20.88299
##
## G-structure: ~species
##
## post.mean l-95% CI u-95% CI eff.samp
## species 6.563 0.8234 14.58 1109
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 1 1 1 0
##
## Location effects: terrestrial ~ slope
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -5.8417 -10.7523 -1.5306 1262 0.0004 ***
## slope 3.0653 0.4462 6.1785 1418 0.0112 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m2a.twi<-MCMCglmm(terrestrial~twi, random=~species, ginvers=list(species=invV),data=dat.scaled, slice=TRUE, nitt=nitt, thin=thin, burnin=burnin,family="categorical", prior=prior, verbose=FALSE)
m2b.twi<-MCMCglmm(terrestrial~twi, random=~species, ginvers=list(species=invV),data=dat.scaled, slice=TRUE, nitt=nitt, thin=thin, burnin=burnin,family="categorical", prior=prior, verbose=FALSE)
chain diagnostics
#check chains:
plot(m2a.twi$Sol)
plot(as.matrix(m2a.twi$Sol[,1]), type="l")
lines(as.matrix(m2b.twi$Sol[,1]), col="red")
## check ESS
as.matrix(effectiveSize(m2a.twi$Sol)) # ESS should be above 200
## [,1]
## (Intercept) 271.6655
## twi 263.1794
## check whether sampling interval (thin) is high enough to avoid autocorrelation
autocorr.plot(m2a.twi$Sol, lag.max=50, auto.layout=T) ### this should show a rappid tapering off to 0
model summary
#check chains:
summary(m2a.twi)
##
## Iterations = 1001:1000901
## Thinning interval = 100
## Sample size = 10000
##
## DIC: 18.97408
##
## G-structure: ~species
##
## post.mean l-95% CI u-95% CI eff.samp
## species 6.79 0.7207 15.03 510.9
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 1 1 1 0
##
## Location effects: terrestrial ~ twi
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -8.410 -16.155 -2.377 271.7 2e-04 ***
## twi -5.728 -11.833 -0.830 263.2 1e-03 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m2a.forest<-MCMCglmm(terrestrial~forest, random=~species, ginvers=list(species=invV),data=dat.scaled, slice=TRUE, nitt=nitt, thin=thin, burnin=burnin,family="categorical", prior=prior, verbose=FALSE)
m2b.forest<-MCMCglmm(terrestrial~forest, random=~species, ginvers=list(species=invV),data=dat.scaled, slice=TRUE, nitt=nitt, thin=thin, burnin=burnin,family="categorical", prior=prior, verbose=FALSE)
chain diagnostics
#check chains:
plot(m2a.forest$Sol)
plot(as.matrix(m2a.forest$Sol[,1]), type="l")
lines(as.matrix(m2b.forest$Sol[,1]), col="red")
## check ESS
as.matrix(effectiveSize(m2a.forest$Sol)) # ESS should be above 200
## [,1]
## (Intercept) 1873.534
## forest 1990.068
## check whether sampling interval (thin) is high enough to avoid autocorrelation
autocorr.plot(m2a.forest$Sol, lag.max=50, auto.layout=T) ### this should show a rappid tapering off to 0
model summary
#check chains:
summary(m2a.forest)
##
## Iterations = 1001:1000901
## Thinning interval = 100
## Sample size = 10000
##
## DIC: 21.36982
##
## G-structure: ~species
##
## post.mean l-95% CI u-95% CI eff.samp
## species 6.679 0.8402 14.36 1319
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 1 1 1 0
##
## Location effects: terrestrial ~ forest
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -5.17902 -9.53072 -0.99711 1874 0.0052 **
## forest 2.73549 0.01994 5.85252 1990 0.0368 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m2a.Q<-MCMCglmm(terrestrial~Q, random=~species, ginvers=list(species=invV),data=dat.scaled, slice=TRUE, nitt=nitt, thin=thin, burnin=burnin,family="categorical", prior=prior, verbose=FALSE)
m2b.Q<-MCMCglmm(terrestrial~Q, random=~species, ginvers=list(species=invV),data=dat.scaled, slice=TRUE, nitt=nitt, thin=thin, burnin=burnin,family="categorical", prior=prior, verbose=FALSE)
chain diagnostics
#check chains:
plot(m2a.Q$Sol)
plot(as.matrix(m2a.Q$Sol[,1]), type="l")
lines(as.matrix(m2b.Q$Sol[,1]), col="red")
## check ESS
as.matrix(effectiveSize(m2a.Q$Sol)) # ESS should be above 200
## [,1]
## (Intercept) 1308.937
## Q 1738.889
## check whether sampling interval (thin) is high enough to avoid autocorrelation
autocorr.plot(m2a.Q$Sol, lag.max=50, auto.layout=T) ### this should show a rappid tapering off to 0
model summary
#check chains:
summary(m2a.Q)
##
## Iterations = 1001:1000901
## Thinning interval = 100
## Sample size = 10000
##
## DIC: 20.19528
##
## G-structure: ~species
##
## post.mean l-95% CI u-95% CI eff.samp
## species 8.272 1.331 16.88 1036
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 1 1 1 0
##
## Location effects: terrestrial ~ Q
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -5.8666 -11.2851 -1.5559 1309 0.0032 **
## Q 2.3707 -0.6499 5.5311 1739 0.0880 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m2a.bio8<-MCMCglmm(terrestrial~bio8, random=~species, ginvers=list(species=invV),data=dat.scaled, slice=TRUE, nitt=nitt, thin=thin, burnin=burnin,family="categorical", prior=prior, verbose=FALSE)
m2b.bio8<-MCMCglmm(terrestrial~bio8, random=~species, ginvers=list(species=invV),data=dat.scaled, slice=TRUE, nitt=nitt, thin=thin, burnin=burnin,family="categorical", prior=prior, verbose=FALSE)
chain diagnostics
#check chains:
plot(m2a.bio8$Sol)
plot(as.matrix(m2a.bio8$Sol[,1]), type="l")
lines(as.matrix(m2b.bio8$Sol[,1]), col="red")
## check ESS
as.matrix(effectiveSize(m2a.bio8$Sol)) # ESS should be above 200
## [,1]
## (Intercept) 1506.324
## bio8 1978.726
## check whether sampling interval (thin) is high enough to avoid autocorrelation
autocorr.plot(m2a.bio8$Sol, lag.max=50, auto.layout=T) ### this should show a rappid tapering off to 0
model summary
#check chains:
summary(m2a.bio8)
##
## Iterations = 1001:1000901
## Thinning interval = 100
## Sample size = 10000
##
## DIC: 21.57676
##
## G-structure: ~species
##
## post.mean l-95% CI u-95% CI eff.samp
## species 7.933 1.397 16.58 1095
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 1 1 1 0
##
## Location effects: terrestrial ~ bio8
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -5.4514 -10.4285 -1.2647 1506 0.0034 **
## bio8 -1.4247 -4.0036 0.8596 1979 0.2070
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m2a.bio16<-MCMCglmm(terrestrial~bio16, random=~species, ginvers=list(species=invV),data=dat.scaled, slice=TRUE, nitt=nitt, thin=thin, burnin=burnin,family="categorical", prior=prior, verbose=FALSE)
m2b.bio16<-MCMCglmm(terrestrial~bio16, random=~species, ginvers=list(species=invV),data=dat.scaled, slice=TRUE, nitt=nitt, thin=thin, burnin=burnin,family="categorical", prior=prior, verbose=FALSE)
chain diagnostics
#check chains:
plot(m2a.bio16$Sol)
plot(as.matrix(m2a.bio16$Sol[,1]), type="l")
lines(as.matrix(m2b.bio16$Sol[,1]), col="red")
## check ESS
as.matrix(effectiveSize(m2a.bio16$Sol)) # ESS should be above 200
## [,1]
## (Intercept) 1725.170
## bio16 2911.496
## check whether sampling interval (thin) is high enough to avoid autocorrelation
autocorr.plot(m2a.bio16$Sol, lag.max=50, auto.layout=T) ### this should show a rappid tapering off to 0
model summary
#check chains:
summary(m2a.bio16)
##
## Iterations = 1001:1000901
## Thinning interval = 100
## Sample size = 10000
##
## DIC: 21.98019
##
## G-structure: ~species
##
## post.mean l-95% CI u-95% CI eff.samp
## species 7.924 1.276 17.11 1054
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 1 1 1 0
##
## Location effects: terrestrial ~ bio16
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -5.1641 -9.8260 -0.9323 1725 0.005 **
## bio16 1.3986 -1.1339 4.1557 2911 0.250
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
To plot the predicted curves for the figure output, first re-run the mcmc on the non-scaled data for variables of interest, then predict points ##### Slope
dat$species <- rownames(dat)
m.slope<-MCMCglmm(terrestrial~slope, random=~species, ginvers=list(species=invV),data=dat, slice=TRUE, nitt=nitt, thin=thin, burnin=burnin,family="categorical", prior=prior, verbose=FALSE)
p.slope<-predict.MCMCglmm(m.slope, interval="confidence")
m.twi<-MCMCglmm(terrestrial~twi, random=~species, ginvers=list(species=invV),data=dat, slice=TRUE, nitt=nitt, thin=thin, burnin=burnin,family="categorical", prior=prior, verbose=FALSE)
p.twi<-predict.MCMCglmm(m.twi, interval="confidence")
Response curves
#divide plot area:
par(mfrow=c(1,2))
par(mar=c(2,4,1,1))
#plot predicted response curves for slope:
plot(NA,xaxt="n", yaxt="n", xlim=c(0,12), ylim=c(0,1), bty="n", xlab=NA, ylab="Probability")
axis(2, at=c(0,0.5,1), las=2)
axis(1, las=1)
polygon(x=c(sort(dat$slope), sort(dat$slope,decreasing=T)), y=c(sort(as.numeric(p.slope[,"lwr"])), sort(as.numeric(p.slope[,"upr"]), decreasing = T)), border = NA, col="grey95")
points(x=dat$slope, y=dat$terrestrial, pch=20, col=ifelse(dat$terrestrial==1, "black", "grey"))
lines(x=sort(dat$slope), y=sort(as.numeric(p.slope[,"fit"])))
#plot predicted response curves for TWI:
plot(NA,xaxt="n", yaxt="n", xlim=c(7,16), ylim=c(0,1), bty="n", xlab=NA, ylab="Probability")
axis(2, at=c(0,0.5,1), las=2)
axis(1, las=1)
polygon(x=c(sort(dat$twi), sort(dat$twi,decreasing=T)), y=c(sort(as.numeric(p.twi[,"lwr"]), decreasing = T), sort(as.numeric(p.twi[,"upr"]), decreasing = F)), border = NA, col="grey95")
points(x=dat$twi, y=dat$terrestrial, pch=20, col=ifelse(dat$terrestrial==1, "black", "grey"))
lines(x=sort(dat$twi), y=sort(as.numeric(p.twi[,"fit"]), decreasing = T))
Phylograms
#divide plot area:
par(mfrow=c(1,2))
par(mar=c(4,4,0,1))
### make a probability matrix of tip states
prob.mat<-data.frame(model.matrix(~as.factor(dat$terrestrial)-1))
colnames(prob.mat)<-c("0","1")
rownames(prob.mat)<-tree$tip
prob.mat<-prob.mat[tree$tip,]
prob.mat<-as.matrix(prob.mat)
mtree<-make.simmap(tree, x=prob.mat, model="ER", nsim=1)
## make.simmap is sampling character histories conditioned on the transition matrix
##
## Q =
## 0 1
## 0 -0.006131104 0.006131104
## 1 0.006131104 -0.006131104
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## 0 1
## 0.5 0.5
cols<-setNames(c("black","grey"),c(1,0))
slope.dat<-dat$slope; names(slope.dat)<-rownames(dat)
twi.dat<-dat$twi; names(twi.dat)<-rownames(dat)
phenogram(mtree,slope.dat, ylim=c(0, 12), fsize=0.0000001, spread.labels=F, colors = cols)
phenogram(mtree,twi.dat, ylim=c(7,16), fsize=0.0000001, spread.labels=F, colors = cols)