Results
-#load(system.file("extdata", "lhreg-results-DE.Rdata", package = "lhreg"))Load results and set some values:
+load(system.file("extdata", "lhreg-results-DE.rda", package = "lhreg"))
+
+Crit <- -0.5*qchisq(0.95, 1)
+ltmp <- seq(0, 1, by=0.0001)
+## red-yl-blue
+Col <- c("#2C7BB6", "#6BAACF", "#ABD9E9", "#D4ECD3", "#FFFFBF", "#FED690",
+ "#FDAE61", "#EA633E", "#D7191C")
+Col1 <- Col[1]
+Col2 <- rgb(171/255, 217/255, 233/255, 0.5) # Col[3]
+Col3 <- Col[9]
+Col4 <- rgb(253/255, 174/255, 97/255, 0.5) # Col[7]
+prt <- exp(pr_tMx)
+prp <- exp(pr_pMlx)
+
+
+
+Table with estimates and MSE
+ff <- function(z) {
+ zz <- z$summary
+
+ 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"))
+ }
+ structure(sapply(1:nrow(zz), function(i)
+ paste0(round(zz[i,1], 3), " (SE +/- ", round(zz[i,2], 3), pcut(zz[i,4]), ")")),
+ names=rownames(zz))
+}
+zzz <- lapply(list(pM0, pMp0, pMl0, pMx, pMpx, pMlx), ff)
+m <- matrix("", length(zzz[[6]]), 6)
+rownames(m) <- names(zzz[[6]])
+colnames(m) <- c("M0", "Mp0", "Ml0", "Mx", "Mpx", "Mlx")
+for (i in 1:6) {
+ j <- match(names(zzz[[i]]), rownames(m))
+ m[j,i] <- zzz[[i]]
+}
+ii <- grep("(SE +/- NANA)", m,fixed=TRUE)
+m[ii] <- gsub("(SE +/- NANA)", "(fixed)", m[ii], fixed=TRUE)
+m[m==""] <- "n/a"
+m1 <- m
+
+zzz <- lapply(list(tM0, tMp0, tMl0, tMx, tMpx, tMlx), ff)
+m <- matrix("", length(zzz[[6]]), 6)
+rownames(m) <- names(zzz[[6]])
+colnames(m) <- c("M0", "Mp0", "Ml0", "Mx", "Mpx", "Mlx")
+for (i in 1:6) {
+ j <- match(names(zzz[[i]]), rownames(m))
+ m[j,i] <- zzz[[i]]
+}
+ii <- grep("(SE +/- NANA)", m,fixed=TRUE)
+m[ii] <- gsub("(SE +/- NANA)", "(fixed)", m[ii], fixed=TRUE)
+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## M0 Mp0
+## beta_Intercept -1.449 (SE +/- 0.033***) -1.755 (SE +/- 0.396***)
+## beta_Mig2W n/a n/a
+## beta_xMaxFreqkHz n/a n/a
+## sigma 0.374 (SE +/- 0.024***) 0.877 (SE +/- 0.062***)
+## lambda 0 (fixed) 1 (fixed)
+## df 2 2
+## dAIC 41.061 38.528
+## XV_MSE 0.207 0.169
+## Ml0 Mx
+## beta_Intercept -1.712 (SE +/- 0.229***) -1.627 (SE +/- 0.088***)
+## beta_Mig2W n/a -0.368 (SE +/- 0.08***)
+## beta_xMaxFreqkHz n/a 0.368 (SE +/- 0.122**)
+## sigma 0.536 (SE +/- 0.066***) 0.328 (SE +/- 0.022***)
+## lambda 0.824 (SE +/- 0.075***) 0 (fixed)
+## df 3 4
+## dAIC 1.74 13.017
+## XV_MSE 0.154 0.178
+## Mpx Mlx
+## beta_Intercept -1.77 (SE +/- 0.397***) -1.722 (SE +/- 0.217***)
+## beta_Mig2W -0.167 (SE +/- 0.131ns) -0.199 (SE +/- 0.104+)
+## beta_xMaxFreqkHz 0.174 (SE +/- 0.116ns) 0.19 (SE +/- 0.126ns)
+## sigma 0.863 (SE +/- 0.061***) 0.493 (SE +/- 0.07***)
+## lambda 1 (fixed) 0.771 (SE +/- 0.106***)
+## df 4 5
+## dAIC 38.739 0
+## XV_MSE 0.169 0.153print.default(m2, quote=FALSE) # DD results## M0 Mp0
+## beta_Intercept -0.196 (SE +/- 0.026***) -0.085 (SE +/- 0.345ns)
+## beta_logmass n/a n/a
+## beta_Mig2W n/a n/a
+## beta_xMaxFreqkHz n/a n/a
+## beta_Hab2Open n/a n/a
+## sigma 0.305 (SE +/- 0.019***) 0.762 (SE +/- 0.048***)
+## lambda 0 (fixed) 1 (fixed)
+## df 2 2
+## dAIC 80.685 82.206
+## XV_MSE 0.098 0.081
+## Ml0 Mx
+## beta_Intercept -0.085 (SE +/- 0.167ns) -0.582 (SE +/- 0.107***)
+## beta_logmass n/a 0.158 (SE +/- 0.022***)
+## beta_Mig2W n/a -0.168 (SE +/- 0.056**)
+## beta_xMaxFreqkHz n/a -0.229 (SE +/- 0.088**)
+## beta_Hab2Open n/a 0.145 (SE +/- 0.043***)
+## sigma 0.416 (SE +/- 0.05***) 0.221 (SE +/- 0.014***)
+## lambda 0.747 (SE +/- 0.097***) 0 (fixed)
+## df 3 6
+## dAIC 54.144 0
+## XV_MSE 0.077 0.057
+## Mpx Mlx
+## beta_Intercept -0.814 (SE +/- 0.379*) -0.581 (SE +/- 0.106***)
+## beta_logmass 0.188 (SE +/- 0.049***) 0.158 (SE +/- 0.022***)
+## beta_Mig2W -0.17 (SE +/- 0.102+) -0.17 (SE +/- 0.057**)
+## beta_xMaxFreqkHz -0.183 (SE +/- 0.089*) -0.23 (SE +/- 0.087**)
+## beta_Hab2Open 0.136 (SE +/- 0.049**) 0.145 (SE +/- 0.044***)
+## sigma 0.69 (SE +/- 0.044***) 0.221 (SE +/- 0.014***)
+## lambda 1 (fixed) 0.344 (SE +/- 0.459ns)
+## df 6 7
+## dAIC 62.718 2.003
+## XV_MSE 0.063 0.054Note: estimates and p-values for lambda taken from profile likelihood results due to issues with boundary condition when traits used for DD model.
+
+
+
+
+
+
+
+
+Profile likelihood figures
+#pdf("Fig1.pdf", width=14, height=6)
+op <- par(mfrow=c(1,2))
+
+Res1 <- pl_pMl0
+Res2 <- pl_pMlx
+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)## Max_Ml0 CI1 CI2
+## 0.835 0.627 0.933round(c(Max_Mlx=ltmp[which.max(tmp2)], CI=range(ltmp[i2])), 3)## Max_Mlx CI1 CI2
+## 0.789 0.447 0.921Res1 <- pl_tMl0
+Res2 <- pl_tMlx
+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)## Max_Ml0 CI1 CI2
+## 0.759 0.504 0.899round(c(Max_Mlx=ltmp[which.max(tmp2)], CI=range(ltmp[i2])), 3)## Max_Mlx CI1 CI2
+## 0.000 0.000 0.421par(op)
+
+#dev.off()
+
+
+
LOO cross validation plots
+#pdf("Fig2.pdf", width=14, height=7)
+
+op <- par(mfrow=c(1,2))
+
+Max <- 0.7
+plot(prp, xaxs = "i", yaxs = "i", type="n",
+ ylim=c(0, Max), xlim=c(0, Max),
+ xlab="Time-removal SR Estimate", ylab="LOO SR Estimate")
+abline(0,1,lty=1, col=1)
+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")
+points(prp, xaxs = "i", yaxs = "i",
+ col=c(Col1,Col3)[as.integer(x$Mig2)],
+ cex=0.2+2*x$MaxFreqkHz/10)
+box()
+legend("topleft", pch=21, col=c(Col1,Col3,1,1),
+ pt.cex=c(1.5,1.5,2,1), bty="n",
+ legend=c("Migrant", "Resident", "Song Picth: High", "Song Pitch: Low"))
+text(0.9*Max, 0.05*Max, expression(M[lambda*beta]-SR))
+
+Max <- 2.1
+plot(prt, xaxs = "i", yaxs = "i", type="n",
+ ylim=c(0, Max), xlim=c(0, Max),
+ xlab="Distance-sampling DD Estimate", ylab="LOO DD Estimate")
+abline(0,1,lty=1, col=1)
+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")
+points(prt, xaxs = "i", yaxs = "i",
+ cex=0.2+2*x$logmass/5,
+ col=c(Col1,Col3)[as.integer(x$Hab2)])
+box()
+legend("topleft", pch=21, col=c(Col1,Col3,1,1), pt.cex=c(1.5,1.5,2,1), bty="n",
+ legend=c("Habitat: Closed", "Habitat: Open", "Body Mass: Large", "Body Mass: Small"))
+text(0.9*Max, 0.05*Max, expression(M[0*beta]-DD))
par(op)
+
+#dev.off()