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.

Import packages and data

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)

Data prep

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

Check intercorrelations

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])

Standardize Variables

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

binaryPGLMM

first we will run each of the models detailed in the manuscript using the binaryPGLMM() function in the ape package.

PC1 and PC2

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
Slope
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
TWI
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
Forest cover
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
Humidity (Q)
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
Temperature (BIO8)
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
Precipitation
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

MCMCglmm

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

PC1 and PC2
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
Slope
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
TWI
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
Forest cover
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
Humidity (Q)
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
Temperature (BIO8)
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
Precipitation (BIO16)
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

Plot response curves

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")
TWI
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")

Make plots

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)