Phylogeny and Species Trait Effects on Detectability

library(lhreg)
data(lhreg_data)
str(lhreg_data)

with(lhreg_data, plot(exp(logphi), exp(logtau),
    cex=logmass*0.5, col=Mig2, pch=c(21, 22)[Hab2]))
legend("topright", bty="n", pch=c(21, 21, 22, 22), col=c(1,2,1,2),
    legend=c("Migratory/Closed", "Resident/Closed",
    "Migratory/Open", "Resident/Open"))

library(ape)
mph <- read.nexus("11960.tre") # 1000 trees with Ericson backbone
lhreg_tree <- consensus(mph)
table(sapply(mph, function(z) length(z$tip.label)))
CORR <- TRUE
vv <- list()
vv[[1]] <- vcv(mph[[1]], corr=CORR)
for (i in 2:length(mph)) {
    v <- vcv(mph[[i]], corr=CORR)
    v <- v[rownames(vv[[1]]), colnames(vv[[1]])]
    vv[[i]] <- v
}
vvv <- v
for (i in 1:length(v)) {
    vvv[i] <- mean(sapply(vv, function(z) z[i]))
}
spp <- intersect(rownames(lhreg_data), rownames(vvv))
vvv <- vvv[spp,spp]
cor_matrix <- as.matrix(nearPD(vvv, corr=TRUE)$mat)

library(ape)
data(cor_matrix)
data(lhreg_tree)

str(cor_matrix)

heatmap(cor_matrix)

library(lhreg)
data(lhreg_data)
data(cor_matrix)

y <- lhreg_data$logphi
m4 <- lm(y ~ logmass + Mig2 + Nesthm + xMaxFreqkHz + Hab4, lhreg_data)
m3 <- lm(y ~ logmass + Mig2 + Nesthm + xMaxFreqkHz + Hab3, lhreg_data)
m2 <- lm(y ~ logmass + Mig2 + Nesthm + xMaxFreqkHz + Hab2, lhreg_data)
m4s <- step(m4, trace=0)
m3s <- step(m3, trace=0)
m2s <- step(m2, trace=0)
AIC(m4,m3,m2,m4s,m3s,m2s)

summary(m2s)

mp <- update(m2s, .~.-Nesthm)

amp <- anova(mp)
round(structure(100 * amp[,"Sum Sq"] / sum(amp[,"Sum Sq"]),
    names=rownames(amp)), 1)

mp0 <- lm(y ~ 1)
summary(mp0) # null model for SR

summary(mp) # full model for SR

## evaluating interactions
mpx <- lm(y ~ logmass + Mig2 + Nesthm + xMaxFreqkHz + Hab2 +
    logmass:xMaxFreqkHz + Hab2:xMaxFreqkHz, lhreg_data)
mpx2 <- step(mpx)

summary(mpx2)

y <- lhreg_data$logtau
m4 <- lm(y ~ logmass + Mig2 + Nesthm + xMaxFreqkHz + Hab4, lhreg_data)
m3 <- lm(y ~ logmass + Mig2 + Nesthm + xMaxFreqkHz + Hab3, lhreg_data)
m2 <- lm(y ~ logmass + Mig2 + Nesthm + xMaxFreqkHz + Hab2, lhreg_data)
m4s <- step(m4, trace=0)
m3s <- step(m3, trace=0)
m2s <- step(m2, trace=0)
AIC(m4,m3,m2,m4s,m3s,m2s)

summary(m2s)

mt <- update(m2s, .~.-Nesthm)

amt <- anova(mt)
round(structure(100 * amt[,"Sum Sq"] / sum(amt[,"Sum Sq"]),
    names=rownames(amt)), 1)

mt0 <- lm(y ~ 1)
summary(mp0) # null model for DD

summary(mp) # full model for DD

## evaluating interactions
mtx <- lm(y ~ logmass + Mig2 + Nesthm + xMaxFreqkHz + Hab2 +
    logmass:xMaxFreqkHz + Hab2:xMaxFreqkHz, lhreg_data)
mtx2 <- step(mtx)

summary(mtx2)

met <- "DE" # can use "SANN", "Nelder-Mead" is quickest
x <- lhreg_data
vc <- cor_matrix

## model matrix definitions
X0 <- matrix(1, nrow(x), 1) # null (for both)
colnames(X0) <- "Intercept"
Xt <- model.matrix(mt) # DD
Xp <- model.matrix(mp) # SR

## fit models for tau

tM00 <- lhreg(Y=x$logtau, X=X0, SE=x$logtauSE, V=vc, lambda=0,
    hessian=TRUE, method=met)
tMl0 <- lhreg(Y=x$logtau, X=X0, SE=x$logtauSE, V=vc, lambda=NA,
    hessian=TRUE, method=met)
tM0b <- lhreg(Y=x$logtau, X=Xt, SE=x$logtauSE, V=vc, lambda=0,
    hessian=TRUE, method=met)
tMlb <- lhreg(Y=x$logtau, X=Xt, SE=x$logtauSE, V=vc, lambda=NA,
    hessian=TRUE, method=met)

## fit models for phi

pM00 <- lhreg(Y=x$logphi, X=X0, SE=x$logphiSE, V=vc, lambda=0,
    hessian=TRUE, method=met)
pMl0 <- lhreg(Y=x$logphi, X=X0, SE=x$logphiSE, V=vc, lambda=NA,
    hessian=TRUE, method=met)
pM0b <- lhreg(Y=x$logphi, X=Xp, SE=x$logphiSE, V=vc, lambda=0,
    hessian=TRUE, method=met)
pMlb <- lhreg(Y=x$logphi, X=Xp, SE=x$logphiSE, V=vc, lambda=NA,
    hessian=TRUE, method=met)

## set up lambda values to evaluate at
lam <- seq(0, 1, by=0.01)

## parallel computation is faster
cl <- makeCluster(4)
## load package on workers
tmp <- clusterEvalQ(cl, library(lhreg))

object <- tMl0
clusterExport(cl, "object")
pl_tMl0 <- pbsapply(lam, function(z, ...) profile_lambda1(object, z, ...),
    cl=cl, method=met)

object <- tMlb
clusterExport(cl, "object")
pl_tMlb <- pbsapply(lam, function(z, ...) profile_lambda1(object, z, ...),
    cl=cl, method=met)

object <- pMl0
clusterExport(cl, "object")
pl_pMl0 <- pbsapply(lam, function(z, ...) profile_lambda1(object, z, ...),
    cl=cl, method=met)

object <- pMlb
clusterExport(cl, "object")
pl_pMlb <- pbsapply(lam, function(z, ...) profile_lambda1(object, z, ...),
    cl=cl, method=met)

stopCluster(cl) # close cluster

n <- nrow(x) # we will do n runs

cl <- makeCluster(4) # parallel if you wish
tmp <- clusterEvalQ(cl, library(lhreg)) # load package

loo_tM00 <- t(pbsapply(1:n, loo1, object=tM00, cl=cl))
#loo_tM00 <- t(pbsapply(1:n, loo2, object=tM00, cl=cl, method=met))
loo_tMl0 <- t(pbsapply(1:n, loo2, object=tMl0, cl=cl, method=met))
loo_tM0b <- t(pbsapply(1:n, loo1, object=tM0b, cl=cl))
#loo_tM0b <- t(pbsapply(1:n, loo2, object=tM0b, cl=cl, method=met))
loo_tMlb <- t(pbsapply(1:n, loo2, object=tMlb, cl=cl, method=met))

loo_pM00 <- t(pbsapply(1:n, loo1, object=pM00, cl=cl))
#loo_pM00 <- t(pbsapply(1:n, loo2, object=pM00, cl=cl, method=met))
loo_pMl0 <- t(pbsapply(1:n, loo2, object=pMl0, cl=cl, method=met))
loo_pM0b <- t(pbsapply(1:n, loo1, object=pM0b, cl=cl))
#loo_pM0b <- t(pbsapply(1:n, loo2, object=pM0b, cl=cl, method=met))
loo_pMlb <- t(pbsapply(1:n, loo2, object=pMlb, cl=cl, method=met))

stopCluster(cl)

nsim <- 999
## simulated data is nice, thus quick optimization works
metQ <- "Nelder-Mead"

cl <- makeCluster(4) # parallel if you wish
tmp <- clusterEvalQ(cl, library(lhreg)) # load package

pb_tM00 <- parametric_bootstrap(tM00, nsim=nsim, method=metQ, cl=cl)
pb_tMl0 <- parametric_bootstrap(tMl0, nsim=nsim, method=metQ, cl=cl)
pb_tM0b <- parametric_bootstrap(tM0b, nsim=nsim, method=metQ, cl=cl)
pb_tMlb <- parametric_bootstrap(tMlb, nsim=nsim, method=metQ, cl=cl)

pb_pM00 <- parametric_bootstrap(pM00, nsim=nsim, method=metQ, cl=cl)
pb_pMl0 <- parametric_bootstrap(pMl0, nsim=nsim, method=metQ, cl=cl)
pb_pM0b <- parametric_bootstrap(pM0b, nsim=nsim, method=metQ, cl=cl)
pb_pMlb <- parametric_bootstrap(pMlb, nsim=nsim, method=metQ, cl=cl)

stopCluster(cl)

cl <- makeCluster(4)
pit <- pred_int(tM0b, pb_tM0b, cl=cl)
pip <- pred_int(pMlb, pb_pMlb, cl=cl)
stopCluster(cl)

## save the results:

save(list=c("cor_matrix", "lam", "lhreg_data", "met", "n",
    "vc", "x", "X0", "Xp", "Xt",
    "amp", "mp", "mp0", "mpx", "mpx2",
    "amt", "mt", "mt0", "mtx", "mtx2",
    "pM00", "pM0b", "pMl0", "pMlb",
    "tM00", "tM0b", "tMl0", "tMlb",
    "pl_pMl0", "pl_pMlb", "pl_tMl0", "pl_tMlb",
    "loo_pM00", "loo_pM0b", "loo_pMl0", "loo_pMlb",
    "loo_tM00", "loo_tM0b", "loo_tMl0", "loo_tMlb",
    "pb_pM00", "pb_pM0b", "pb_pMl0", "pb_pMlb",
    "pb_tM00", "pb_tM0b", "pb_tMl0", "pb_tMlb",
    "pit", "pip"),
    file="lhreg-results-DE2.rda")

get_CI <- function(x, level=0.95)
    t(apply(x$estimates, 2, quantile, c((1-level)/2, 1-(1-level)/2)))
sf <- function(z, loo, pb, type=c("wald", "jack", "boot"), level=0.95) {
    type <- match.arg(type)
    zz <- z$summary
    a <- c((1-level)/2, 1-(1-level)/2)
    dig <- 2
    if (type == "wald") {
        pcut <- function(p) {
            factor(c("***", "**", "*", "+", "ns")[as.integer(cut(p,
                c(1, 0.1, 0.05, 0.01, 0.001, 0),
                include.lowest=TRUE, right=FALSE))],
                levels=c("***", "**", "*", "+", "ns"))
        }
        out <- structure(sapply(1:nrow(zz), function(i)
            paste0(round(zz[i,1], dig), " (SE +/- ", round(zz[i,2], dig),
            pcut(zz[i,4]), ")")), names=rownames(zz))
    }
    if (type == "jack") {
        CI <- t(apply(rbind(zz[,1], loo[,3:ncol(loo),drop=FALSE]), 2,
            quantile, a))
        out <- structure(sapply(1:nrow(zz), function(i)
            paste0(round(zz[i,1], dig), " (", round(CI[i,1], dig), ", ",
            round(CI[i,2], dig), ")")), names=rownames(zz))
    }
    if (type == "boot") {
        CI <- get_CI(pb, level)
        out <- structure(sapply(1:nrow(zz), function(i)
            paste0(round(zz[i,1], dig), " (", round(CI[i,1], dig), ", ",
            round(CI[i,2], dig), ")")), names=rownames(zz))
    }
    out
}
## AIC tables
aict <- AIC(tM00, tMl0, tM0b, tMlb)
aict$dAIC <- aict$AIC-min(aict$AIC)

aicp <- AIC(pM00, pMl0, pM0b, pMlb)
aicp$dAIC <- aicp$AIC-min(aicp$AIC)

cbind(aict, t(sapply(list(tM00, tMl0, tM0b, tMlb),
    function(z) z$summary[c("sigma","lambda"), 1])))

cbind(aicp, t(sapply(list(pM00, pMl0, pM0b, pMlb),
    function(z) z$summary[c("sigma","lambda"), 1])))

## MSE
SSEt <- c(
    tM00 = sum((loo_tM00[,"pred"] - tM00$Y)^2),
    tMl0 = sum((loo_tMl0[,"pred"] - tMl0$Y)^2),
    tM0b = sum((loo_tM0b[,"pred"] - tM0b$Y)^2),
    tMlb = sum((loo_tMlb[,"pred"] - tMlb$Y)^2))
SSEp <- c(
    pM00 = sum((loo_pM00[,"pred"] - pM00$Y)^2),
    pMl0 = sum((loo_pMl0[,"pred"] - pMl0$Y)^2),
    pM0b = sum((loo_pM0b[,"pred"] - pM0b$Y)^2),
    pMlb = sum((loo_pMlb[,"pred"] - pMlb$Y)^2))
MSEt <- SSEt / n
MSEp <- SSEp / n

Type <- "boot" # can be wald, jack, boot
zzz <- list(
    M00=sf(pM00, loo_pM00, pb_pM00, Type, Level),
    Ml0=sf(pMl0, loo_pMl0, pb_pMl0, Type, Level),
    M0b=sf(pM0b, loo_pM0b, pb_pM0b, Type, Level),
    Mlb=sf(pMlb, loo_pMlb, pb_pMlb, Type, Level))
m <- matrix("", length(zzz[[4]]), 4)
rownames(m) <- names(zzz[[4]])
colnames(m) <- c("M00", "Ml0", "M0b", "Mlb")
for (i in 1:4) {
    j <- match(names(zzz[[i]]), rownames(m))
    m[j,i] <- zzz[[i]]
}
m["lambda", c("M00", "M0b")] <- "0 (fixed)"
#m[m==""] <- "n/a"
m1 <- m

zzz <- list(
    M00=sf(tM00, loo_tM00, pb_tM00, Type, Level),
    Ml0=sf(tMl0, loo_tMl0, pb_tMl0, Type, Level),
    M0b=sf(tM0b, loo_tM0b, pb_tM0b, Type, Level),
    Mlb=sf(tMlb, loo_tMlb, pb_tMlb, Type, Level))
m <- matrix("", length(zzz[[4]]), 4)
rownames(m) <- names(zzz[[4]])
colnames(m) <- c("M00", "Ml0", "M0b", "Mlb")
for (i in 1:4) {
    j <- match(names(zzz[[i]]), rownames(m))
    m[j,i] <- zzz[[i]]
}
m["lambda", c("M00", "M0b")] <- "0 (fixed)"
#m[m==""] <- "n/a"
m2 <- m

m1 <- rbind(m1, df=aicp$df, dAIC=round(aicp$dAIC, 3),
    XV_MSE=round(MSEp, 3))
m2 <- rbind(m2, df=aict$df, dAIC=round(aict$dAIC, 3),
    XV_MSE=round(MSEt, 3))

print.default(m1, quote=FALSE) # SR results

print.default(m2, quote=FALSE) # DD results

Var_fun <- function(MSE) {
    Var <- 100*(MSE[1]-MSE[-1])/MSE[1]
    Var <- c(Var, Var[1]+Var[2]-Var[3])
    names(Var) <- c("phylo", "traits", "both", "shared")
    Var
}
## SR
round(Var_fun(MSEp), 1)

## DD
round(Var_fun(MSEt), 1)

library(ape)
load(system.file("extdata", "mph.rda", package = "lhreg"))
tre <- mph[[1000]] # pick one tree

ii <- match(tre$tip.label, rownames(lhreg_data))
d2 <- lhreg_data[ii,]
NAMES <- paste0(d2$common_name, " (", d2$spp, ")")
tre$tip.label <- NAMES
xy <- data.frame(
    x=node.depth.edgelength(tre)[1:length(tre$tip.label)],
    y=node.height(tre)[1:length(tre$tip.label)])
phy_pts_size <- function(z, vari, ...)
    points(xy[,1] + z, xy[,2], pch=19, cex=size_fun(vari, Min=0.3, Max=1.2), ...)
phy_pts_col <- function(z, vari, col=1, ...)
    points(xy[,1] + z, xy[,2], pch=19,
        cex=size_fun(vari, Min=0.3, Max=1.2),
        col=size_col_fun(vari, colorRampPalette(c("#000000", col)), Min=0.3), ...)
phy_pts_2 <- function(z, vari, cex=0.6, col="#000000", ...) {
    points(xy[,1] + z, xy[,2], pch=19,
        col=c(col, "#FFFFFF")[as.integer(vari)], cex=cex, ...)
    points(xy[,1] + z, xy[,2], pch=21, col=col, cex=cex)
}

#pdf("Fig1.pdf", width=10, height=15)
op <- par(mar=c(1,1,1,1))

plot(tre, cex=0.6, label.offset=22, font=1)
segments(xy[,1], xy[,2], xy[,1]+21, xy[,2], lwd=1, col=1)
phy_pts_size(2, exp(d2$logphi), col=Col1)
phy_pts_size(5, exp(d2$logtau), col=Col3)
phy_pts_2(10, d2$Mig2, col=1)
phy_pts_size(13, d2$MaxFreqkHz, col="#4daf4a")
phy_pts_size(16, sqrt(exp(d2$logmass)), col="#984ea3")
phy_pts_2(19, d2$Hab2, col="#ff7f00")
text(xy[1,1]+c(2,5,10,13,16,19), rep(142, 6),
    c("SR", "DD", "Migr", "Pitch", "Mass", "Habitat"),
    pos=4, srt=90, cex=0.65, offset=0)

par(op)
#dev.off()

#pdf("FigPL.pdf", width=14, height=6)
op <- par(mfrow=c(1,2))

Res1 <- pl_pMl0
Res2 <- pl_pMlb
yv <- exp(Res1-max(Res1))
plot(lam, yv, type="n", lwd=1, ylim=exp(c(2*Crit, 0)),
    xlab=expression(lambda), ylab="Profile Likelihood Ratio")
tmp1 <- splinefun(lam, yv)(ltmp)
i1 <- tmp1 > exp(Crit)
yv <- exp(Res2-max(Res2))
tmp2 <- splinefun(lam, yv)(ltmp)
i2 <- tmp2 > exp(Crit)
polygon(c(ltmp[i1], rev(ltmp[i1])),
    c(tmp1[i1], rep(-1, sum(i1))),
    col=Col2, border=NA)
polygon(c(ltmp[i2], rev(ltmp[i2])),
    c(tmp2[i2], rep(-1, sum(i2))),
    col=Col4, border=NA)
abline(h=exp(Crit), col="grey")
lines(ltmp, tmp1, col=Col1, lwd=1)
lines(ltmp[i1], tmp1[i1], col=Col1, lwd=3)
lines(rep(ltmp[which.max(tmp1)], 2), c(1, -1), col=Col1)
lines(ltmp, tmp2, col=Col3, lwd=1)
lines(ltmp[i2], tmp2[i2], col=Col3, lwd=3)
lines(rep(ltmp[which.max(tmp2)], 2), c(1, -1), col=Col3)
box()
legend(0.1, 1, bty="n", border=c(Col1, Col3), fill=c(Col2, Col4), pt.cex=2, pt.lwd=2,
    legend=c(expression(M[lambda*0]-SR), expression(M[lambda*beta]-SR)))
round(c(Max_Ml0=ltmp[which.max(tmp1)], CI=range(ltmp[i1])), 3)

round(c(Max_Mlx=ltmp[which.max(tmp2)], CI=range(ltmp[i2])), 3)

Res1 <- pl_tMl0
Res2 <- pl_tMlb
yv <- exp(Res1-max(Res1))
plot(lam, yv, type="n", lwd=1, ylim=exp(c(2*Crit, 0)),
    xlab=expression(lambda), ylab="Profile Likelihood Ratio")
tmp1 <- splinefun(lam, yv)(ltmp)
i1 <- tmp1 > exp(Crit)
yv <- exp(Res2-max(Res2))
tmp2 <- splinefun(lam, yv)(ltmp)
i2 <- tmp2 > exp(Crit)
polygon(c(ltmp[i1], rev(ltmp[i1])),
    c(tmp1[i1], rep(-1, sum(i1))),
    col=Col2, border=NA)
polygon(c(ltmp[i2], rev(ltmp[i2])),
    c(tmp2[i2], rep(-1, sum(i2))),
    col=Col4, border=NA)
abline(h=exp(Crit), col="grey")
lines(ltmp, tmp1, col=Col1, lwd=1)
lines(ltmp[i1], tmp1[i1], col=Col1, lwd=3)
lines(rep(ltmp[which.max(tmp1)], 2), c(1, -1), col=Col1)
lines(ltmp, tmp2, col=Col3, lwd=1)
lines(ltmp[i2], tmp2[i2], col=Col3, lwd=3)
lines(rep(ltmp[which.max(tmp2)], 2), c(1, -1), col=Col3)
box()
legend(0.1, 1, bty="n", border=c(Col1, Col3), fill=c(Col2, Col4), pt.cex=2, pt.lwd=2,
    legend=c(expression(M[lambda*0]-DD), expression(M[lambda*beta]-DD)))

round(c(Max_Ml0=ltmp[which.max(tmp1)], CI=range(ltmp[i1])), 3)

round(c(Max_Mlx=ltmp[which.max(tmp2)], CI=range(ltmp[i2])), 3)

par(op)
#dev.off()

#pdf("FigPB.pdf", width=10, height=10)
op <- par(mfrow=c(2,2))
plot(density(pb_pMl0$estimates[,"lambda"]),
    xlim=c(0,1), main="SR Ml0", col=Col1)
plot(density(pb_pMlb$estimates[,"lambda"]),
    xlim=c(0,1), main="SR Mlb", col=Col1)
plot(density(pb_tMl0$estimates[,"lambda"]),
    xlim=c(0,1), main="DD Ml0", col=Col3)
plot(density(pb_tMlb$estimates[,"lambda"]),
    xlim=c(0,1), main="DD Mlb", col=Col3)

par(op)
#dev.off()

#pdf("Fig2.pdf", width=12.75, height=7)
op <- par(mfrow=c(1,2))

Max <- 0.7
D <- as.matrix(dist(prp))
diag(D) <- Inf
D <- apply(D, 1, min)
plot(prp, xaxs = "i", yaxs = "i", type="n", asp=1,
    ylim=c(0, Max), xlim=c(0, Max),
    xlab="Time-removal SR Estimate (/min)",
    ylab="LOO SR Estimate (/min)")
abline(0,2,lty=2, col="grey")
abline(0,1/2,lty=2, col="grey")
abline(0,1.5,lty=2, col="grey")
abline(0,1/1.5,lty=2, col="grey")
abline(0,1,lty=1, col=1)
r0 <- size_fun(x$MaxFreqkHz, 0.2*Max/50, Max/50)
draw.ellipse(prp[,1], prp[,2], r0, r0,
    border=c(Col1,Col3)[as.integer(x$Mig2)])
segments(prp[,1], PIp[,1], prp[,1], PIp[,2],
    col=c(Col1,Col3)[as.integer(x$Mig2)], lwd=1)
box()
legend("topleft", pch=21, col=c(Col1,Col3),
    pt.cex=c(2,2), bty="n",
    legend=c("Migrant", "Resident"))
r <- seq(0.2*Max/50, Max/50, len=5)
draw.ellipse(rep(0.06, 5), 0.55+r-max(r), r, r, border=1)
text(c(0.09, 0.06, 0.06), c(0.61, 0.52, 0.58),
    c("Song Pitch (kHz)",round(range(x$MaxFreqkHz),1)),
    cex=c(1,0.75,0.75))
text(0.9*Max, 0.05*Max, expression(M[lambda*beta]-SR))
#Ti <- D > 0.03
#Ti <- D > 0.03 | (prp[,1]/prp[,2] > 2 | prp[,1]/prp[,2] < 1/2)
Ti <- D > 0.03 | (prp[,1]/prp[,2] > 3.2 | prp[,1]/prp[,2] < 1/3.2)
text(prp[,1], prp[,2], ifelse(Ti, as.character(x$spp), NA),
    cex=0.6, pos=3, col=1)

Max <- 100*2.1
D <- as.matrix(dist(prt))
diag(D) <- Inf
D <- apply(D, 1, min)
plot(100*prt, xaxs = "i", yaxs = "i", type="n", asp=1,
    ylim=c(0, Max), xlim=c(0, Max),
    xlab="Distance-sampling DD Estimate (m)",
    ylab="LOO DD Estimate (m)")
abline(0,2,lty=2, col="grey")
abline(0,1/2,lty=2, col="grey")
abline(0,1.5,lty=2, col="grey")
abline(0,1/1.5,lty=2, col="grey")
abline(0,1,lty=1, col=1)
r0 <- size_fun(x$logmass, 0.2*Max/50, Max/50)
draw.ellipse(100*prt[,1], 100*prt[,2], r0, r0,
    border=c(Col1,Col3)[as.integer(x$Hab2)])
segments(100*prt[,1], 100*PIt[,1], 100*prt[,1], 100*PIt[,2],
    col=c(Col1,Col3)[as.integer(x$Hab2)], lwd=1)
box()
legend("topleft", pch=21, col=c(Col1,Col3), pt.cex=c(2,2), bty="n",
    legend=c("Habitat: Closed", "Habitat: Open"))
r <- seq(0.2*Max/50, Max/50, len=5)
S <- Max/0.7
draw.ellipse(S*rep(0.06, 5), S*0.55+r-max(r), r, r, border=1)
text(S*c(0.09, 0.06, 0.06), S*c(0.61, 0.52, 0.58),
    c("Body Mass (g)",round(range(x$logmass),1)),
    cex=c(1,0.75,0.75))
text(0.9*Max, 0.05*Max, expression(M[0*beta]-DD))
Ti <- D > 0.09 | (prt[,1]/prt[,2] > 2 | prt[,1]/prt[,2] < 1/2)
text(100*prt[,1], 100*prt[,2], ifelse(Ti, as.character(x$spp), NA),
    cex=0.6, pos=3, col=1)

par(op)
#dev.off()

library(intrval)
## SR
## # species with overlapping PI
sum(prp[,"obs"] %[]% PIp)

## # species with >150% or <66% of original estimates
sum((prp[,"pred"]/prp[,"obs"]) %)(% c(3/2, 2/3))

## # species with >200% or <50% of original estimates
sum((prp[,"pred"]/prp[,"obs"]) %)(% c(2, 1/2))

## DD
## % species with overlapping PI
sum(prt[,"obs"] %[]% PIt)

## # species with >150% or <66% of original estimates
sum((prt[,"pred"]/prt[,"obs"]) %)(% c(3/2, 2/3))

## # species with >200% or <50% of original estimates
sum((prt[,"pred"]/prt[,"obs"]) %)(% c(2, 1/2))

library(phytools)

contMap2 <-
function (tree, x, res = 100, fsize = NULL, ftype = NULL, lwd = 4,
    legend = NULL, lims = NULL, outline = TRUE, sig = 3, type = "phylogram",
    direction = "rightwards", plot = TRUE, col_fun=rainbow, ...)
{
    if (!inherits(tree, "phylo"))
        stop("tree should be an object of class \"phylo\".")
    if (hasArg(mar))
        mar <- list(...)$mar
    else mar <- rep(0.3, 4)
    if (hasArg(offset))
        offset <- list(...)$offset
    else offset <- NULL
    if (hasArg(method))
        method <- list(...)$method
    else method <- "fastAnc"
    if (hasArg(hold))
        hold <- list(...)$hold
    else hold <- TRUE
    if (hasArg(leg.txt))
        leg.txt <- list(...)$leg.txt
    else leg.txt <- "trait value"
    h <- max(nodeHeights(tree))
    steps <- 0:res/res * max(h)
    H <- nodeHeights(tree)
    if (method == "fastAnc")
        a <- fastAnc(tree, x)
    else if (method == "anc.ML") {
        fit <- anc.ML(tree, x)
        a <- fit$ace
        if (!is.null(fit$missing.x))
            x <- c(x, fit$missing.x)
    }
    else if (method == "user") {
        if (hasArg(anc.states))
            anc.states <- list(...)$anc.states
        else {
            cat("No ancestral states have been provided. Using states estimated with fastAnc.\n\n")
            a <- fastAnc(tree, x)
        }
        if (length(anc.states) < tree$Nnode) {
            nodes <- as.numeric(names(anc.states))
            tt <- tree
            for (i in 1:length(nodes)) {
                M <- matchNodes(tt, tree, method = "distances")
                ii <- M[which(M[, 2] == nodes[i]), 1]
                tt <- bind.tip(tt, nodes[i], edge.length = 0,
                  where = ii)
            }
            a <- fastAnc(tt, c(x, anc.states))
            M <- matchNodes(tree, tt, method = "distances")
            a <- a[as.character(M[, 2])]
            names(a) <- M[, 1]
        }
        else {
            if (is.null(names(anc.states)))
                names(anc.states) <- 1:tree$Nnode + Ntip(tree)
            a <- anc.states[as.character(1:tree$Nnode + Ntip(tree))]
        }
    }
    y <- c(a, x[tree$tip.label])
    names(y)[1:Ntip(tree) + tree$Nnode] <- 1:Ntip(tree)
    A <- matrix(y[as.character(tree$edge)], nrow(tree$edge),
        ncol(tree$edge))
    cols <- col_fun(1001)
    names(cols) <- 0:1000
    if (is.null(lims))
        lims <- c(min(y), max(y))
    trans <- 0:1000/1000 * (lims[2] - lims[1]) + lims[1]
    names(trans) <- 0:1000
    tree$maps <- vector(mode = "list", length = nrow(tree$edge))
    for (i in 1:nrow(tree$edge)) {
        XX <- cbind(c(H[i, 1], steps[intersect(which(steps >
            H[i, 1]), which(steps < H[i, 2]))]), c(steps[intersect(which(steps >
            H[i, 1]), which(steps < H[i, 2]))], H[i, 2])) - H[i,
            1]
        YY <- rowMeans(XX)
        if (!all(YY == 0)) {
            b <- vector()
            for (j in 1:length(YY)) b[j] <- (A[i, 1]/YY[j] +
                A[i, 2]/(max(XX) - YY[j]))/(1/YY[j] + 1/(max(XX) -
                YY[j]))
        }
        else b <- A[i, 1]
        d <- sapply(b, phytools:::getState, trans = trans)
        tree$maps[[i]] <- XX[, 2] - XX[, 1]
        names(tree$maps[[i]]) <- d
    }
    tree$mapped.edge <- phytools:::makeMappedEdge(tree$edge, tree$maps)
    tree$mapped.edge <- tree$mapped.edge[, order(as.numeric(colnames(tree$mapped.edge)))]
    class(tree) <- c("simmap", setdiff(class(tree), "simmap"))
    xx <- list(tree = tree, cols = cols, lims = lims)
    class(xx) <- "contMap"
    if (plot)
        phytools:::plot.contMap(xx, fsize = fsize, ftype = ftype, lwd = lwd,
            legend = legend, outline = outline, sig = sig, type = type,
            mar = mar, direction = direction, offset = offset,
            hold = hold, leg.txt = leg.txt)
    invisible(xx)
}

## log singing rates
contMap2(tre, structure(d2$logphi, .Names=NAMES),
    fsize=0.5, outline=FALSE, col_fun=col_fun)

## log detection distances
contMap2(tre, structure(d2$logtau, .Names=NAMES),
    fsize=0.5, outline=FALSE, col_fun=col_fun)

## log body mass
contMap2(tre, structure(d2$logmass, .Names=NAMES),
    fsize=0.5, outline=FALSE, col_fun=col_fun)

## pitch
contMap2(tre, structure(d2$MaxFreqkHz, .Names=NAMES),
    fsize=0.5, outline=FALSE, col_fun=col_fun)

## migratory status
contMap2(tre, structure(d2$Mig2, .Names=NAMES),
    fsize=0.5, outline=FALSE, col_fun=col_fun)

## habitat associations
contMap2(tre, structure(d2$Hab2, .Names=NAMES),
    fsize=0.5, outline=FALSE, col_fun=col_fun)