R Under development (unstable) (2025-08-31 r88749 ucrt) -- "Unsuffered Consequences" Copyright (C) 2025 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > #devtools::install_github("psolymos/detect") > library(detect) Loading required package: Formula Loading required package: stats4 Loading required package: pbapply detect 0.5-0 2025-08-31 > > ## --- run examples with \dontrun sections --- > > help_pages <- c(#"bootstrap", + "cmulti", "convertEDR", "databu", "datocc", + "hbootindex", "oven", "svabu", "svocc") > > for (i in help_pages) { + cat("\n\n---------- detect example:", i, "----------\n\n") + eval(parse(text=paste0("example(", i, + ", package = 'detect', run.dontrun = TRUE)"))) + } ---------- detect example: cmulti ---------- cmulti> simfun1 <- function(n = 10, phi = 0.1, c=1, tau=0.8, type="rem") { cmulti+ if (type=="dis") { cmulti+ Dparts <- matrix(c(0.5, 1, NA, cmulti+ 0.5, 1, Inf, cmulti+ 1, Inf, NA), 3, 3, byrow=TRUE) cmulti+ D <- Dparts[sample.int(3, n, replace=TRUE),] cmulti+ CP <- 1-exp(-(D/tau)^2) cmulti+ } else { cmulti+ Dparts <- matrix(c(5, 10, NA, cmulti+ 3, 5, 10, cmulti+ 3, 5, NA), 3, 3, byrow=TRUE) cmulti+ D <- Dparts[sample.int(3, n, replace=TRUE),] cmulti+ CP <- 1-c*exp(-D*phi) cmulti+ } cmulti+ k <- ncol(D) cmulti+ P <- CP - cbind(0, CP[, -k, drop=FALSE]) cmulti+ Psum <- rowSums(P, na.rm=TRUE) cmulti+ PPsum <- P / Psum cmulti+ Pok <- !is.na(PPsum) cmulti+ N <- rpois(n, 10) cmulti+ Y <- matrix(NA, ncol(PPsum), nrow(PPsum)) cmulti+ Ypre <- sapply(1:n, function(i) rmultinom(1, N, PPsum[i,Pok[i,]])) cmulti+ Y[t(Pok)] <- unlist(Ypre) cmulti+ Y <- t(Y) cmulti+ list(Y=Y, D=D) cmulti+ } cmulti> n <- 200 cmulti> x <- rnorm(n) cmulti> X <- cbind(1, x) cmulti> ## removal, constant cmulti> vv <- simfun1(n=n, phi=exp(-1.5)) cmulti> m1 <- cmulti(vv$Y | vv$D ~ 1, type="rem") cmulti> coef(m1) log.phi_(Intercept) -1.471258 cmulti> ## mixture, constant (mix and fmix are identical) cmulti> vv <- simfun1(n=n, phi=exp(-1.5), c=plogis(0.8)) cmulti> m2 <- cmulti(vv$Y | vv$D ~ 1, type="mix") cmulti> coef(m2) log.phi logit.c_(Intercept) -1.4800450 0.9944858 cmulti> m2f <- cmulti(vv$Y | vv$D ~ 1, type="fmix") cmulti> coef(m2f) log.phi_(Intercept) logit.c -1.4800450 0.9944858 cmulti> ## dist, constant cmulti> vv <- simfun1(n=n, tau=exp(-0.2), type="dis") cmulti> m3 <- cmulti(vv$Y | vv$D ~ 1, type="dis") cmulti> coef(m3) log.tau_(Intercept) -0.1976422 cmulti> ## removal, not constant cmulti> log.phi <- crossprod(t(X), c(-2,-1)) cmulti> vv <- simfun1(n=n, phi=exp(cbind(log.phi, log.phi, log.phi))) cmulti> m1 <- cmulti(vv$Y | vv$D ~ x, type="rem") cmulti> coef(m1) log.phi_(Intercept) log.phi_x -2.154161 -1.238920 cmulti> ## mixture, fixed phi, varying c cmulti> logit.c <- crossprod(t(X), c(-2,1)) cmulti> vv <- simfun1(n=n, phi=exp(-1.5), c=plogis(cbind(logit.c, logit.c, logit.c))) cmulti> m2 <- cmulti(vv$Y | vv$D ~ x, type="mix") cmulti> coef(m2) log.phi logit.c_(Intercept) logit.c_x -1.353381 -1.709030 1.066189 cmulti> ## mixture, varying phi, fixed c cmulti> log.phi <- crossprod(t(X), c(-2,-1)) cmulti> vv <- simfun1(n=n, phi=exp(cbind(log.phi, log.phi, log.phi)), c=plogis(0.8)) cmulti> m2f <- cmulti(vv$Y | vv$D ~ x, type="fmix") cmulti> coef(m2f) log.phi_(Intercept) log.phi_x logit.c -2.0405233 -1.0961060 0.7000887 cmulti> ## dist, not constant cmulti> log.tau <- crossprod(t(X), c(-0.5,-0.2)) cmulti> vv <- simfun1(n=n, tau=exp(cbind(log.tau, log.tau, log.tau)), type="dis") cmulti> m3 <- cmulti(vv$Y | vv$D ~ x, type="dis") cmulti> coef(m3) log.tau_(Intercept) log.tau_x -0.5385985 -0.1992466 cmulti> summary(m3) Call: cmulti(formula = vv$Y | vv$D ~ x, type = "dis") Distance Sampling (half-normal, circular area) Conditional Maximum Likelihood estimates Coefficients: Estimate Std. Error z value Pr(>|z|) log.tau_(Intercept) -0.53860 0.01958 -27.513 <2e-16 *** log.tau_x -0.19925 0.02003 -9.947 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Log-likelihood: -314.6 BIC = 639.7 cmulti> coef(m3) log.tau_(Intercept) log.tau_x -0.5385985 -0.1992466 cmulti> vcov(m3) log.tau_(Intercept) log.tau_x log.tau_(Intercept) 3.832226e-04 5.832549e-06 log.tau_x 5.832549e-06 4.011977e-04 cmulti> AIC(m3) [1] 633.1022 cmulti> confint(m3) 2.5 % 97.5 % log.tau_(Intercept) -0.5769669 -0.5002301 log.tau_x -0.2385045 -0.1599887 cmulti> logLik(m3) 'log Lik.' -314.5511 (df=2) cmulti> ## fitted values cmulti> plot(exp(log.tau), fitted(m3)) cmulti> ## prediction for new locations (type = 'rem') cmulti> ndf <- data.frame(x=seq(-1, 1, by=0.1)) cmulti> summary(predict(m1, newdata=ndf, type="link")) Min. 1st Qu. Median Mean 3rd Qu. Max. -3.3931 -2.7736 -2.1542 -2.1542 -1.5347 -0.9152 cmulti> summary(pr1 <- predict(m1, newdata=ndf, type="response")) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.03360 0.06244 0.11600 0.15150 0.21552 0.40042 cmulti> ## turing singing rates into probabilities requires total duration cmulti> ## 5 minutes used here cmulti> psing <- 1-exp(-5*pr1) cmulti> plot(ndf$x, psing, type="l", ylim=c(0,1)) cmulti> ## prediction for new locations (type = 'dis') cmulti> summary(predict(m3, newdata=ndf, type="link")) Min. 1st Qu. Median Mean 3rd Qu. Max. -0.7378 -0.6382 -0.5386 -0.5386 -0.4390 -0.3394 cmulti> summary(pr3 <- predict(m3, newdata=ndf, type="response")) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.4781 0.5282 0.5836 0.5878 0.6447 0.7122 cmulti> ## turing EDR into probabilities requires finite truncation distances cmulti> ## r=0.5 used here (50 m) cmulti> r <- 0.5 cmulti> pdet <- pr3^2*(1-exp(-r^2/pr3^2))/r^2 cmulti> plot(ndf$x, pdet, type="l", ylim=c(0,1)) cmulti> ## joint removal-distance estimation cmulti> ## is not different from 2 orthogonal estimations cmulti> cmulti> simfun12 <- function(n = 10, phi = 0.1, c=1, tau=0.8, type="rem") { cmulti+ Flat <- function(x, DIM, dur=TRUE) { cmulti+ x <- array(x, DIM) cmulti+ if (!dur) { cmulti+ x <- aperm(x,c(1,3,2)) cmulti+ } cmulti+ dim(x) <- c(DIM[1], DIM[2]*DIM[3]) cmulti+ x cmulti+ } cmulti+ Dparts1 <- matrix(c(5, 10, NA, cmulti+ 3, 5, 10, cmulti+ 3, 5, NA), 3, 3, byrow=TRUE) cmulti+ D1 <- Dparts1[sample.int(3, n, replace=TRUE),] cmulti+ CP1 <- 1-c*exp(-D1*phi) cmulti+ Dparts2 <- matrix(c(0.5, 1, NA, cmulti+ 0.5, 1, Inf, cmulti+ 1, Inf, NA), 3, 3, byrow=TRUE) cmulti+ D2 <- Dparts2[sample.int(3, n, replace=TRUE),] cmulti+ CP2 <- 1-exp(-(D2/tau)^2) cmulti+ k1 <- ncol(D1) cmulti+ k2 <- ncol(D2) cmulti+ DIM <- c(n, k1, k2) cmulti+ P1 <- CP1 - cbind(0, CP1[, -k1, drop=FALSE]) cmulti+ P2 <- CP2 - cbind(0, CP2[, -k2, drop=FALSE]) cmulti+ Psum1 <- rowSums(P1, na.rm=TRUE) cmulti+ Psum2 <- rowSums(P2, na.rm=TRUE) cmulti+ Pflat <- Flat(P1, DIM, dur=TRUE) * Flat(P2, DIM, dur=FALSE) cmulti+ PsumFlat <- Psum1 * Psum2 cmulti+ PPsumFlat <- Pflat / PsumFlat cmulti+ PokFlat <- !is.na(PPsumFlat) cmulti+ N <- rpois(n, 10) cmulti+ Yflat <- matrix(NA, ncol(PPsumFlat), nrow(PPsumFlat)) cmulti+ YpreFlat <- sapply(1:n, function(i) rmultinom(1, N, PPsumFlat[i,PokFlat[i,]])) cmulti+ Yflat[t(PokFlat)] <- unlist(YpreFlat) cmulti+ Yflat <- t(Yflat) cmulti+ Y <- array(Yflat, DIM) cmulti+ k1 <- dim(Y)[2] cmulti+ k2 <- dim(Y)[3] cmulti+ Y1 <- t(sapply(1:n, function(i) { cmulti+ count <- rowSums(Y[i,,], na.rm=TRUE) cmulti+ nas <- rowSums(is.na(Y[i,,])) cmulti+ count[nas == k2] <- NA cmulti+ count cmulti+ })) cmulti+ Y2 <- t(sapply(1:n, function(i) { cmulti+ count <- colSums(Y[i,,], na.rm=TRUE) cmulti+ nas <- colSums(is.na(Y[i,,])) cmulti+ count[nas == k2] <- NA cmulti+ count cmulti+ })) cmulti+ list(Y=Y, D1=D1, D2=D2, Y1=Y1, Y2=Y2) cmulti+ } cmulti> ## removal and distance, constant cmulti> vv <- simfun12(n=n, phi=exp(-1.5), tau=exp(-0.2)) cmulti> res <- cmulti2.fit(vv$Y, vv$D1, vv$D2) cmulti> res1 <- cmulti.fit(vv$Y1, vv$D1, NULL, "rem") cmulti> res2 <- cmulti.fit(vv$Y2, vv$D2, NULL, "dis") cmulti> ## points estimates are identical cmulti> cbind(res$coef, c(res1$coef, res2$coef)) [,1] [,2] [1,] -1.5315358 -1.5315358 [2,] -0.2335382 -0.2335382 cmulti> ## standard errors are identical cmulti> cbind(sqrt(diag(res$vcov)), cmulti+ c(sqrt(diag(res1$vcov)),sqrt(diag(res2$vcov)))) [,1] [,2] [1,] 0.07224423 0.07224423 [2,] 0.02204826 0.02204826 cmulti> ## removal and distance, not constant cmulti> vv <- simfun12(n=n, cmulti+ phi=exp(cbind(log.phi, log.phi, log.phi)), cmulti+ tau=exp(cbind(log.tau, log.tau, log.tau))) cmulti> res <- cmulti2.fit(vv$Y, vv$D1, vv$D2, X1=X, X2=X) cmulti> res1 <- cmulti.fit(vv$Y1, vv$D1, X, "rem") cmulti> res2 <- cmulti.fit(vv$Y2, vv$D2, X, "dis") cmulti> ## points estimates are identical cmulti> cbind(res$coef, c(res1$coef, res2$coef)) [,1] [,2] [1,] -2.0135266 -2.0131529 [2,] -0.9446835 -0.9443971 [3,] -0.4884864 -0.4885030 [4,] -0.1712155 -0.1712501 cmulti> ## standard errors are identical cmulti> cbind(sqrt(diag(res$vcov)), cmulti+ c(sqrt(diag(res1$vcov)),sqrt(diag(res2$vcov)))) [,1] [,2] [1,] 0.10200120 0.10194772 [2,] 0.09519510 0.09515732 [3,] 0.01650040 0.01650027 [4,] 0.01680185 0.01680166 ---------- detect example: convertEDR ---------- cnvEDR> convertEDR(1, 0.5, truncated=FALSE) [1] 0.4703182 cnvEDR> ## should be close to 1 cnvEDR> convertEDR(convertEDR(1, 0.5, truncated=FALSE), 0.5, truncated=TRUE) [1] 1.000002 ---------- detect example: databu ---------- databu> data(databu) databu> str(databu) 'data.frame': 1000 obs. of 12 variables: $ N : num 8 12 3 3 0 0 11 8 9 5 ... $ Y : num 7 2 3 2 0 0 0 8 6 4 ... $ x1 : num 0.114 0.622 0.609 0.623 0.861 ... $ x2 : num 0.985 -1.225 0.71 -0.109 1.783 ... $ x3 : num 0.927 -0.586 -0.828 -0.568 -0.521 ... $ x4 : num -0.6699 0.4302 -0.0794 0.2769 -0.0882 ... $ x5 : num 1 1 0 0 1 0 1 1 1 0 ... $ x6 : num 0 0 1 0 1 0 1 1 1 1 ... $ p : num 0.922 0.125 0.918 0.686 0.983 ... $ lambda: num 11.12 7.41 4.54 4.49 6.12 ... $ A : num 1 1 1 1 0 0 1 1 1 1 ... $ phi : num 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 ... databu> ## simulation databu> n <- 1000 databu> set.seed(1234) databu> x1 <- runif(n,0,1) databu> x2 <- rnorm(n,0,1) databu> x3 <- runif(n,-1,1) databu> x4 <- runif(n,-1,1) databu> x5 <- rbinom(n,1,0.6) databu> x6 <- rbinom(n,1,0.4) databu> x7 <- rnorm(n,0,1) databu> X <- model.matrix(~ x1 + x5) databu> Z <- model.matrix(~ x2 + x5) databu> Q <- model.matrix(~ x7) databu> beta <- c(2,-0.8,0.5) databu> theta <- c(1, 2, -0.5) databu> phi <- 0.3 databu> p <- drop(binomial("logit")$linkinv(Z %*% theta)) databu> lambda <- drop(exp(X %*% beta)) databu> A <- rbinom(n, 1, 1-phi) databu> N <- rpois(n, lambda * A) databu> Y <- rbinom(n, N, p) databu> databu <- data.frame(N=N, Y=Y, x1, x2, x3, x4, x5, x6, p=p, lambda=lambda, A, phi) ---------- detect example: datocc ---------- datocc> data(datocc) datocc> str(datocc) 'data.frame': 1000 obs. of 8 variables: $ Y : num 1 1 1 1 1 1 0 0 0 1 ... $ W : num 1 0 0 0 1 1 0 0 0 1 ... $ x1 : num -0.773 0.245 0.219 0.247 0.722 ... $ x2 : Factor w/ 2 levels "0","1": 2 1 1 1 2 1 1 1 2 2 ... $ x3 : num -1.205 0.301 -1.539 0.635 0.703 ... $ x4 : num -0.9738 -0.0996 -0.1107 1.1922 -1.6559 ... $ p.occ: num 0.71 0.872 0.869 0.873 0.927 ... $ p.det: num 0.605 0.591 0.457 0.615 0.562 ... datocc> ## simulation datocc> n <- 1000 datocc> set.seed(1234) datocc> x1 <- runif(n, -1, 1) datocc> x2 <- as.factor(rbinom(n, 1, 0.5)) datocc> x3 <- rnorm(n) datocc> x4 <- rnorm(n) datocc> beta <- c(0.6, 0.5) datocc> theta <- c(0.4, -0.5, 0.3) datocc> X <- model.matrix(~ x1) datocc> Z <- model.matrix(~ x1 + x3) datocc> mu <- drop(X %*% beta) datocc> nu <- drop(Z %*% theta) datocc> p.occ <- binomial("cloglog")$linkinv(mu) datocc> p.det <- binomial("logit")$linkinv(nu) datocc> Y <- rbinom(n, 1, p.occ) datocc> W <- rbinom(n, 1, Y * p.det) datocc> datocc <- data.frame(Y, W, x1, x2, x3, x4, p.occ, p.det) ---------- detect example: hbootindex ---------- hbtndx> ## equal group sizes hbtndx> groups <- rep(1:4, each=5) hbtndx> strata <- rep(1:2, each=10) hbtndx> hbootindex(groups, strata, 3) [,1] [,2] [,3] [,4] [1,] 5 1 4 12 [2,] 7 9 17 6 [3,] 3 3 5 3 [4,] 19 7 6 8 [5,] 6 9 13 5 [6,] 17 18 1 12 [7,] 4 11 1 20 [8,] 20 8 7 19 [9,] 2 13 6 16 [10,] 13 18 5 20 [11,] 12 5 12 3 [12,] 16 11 11 17 [13,] 1 4 11 3 [14,] 14 13 20 6 [15,] 9 16 17 8 [16,] 10 18 13 14 [17,] 18 10 19 9 [18,] 11 17 18 19 [19,] 15 13 3 19 [20,] 8 8 7 1 hbtndx> ## unequal group sizes hbtndx> groups <- groups[-c(5,9,10,11)] hbtndx> strata <- strata[-c(5,9,10,11)] hbtndx> hbootindex(groups, strata, 3) [,1] [,2] [,3] [,4] [1,] 15 14 13 8 [2,] 3 1 11 15 [3,] 7 16 1 16 [4,] 9 8 15 2 [5,] 16 1 7 16 [6,] 5 16 15 2 [7,] 1 7 11 4 [8,] 6 1 4 7 [9,] 8 1 8 16 [10,] 12 16 4 8 [11,] 13 15 8 2 [12,] 11 5 14 15 [13,] 10 14 12 6 [14,] 14 16 2 16 [15,] 4 16 7 7 [16,] 2 3 5 8 ---------- detect example: oven ---------- oven> data(oven) oven> str(oven) 'data.frame': 891 obs. of 11 variables: $ count : int 1 0 0 1 0 0 0 0 0 0 ... $ route : int 2 2 2 2 2 2 2 2 2 2 ... $ stop : int 2 4 6 8 10 12 14 16 18 20 ... $ pforest: num 0.947 0.903 0.814 0.89 0.542 ... $ pdecid : num 0.575 0.562 0.549 0.679 0.344 ... $ pagri : num 0 0 0 0 0.414 ... $ long : num 609343 608556 607738 607680 607944 ... $ lat : num 5949071 5947735 5946301 5944720 5943088 ... $ observ : Factor w/ 4 levels "ARS","DW","RDW",..: 4 4 4 4 4 4 4 4 4 4 ... $ julian : int 181 181 181 181 181 181 181 181 181 181 ... $ timeday: int 2 4 6 8 10 12 14 16 18 20 ... ---------- detect example: svabu ---------- svabu> data(databu) svabu> ## fit BZIP and BP models svabu> m00 <- svabu(Y ~ x1 + x5 | x2 + x5, databu[1:200,]) svabu> ## print method svabu> m00 Call: svabu(formula = Y ~ x1 + x5 | x2 + x5, data = databu[1:200, ]) Single visit Binomial - Zero Inflated Poisson model Conditional Maximum Likelihood estimates Coefficients for abundance (log link): (Intercept) x1 x5 2.0318 -0.8921 0.5385 Coefficients for detection (logit link): (Intercept) x2 x5 0.8124 1.6873 -0.4561 Coefficients for zero inflation (logit link): (Intercept) -0.7758 svabu> ## summary: CMLE svabu> summary(m00) Call: svabu(formula = Y ~ x1 + x5 | x2 + x5, data = databu[1:200, ]) Single visit Binomial - Zero Inflated Poisson model Conditional Maximum Likelihood estimates Coefficients for abundance (log link): Estimate Std. Error z value Pr(>|z|) (Intercept) 2.0318 0.1901 10.686 < 2e-16 *** x1 -0.8921 0.1701 -5.244 1.57e-07 *** x5 0.5385 0.1720 3.131 0.00174 ** Coefficients for detection (logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) 0.8124 0.6837 1.188 0.235 x2 1.6873 0.4040 4.177 2.96e-05 *** x5 -0.4561 0.5659 -0.806 0.420 Coefficients for zero inflation (logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) -0.7758 0.1705 -4.551 5.33e-06 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Log-likelihood: -300.1 on 7 Df AIC = 614.2 svabu> ## coef svabu> coef(m00) sta_(Intercept) sta_x1 sta_x5 det_(Intercept) det_x2 2.0317779 -0.8921099 0.5384850 0.8123504 1.6872739 det_x5 zif_(Intercept) -0.4560919 -0.7758071 svabu> coef(m00, model="sta") ## state (abundance) (Intercept) x1 x5 2.0317779 -0.8921099 0.5384850 svabu> coef(m00, model="det") ## detection (Intercept) x2 x5 0.8123504 1.6872739 -0.4560919 svabu> coef(m00, model="zif") ## zero inflation (this is part of the 'true state'!) (Intercept) -0.7758071 svabu> ## Diagnostics and model comparison svabu> svabu> m01 <- svabu(Y ~ x1 + x5 | x2 + x5, databu[1:200,], zeroinfl=FALSE) svabu> ## compare estimates (note, zero inflation is on the logit scale!) svabu> cbind(truth=c(2,-0.8,0.5, 1,2,-0.5, plogis(0.3)), svabu+ "B-ZIP"=coef(m00), "B-P"=c(coef(m01), NA)) truth B-ZIP B-P sta_(Intercept) 2.0000000 2.0317779 1.4824442 sta_x1 -0.8000000 -0.8921099 -0.8740214 sta_x5 0.5000000 0.5384850 0.4810412 det_(Intercept) 1.0000000 0.8123504 2.1676241 det_x2 2.0000000 1.6872739 2.6505899 det_x5 -0.5000000 -0.4560919 -0.7543411 zif_(Intercept) 0.5744425 -0.7758071 NA svabu> ## fitted svabu> plot(fitted(m00), fitted(m01)) svabu> abline(0,1) svabu> ## compare models svabu> AIC(m00, m01) df AIC m00 7 614.2221 m01 6 947.3850 svabu> BIC(m00, m01) df BIC m00 7 637.3103 m01 6 967.1749 svabu> logLik(m00) 'log Lik.' -300.111 (df=7) svabu> logLik(m01) 'log Lik.' -467.6925 (df=6) svabu> ## diagnostic plot svabu> plot(m00) svabu> plot(m01) svabu> ## Bootstrap svabu> svabu> ## non parametric bootstrap svabu> ## - initial values are the estimates svabu> m02 <- bootstrap(m00, B=25) svabu> attr(m02, "bootstrap") cfs sta_(Intercept) 2.0317779 2.1383662 1.8796803 1.6819710 1.5861799 sta_x1 -0.8921099 -0.9427072 -0.6598919 -0.6442896 -1.1359465 sta_x5 0.5384850 0.3942039 0.4367155 0.7383676 1.1267804 det_(Intercept) 0.8123504 0.9834574 2.2635261 1.8075570 10.4389538 det_x2 1.6872739 1.6314210 2.2257101 2.1404453 1.6394783 det_x5 -0.4560919 -0.6533823 -1.5243790 -1.0736859 -10.2329269 zif_(Intercept) -0.7758071 -0.8721769 -0.9018758 -0.6327677 -0.7527754 sta_(Intercept) 1.9995871 1.7272861 2.0611453 1.7662203 1.9494487 sta_x1 -1.0184141 -1.0096868 -1.1392094 -0.8014978 -0.6421667 sta_x5 0.5535334 0.9690049 0.6720722 0.6413081 0.4455286 det_(Intercept) 1.1898221 1.9102352 1.3400181 2.7275531 1.1698849 det_x2 2.2020826 1.7481780 2.1859442 2.3511040 2.0022592 det_x5 -0.5211715 -1.5791469 -0.8640776 -2.0106466 -0.4654308 zif_(Intercept) -1.0861300 -0.8976331 -0.8939965 -0.7506540 -0.9912754 sta_(Intercept) 1.7449877 1.6972984 2.2009132 1.6861153 1.9477495 sta_x1 -0.7891871 -0.4998508 -1.1388431 -0.7905349 -0.8207768 sta_x5 0.5516090 0.7031281 0.5419532 0.7199541 0.4689128 det_(Intercept) 2.5885119 1.5349872 0.9258446 1.7829465 1.4229612 det_x2 2.6277958 1.7666033 1.6604431 1.8111820 2.1021413 det_x5 -1.1682117 -1.4348129 -0.7357502 -1.1625406 -0.5481337 zif_(Intercept) -0.7539876 -0.6527690 -0.7124697 -0.8952087 -1.0709775 sta_(Intercept) 1.9781782 2.1544208 1.7401016 1.8560490 1.9858931 sta_x1 -0.8845973 -0.8656764 -0.5393025 -0.7888543 -1.0647061 sta_x5 0.6042620 0.5859150 0.6072531 0.5510574 0.5257814 det_(Intercept) 1.1839014 0.4869570 2.1344891 1.3112244 1.2139582 det_x2 1.8342923 1.5400580 2.2910399 2.0773597 2.2423907 det_x5 -0.5899616 -0.5141425 -1.3170923 -0.3503465 -0.2570269 zif_(Intercept) -0.7097423 -1.0285505 -1.0224896 -0.7630790 -0.7721705 sta_(Intercept) 2.1168678 2.1203771 1.7808835 1.9028873 2.02257546 sta_x1 -0.9846326 -0.9296834 -0.6105223 -1.0088491 -1.03486008 sta_x5 0.7337335 0.5581679 0.7246375 0.7044257 0.34033557 det_(Intercept) 0.8169417 0.8600649 1.7573975 2.1915883 1.34243453 det_x2 1.4769239 1.6025528 1.9637932 2.1248886 2.30982455 det_x5 -0.9494331 -0.7676649 -1.4022372 -1.7115860 -0.03189795 zif_(Intercept) -0.8767227 -0.7785346 -0.8479329 -0.7545937 -1.03946033 sta_(Intercept) 1.9831912 sta_x1 -0.9836951 sta_x5 0.5180053 det_(Intercept) 1.8224100 det_x2 2.4151031 det_x5 -1.0040934 zif_(Intercept) -1.1015855 attr(,"type") [1] "nonpar" attr(,"ini") sta_(Intercept) sta_x1 sta_x5 det_(Intercept) det_x2 2.0317779 -0.8921099 0.5384850 0.8123504 1.6872739 det_x5 zif_(Intercept) -0.4560919 -0.7758071 svabu> extractBOOT(m02) $coefficients sta_(Intercept) sta_x1 sta_x5 det_(Intercept) det_x2 1.9130828 -0.8700189 0.6136589 1.8469222 1.9869342 det_x5 zif_(Intercept) -1.2817643 -0.8590525 $std.error sta_(Intercept) sta_x1 sta_x5 det_(Intercept) det_x2 0.1724384 0.1859318 0.1685220 1.8427828 0.3103016 det_x5 zif_(Intercept) 1.8913275 0.1385545 $vcov sta_(Intercept) sta_x1 sta_x5 det_(Intercept) sta_(Intercept) 0.029735000 -0.015099707 -0.017148505 -0.18603591 sta_x1 -0.015099707 0.034570631 -0.004016680 -0.05811423 sta_x5 -0.017148505 -0.004016680 0.028399650 0.20154432 det_(Intercept) -0.186035910 -0.058114229 0.201544316 3.39584849 det_x2 -0.015592927 0.009939520 -0.017052866 -0.00190986 det_x5 0.168945862 0.067831204 -0.229882422 -3.41255910 zif_(Intercept) -0.006578008 0.001907534 0.008227144 0.04722369 det_x2 det_x5 zif_(Intercept) sta_(Intercept) -0.015592927 0.16894586 -0.006578008 sta_x1 0.009939520 0.06783120 0.001907534 sta_x5 -0.017052866 -0.22988242 0.008227144 det_(Intercept) -0.001909860 -3.41255910 0.047223688 det_x2 0.096287098 0.10621726 -0.008865488 det_x5 0.106217262 3.57711975 -0.057667366 zif_(Intercept) -0.008865488 -0.05766737 0.019197360 $cor sta_(Intercept) sta_x1 sta_x5 det_(Intercept) sta_(Intercept) 1.0000000 -0.47095661 -0.5901139 -0.585448447 sta_x1 -0.4709566 1.00000000 -0.1281909 -0.169611258 sta_x5 -0.5901139 -0.12819086 1.0000000 0.648992848 det_(Intercept) -0.5854484 -0.16961126 0.6489928 1.000000000 det_x2 -0.2914135 0.17227716 -0.3261045 -0.003339976 det_x5 0.5180204 0.19288974 -0.7212445 -0.979127593 zif_(Intercept) -0.2753212 0.07404535 0.3523480 0.184954539 det_x2 det_x5 zif_(Intercept) sta_(Intercept) -0.291413469 0.5180204 -0.27532122 sta_x1 0.172277157 0.1928897 0.07404535 sta_x5 -0.326104538 -0.7212445 0.35234802 det_(Intercept) -0.003339976 -0.9791276 0.18495454 det_x2 1.000000000 0.1809857 -0.20620436 det_x5 0.180985722 1.0000000 -0.22006077 zif_(Intercept) -0.206204357 -0.2200608 1.00000000 $B [1] 25 $type [1] "nonpar" $model [1] "full" svabu> summary(m02) Call: svabu(formula = Y ~ x1 + x5 | x2 + x5, data = databu[1:200, ]) Single visit Binomial - Zero Inflated Poisson model Conditional Maximum Likelihood estimates with Nonparametric Bootstrap standard errors (B = 25) Coefficients for abundance (log link): Estimate Std. Error z value Pr(>|z|) (Intercept) 2.0318 0.1724 11.783 < 2e-16 *** x1 -0.8921 0.1859 -4.798 1.6e-06 *** x5 0.5385 0.1685 3.195 0.0014 ** Coefficients for detection (logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) 0.8124 1.8428 0.441 0.659 x2 1.6873 0.3103 5.438 5.4e-08 *** x5 -0.4561 1.8913 -0.241 0.809 Coefficients for zero inflation (logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) -0.7758 0.1386 -5.599 2.15e-08 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Log-likelihood: -300.1 on 7 Df AIC = 614.2 svabu> summary(m02, type="cmle") Call: svabu(formula = Y ~ x1 + x5 | x2 + x5, data = databu[1:200, ]) Single visit Binomial - Zero Inflated Poisson model Conditional Maximum Likelihood estimates Coefficients for abundance (log link): Estimate Std. Error z value Pr(>|z|) (Intercept) 2.0318 0.1901 10.686 < 2e-16 *** x1 -0.8921 0.1701 -5.244 1.57e-07 *** x5 0.5385 0.1720 3.131 0.00174 ** Coefficients for detection (logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) 0.8124 0.6837 1.188 0.235 x2 1.6873 0.4040 4.177 2.96e-05 *** x5 -0.4561 0.5659 -0.806 0.420 Coefficients for zero inflation (logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) -0.7758 0.1705 -4.551 5.33e-06 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Log-likelihood: -300.1 on 7 Df AIC = 614.2 svabu> summary(m02, type="boot") Call: svabu(formula = Y ~ x1 + x5 | x2 + x5, data = databu[1:200, ]) Single visit Binomial - Zero Inflated Poisson model Conditional Maximum Likelihood estimates with Nonparametric Bootstrap standard errors (B = 25) Coefficients for abundance (log link): Estimate Std. Error z value Pr(>|z|) (Intercept) 2.0318 0.1724 11.783 < 2e-16 *** x1 -0.8921 0.1859 -4.798 1.6e-06 *** x5 0.5385 0.1685 3.195 0.0014 ** Coefficients for detection (logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) 0.8124 1.8428 0.441 0.659 x2 1.6873 0.3103 5.438 5.4e-08 *** x5 -0.4561 1.8913 -0.241 0.809 Coefficients for zero inflation (logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) -0.7758 0.1386 -5.599 2.15e-08 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Log-likelihood: -300.1 on 7 Df AIC = 614.2 svabu> ## vcov svabu> vcov(m02, type="cmle") sta_(Intercept) sta_x1 sta_x5 det_(Intercept) sta_(Intercept) 0.03615400 -0.010922993 -0.020094384 -0.108651587 sta_x1 -0.01092299 0.028937607 -0.002478626 -0.002865438 sta_x5 -0.02009438 -0.002478626 0.029584244 0.062541639 det_(Intercept) -0.10865159 -0.002865438 0.062541639 0.467412118 det_x2 -0.04641158 0.003255012 0.005259598 0.194510127 det_x5 0.06857974 0.008270576 -0.078967102 -0.303952040 zif_(Intercept) NA NA NA NA det_x2 det_x5 zif_(Intercept) sta_(Intercept) -0.046411576 0.068579744 NA sta_x1 0.003255012 0.008270576 NA sta_x5 0.005259598 -0.078967102 NA det_(Intercept) 0.194510127 -0.303952040 NA det_x2 0.163193431 -0.056846362 NA det_x5 -0.056846362 0.320236448 NA zif_(Intercept) NA NA 0.02905757 svabu> vcov(m02, type="boot") sta_(Intercept) sta_x1 sta_x5 det_(Intercept) sta_(Intercept) 0.029735000 -0.015099707 -0.017148505 -0.18603591 sta_x1 -0.015099707 0.034570631 -0.004016680 -0.05811423 sta_x5 -0.017148505 -0.004016680 0.028399650 0.20154432 det_(Intercept) -0.186035910 -0.058114229 0.201544316 3.39584849 det_x2 -0.015592927 0.009939520 -0.017052866 -0.00190986 det_x5 0.168945862 0.067831204 -0.229882422 -3.41255910 zif_(Intercept) -0.006578008 0.001907534 0.008227144 0.04722369 det_x2 det_x5 zif_(Intercept) sta_(Intercept) -0.015592927 0.16894586 -0.006578008 sta_x1 0.009939520 0.06783120 0.001907534 sta_x5 -0.017052866 -0.22988242 0.008227144 det_(Intercept) -0.001909860 -3.41255910 0.047223688 det_x2 0.096287098 0.10621726 -0.008865488 det_x5 0.106217262 3.57711975 -0.057667366 zif_(Intercept) -0.008865488 -0.05766737 0.019197360 svabu> vcov(m02, model="sta") (Intercept) x1 x5 (Intercept) 0.02973500 -0.01509971 -0.01714851 x1 -0.01509971 0.03457063 -0.00401668 x5 -0.01714851 -0.00401668 0.02839965 svabu> vcov(m02, model="det") (Intercept) x2 x5 (Intercept) 3.39584849 -0.00190986 -3.4125591 x2 -0.00190986 0.09628710 0.1062173 x5 -3.41255910 0.10621726 3.5771197 svabu> ## confint svabu> confint(m02, type="cmle") ## Wald-type 2.5% 97.5% sta_(Intercept) 1.6591063 2.4044494 sta_x1 -1.2255205 -0.5586993 sta_x5 0.2013698 0.8756002 det_(Intercept) -0.5276289 2.1523296 det_x2 0.8955032 2.4790446 det_x5 -1.5652245 0.6530407 zif_(Intercept) -1.1099081 -0.4417062 svabu> confint(m02, type="boot") ## quantile based 2.5% 97.5% sta_(Intercept) 1.6460494 2.1718555 sta_x1 -1.1389805 -0.5245081 sta_x5 0.3740033 1.0281707 det_(Intercept) 0.6903278 5.6193283 det_x2 1.5163827 2.4948629 det_x5 -5.0940017 -0.1726035 zif_(Intercept) -1.0919258 -0.6452685 svabu> ## parametric bootstrap svabu> simulate(m00, 5) sim_1 sim_2 sim_3 sim_4 sim_5 1 10 16 0 7 9 2 0 1 0 2 1 3 0 3 6 7 3 4 0 0 4 4 0 5 6 5 0 3 0 6 0 4 0 2 0 7 0 0 0 4 1 8 7 8 5 0 0 9 5 3 0 9 4 10 2 3 3 0 2 11 0 0 3 1 0 12 3 0 2 2 1 13 7 0 0 4 4 14 8 0 5 0 1 15 1 2 3 0 2 16 1 2 0 0 0 17 0 0 2 0 0 18 5 7 9 9 8 19 5 6 0 0 0 20 11 0 8 10 11 21 7 0 8 7 5 22 0 0 7 7 3 23 0 7 3 0 13 24 0 6 5 7 0 25 4 2 0 5 5 26 1 3 0 1 2 27 4 0 7 2 4 28 0 0 0 1 2 29 2 5 4 4 4 30 5 6 6 0 7 31 3 0 7 3 6 32 12 11 6 0 0 33 6 4 4 3 8 34 13 10 8 0 0 35 0 7 6 9 8 36 1 0 0 0 2 37 0 2 4 0 4 38 0 0 7 0 11 39 7 1 0 0 1 40 0 0 0 4 0 41 2 5 5 0 5 42 5 3 8 13 5 43 3 6 2 4 4 44 5 4 4 4 5 45 0 0 0 3 0 46 0 0 0 0 0 47 0 0 0 0 1 48 0 6 4 4 0 49 0 13 7 0 9 50 5 8 5 0 5 51 0 0 0 2 0 52 0 0 0 0 0 53 1 1 0 0 0 54 6 1 2 0 1 55 7 2 6 9 7 56 4 0 0 7 6 57 1 1 1 2 3 58 0 0 0 0 0 59 6 10 0 13 0 60 3 4 3 4 3 61 2 2 0 1 0 62 0 2 2 1 5 63 2 7 5 0 0 64 0 0 8 5 5 65 4 6 5 5 7 66 7 3 3 7 7 67 0 0 4 3 0 68 4 4 0 0 5 69 0 8 8 8 9 70 0 4 2 0 0 71 4 4 0 0 4 72 0 0 0 0 1 73 11 4 6 6 6 74 0 2 3 0 4 75 2 2 0 3 0 76 0 1 4 1 0 77 0 6 5 0 5 78 0 2 0 2 1 79 0 0 0 6 4 80 0 3 0 2 0 81 1 0 0 0 0 82 6 8 5 0 0 83 7 0 1 4 0 84 3 4 0 1 3 85 4 4 8 0 7 86 1 0 1 0 0 87 0 0 0 0 0 88 0 5 14 6 5 89 3 6 6 7 8 90 2 0 1 0 1 91 0 7 7 6 8 92 3 0 0 0 4 93 5 3 2 0 8 94 0 0 12 0 0 95 12 0 0 0 11 96 0 3 3 0 0 97 4 2 5 6 2 98 0 0 0 0 7 99 3 0 0 3 5 100 1 0 1 2 0 101 13 4 10 12 6 102 6 0 0 0 2 103 0 0 0 0 0 104 0 5 5 5 7 105 0 2 3 7 4 106 3 2 3 5 2 107 5 12 9 6 8 108 4 2 0 0 7 109 1 0 0 1 0 110 10 13 0 15 12 111 3 2 0 2 2 112 6 0 4 0 8 113 1 0 4 0 2 114 0 5 0 0 6 115 0 0 1 0 0 116 3 2 0 1 0 117 7 0 0 0 8 118 0 3 1 0 0 119 3 2 3 0 4 120 0 2 0 2 1 121 4 3 0 0 4 122 4 1 2 1 3 123 0 0 0 0 0 124 2 5 0 0 0 125 0 0 7 0 4 126 0 5 11 13 11 127 1 0 1 1 1 128 0 3 0 0 0 129 4 3 4 5 4 130 0 0 1 3 2 131 0 2 0 4 6 132 5 0 3 3 4 133 2 2 0 1 0 134 0 8 4 0 7 135 0 0 0 1 1 136 0 0 1 0 0 137 0 2 6 0 0 138 2 3 3 5 1 139 5 0 0 2 0 140 7 5 5 6 4 141 4 0 3 4 0 142 2 2 0 2 6 143 3 7 0 11 0 144 0 0 3 0 1 145 0 2 0 4 7 146 2 0 0 2 1 147 0 2 3 1 7 148 4 2 0 3 0 149 2 2 1 0 2 150 9 2 3 4 4 151 0 2 1 4 2 152 2 3 2 5 2 153 5 0 4 4 0 154 0 1 0 3 0 155 5 0 9 0 5 156 6 3 2 0 3 157 0 4 4 12 0 158 3 2 0 2 0 159 2 3 7 4 0 160 4 0 0 4 0 161 0 0 0 0 0 162 0 10 4 10 5 163 4 3 0 2 1 164 3 1 2 0 2 165 0 0 0 0 0 166 1 0 1 0 0 167 0 2 0 0 0 168 4 0 3 1 2 169 5 3 0 3 3 170 1 1 1 3 1 171 5 0 0 12 7 172 3 0 0 0 7 173 1 2 0 0 1 174 0 2 5 7 3 175 5 2 4 0 3 176 0 3 0 7 0 177 1 6 5 0 0 178 3 4 0 3 3 179 7 4 9 0 5 180 0 0 0 0 0 181 4 4 4 0 2 182 0 1 1 1 0 183 0 8 6 10 7 184 0 3 4 4 7 185 0 1 0 0 0 186 4 3 0 4 0 187 0 4 6 0 0 188 0 0 0 1 0 189 0 4 4 0 0 190 0 3 0 0 4 191 0 5 0 0 6 192 7 6 2 3 5 193 9 9 10 9 19 194 1 0 8 5 0 195 7 3 5 1 2 196 0 0 5 1 6 197 4 4 7 0 0 198 2 2 3 0 1 199 0 0 0 3 2 200 4 6 9 2 3 svabu> m03 <- bootstrap(m00, B=5, type="param") svabu> extractBOOT(m03) $coefficients sta_(Intercept) sta_x1 sta_x5 det_(Intercept) det_x2 2.0934407 -0.9086257 0.4720519 0.6261731 1.7255167 det_x5 zif_(Intercept) -0.1601880 -0.7818176 $std.error sta_(Intercept) sta_x1 sta_x5 det_(Intercept) det_x2 0.18579853 0.04283707 0.18237568 0.72213938 0.33983140 det_x5 zif_(Intercept) 0.57490371 0.09053741 $vcov sta_(Intercept) sta_x1 sta_x5 det_(Intercept) sta_(Intercept) 0.034521093 0.005370583 -0.026364415 -0.12868296 sta_x1 0.005370583 0.001835014 -0.003838381 -0.02024922 sta_x5 -0.026364415 -0.003838381 0.033260887 0.08710034 det_(Intercept) -0.128682960 -0.020249222 0.087100345 0.52148529 det_x2 -0.037302112 -0.009820333 0.001512849 0.14643272 det_x5 0.080818915 0.007767725 -0.090503087 -0.32067756 zif_(Intercept) 0.009392050 0.003510249 -0.004462685 -0.04200064 det_x2 det_x5 zif_(Intercept) sta_(Intercept) -0.037302112 0.080818915 0.009392050 sta_x1 -0.009820333 0.007767725 0.003510249 sta_x5 0.001512849 -0.090503087 -0.004462685 det_(Intercept) 0.146432717 -0.320677563 -0.042000637 det_x2 0.115485378 0.007717497 -0.020563069 det_x5 0.007717497 0.330514274 0.014323815 zif_(Intercept) -0.020563069 0.014323815 0.008197022 $cor sta_(Intercept) sta_x1 sta_x5 det_(Intercept) sta_(Intercept) 1.0000000 0.6747756 -0.77805260 -0.9590864 sta_x1 0.6747756 1.0000000 -0.49131660 -0.6545873 sta_x5 -0.7780526 -0.4913166 1.00000000 0.6613509 det_(Intercept) -0.9590864 -0.6545873 0.66135087 1.0000000 det_x2 -0.5907825 -0.6745948 0.02440985 0.5966966 det_x5 0.7566162 0.3154126 -0.86318000 -0.7724181 zif_(Intercept) 0.5583289 0.9050867 -0.27027214 -0.6424019 det_x2 det_x5 zif_(Intercept) sta_(Intercept) -0.59078253 0.75661619 0.5583289 sta_x1 -0.67459477 0.31541256 0.9050867 sta_x5 0.02440985 -0.86318000 -0.2702721 det_(Intercept) 0.59669664 -0.77241809 -0.6424019 det_x2 1.00000000 0.03950189 -0.6683383 det_x5 0.03950189 1.00000000 0.2751918 zif_(Intercept) -0.66833835 0.27519183 1.0000000 $B [1] 5 $type [1] "param" $model [1] "full" svabu> summary(m03) Call: svabu(formula = Y ~ x1 + x5 | x2 + x5, data = databu[1:200, ]) Single visit Binomial - Zero Inflated Poisson model Conditional Maximum Likelihood estimates with Parametric Bootstrap standard errors (B = 5) Coefficients for abundance (log link): Estimate Std. Error z value Pr(>|z|) (Intercept) 2.03178 0.18580 10.935 < 2e-16 *** x1 -0.89211 0.04284 -20.826 < 2e-16 *** x5 0.53848 0.18238 2.953 0.00315 ** Coefficients for detection (logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) 0.8124 0.7221 1.125 0.261 x2 1.6873 0.3398 4.965 6.87e-07 *** x5 -0.4561 0.5749 -0.793 0.428 Coefficients for zero inflation (logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) -0.77581 0.09054 -8.569 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Log-likelihood: -300.1 on 7 Df AIC = 614.2 svabu> ## Model selection svabu> svabu> m04 <- svabu(Y ~ x1 + x5 | x2 + x5 + x3, databu[1:200,], phi.boot=0) svabu> m05 <- drop1(m04, model="det") svabu> m05 Single visit abundance model Single term deletion on the detection side Model: Y ~ x1 + x5 | x2 + x5 + x3 Dispersion parameter is taken to be 1 Df AIC Delta.AIC 620.15 0.000 x2 1 744.43 124.281 x5 1 617.89 -2.264 x3 1 614.22 -5.927 svabu> m06 <- svabu.step(m04, model="det") Start: AIC=620.15 Y ~ x1 + x5 | x2 + x5 + x3 Df AIC Delta.AIC - x3 1 614.22 -5.927 - x5 1 617.89 -2.264 620.15 0.000 - x2 1 744.43 124.281 Step: AIC=614.22 Y ~ x1 + x5 | x2 + x5 Df AIC Delta.AIC 614.22 0.0000 - x5 1 617.94 3.7165 - x2 svabu> summary(m06) Call: svabu(formula = Y ~ x1 + x5 | x2 + x5, data = databu[1:200, ], phi.boot = 0) Single visit Binomial - Zero Inflated Poisson model Conditional Maximum Likelihood estimates Coefficients for abundance (log link): Estimate Std. Error z value Pr(>|z|) (Intercept) 2.0318 0.1901 10.686 < 2e-16 *** x1 -0.8921 0.1701 -5.244 1.57e-07 *** x5 0.5385 0.1720 3.131 0.00174 ** Coefficients for detection (logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) 0.8124 0.6837 1.188 0.235 x2 1.6873 0.4040 4.177 2.96e-05 *** x5 -0.4561 0.5659 -0.806 0.420 Coefficients for zero inflation (logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) -0.7758 0.1705 -4.551 5.33e-06 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Log-likelihood: -300.1 on 7 Df AIC = 614.2 svabu> m07 <- update(m04, . ~ . | . - x3) svabu> m07 Call: svabu(formula = Y ~ x1 + x5 | x2 + x5, data = databu[1:200, ], phi.boot = 0) Single visit Binomial - Zero Inflated Poisson model Conditional Maximum Likelihood estimates Coefficients for abundance (log link): (Intercept) x1 x5 2.0318 -0.8921 0.5385 Coefficients for detection (logit link): (Intercept) x2 x5 0.8124 1.6873 -0.4561 Coefficients for zero inflation (logit link): (Intercept) -0.7758 svabu> ## Controls svabu> svabu> m00$control $optim.control $optim.control$maxit [1] 20000 $optim.control$fnscale [1] 1 $optim.method [1] "Nelder-Mead" svabu> getOption("detect.optim.control") $maxit [1] 20000 svabu> getOption("detect.optim.method") NULL svabu> options("detect.optim.method"="BFGS") svabu> m08 <- svabu(Y ~ x1 + x5 | x2 + x5, databu[1:100,]) svabu> m08$control ## but original optim method is retained during model selection and bootstrap $optim.control $optim.control$maxit [1] 20000 $optim.control$fnscale [1] 1 $optim.method [1] "BFGS" svabu> ## fitted models can be used to provide initial values svabu> options("detect.optim.method"="Nelder-Mead") svabu> m09 <- svabu(Y ~ x1 + x5 | x2 + x5, databu[1:100,], inits=coef(m08)) svabu> ## Ovenbirds dataset svabu> svabu> data(oven) svabu> ovenc <- oven svabu> ovenc[, c(4:8,10:11)][] <- lapply(ovenc[, c(4:8,10:11)], scale) svabu> moven <- svabu(count ~ pforest | observ + pforest + julian + timeday, ovenc) svabu> summary(moven) Call: svabu(formula = count ~ pforest | observ + pforest + julian + timeday, data = ovenc) Single visit Binomial - Zero Inflated Poisson model Conditional Maximum Likelihood estimates Coefficients for abundance (log link): Estimate Std. Error z value Pr(>|z|) (Intercept) 0.1683 0.2152 0.782 0.434 pforest 1.4556 0.1823 7.983 1.43e-15 *** Coefficients for detection (logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) -0.96016 0.80160 -1.198 0.23099 observDW 1.10912 1.45153 0.764 0.44480 observRDW 1.71846 0.68216 2.519 0.01176 * observSVW 2.21138 0.68903 3.209 0.00133 ** pforest -1.87652 0.38730 -4.845 1.27e-06 *** julian -0.31697 0.08069 -3.928 8.56e-05 *** timeday 0.01570 0.09101 0.173 0.86302 Coefficients for zero inflation (logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) -0.7846 0.1492 -5.259 1.45e-07 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Log-likelihood: -511.5 on 10 Df AIC = 1043 svabu> drop1(moven, model="det") Single visit abundance model Single term deletion on the detection side Model: count ~ pforest | observ + pforest + julian + timeday Dispersion parameter is taken to be 1 Df AIC Delta.AIC 1043.10 0.00 observ 3 518.61 -524.48 pforest 1 150.60 -892.49 julian 1 1046.18 3.08 timeday 1 949.87 -93.23 svabu> moven2 <- update(moven, . ~ . | . - timeday) svabu> summary(moven2) Call: svabu(formula = count ~ pforest | observ + pforest + julian, data = ovenc) Single visit Binomial - Zero Inflated Poisson model Conditional Maximum Likelihood estimates Coefficients for abundance (log link): Estimate Std. Error z value Pr(>|z|) (Intercept) 0.3220 0.2350 1.370 0.171 pforest 1.3870 0.1704 8.139 3.99e-16 *** Coefficients for detection (logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) -0.87617 0.76057 -1.152 0.24932 observDW -0.24569 1.47443 -0.167 0.86766 observRDW 1.34481 0.59672 2.254 0.02422 * observSVW 1.77861 0.59878 2.970 0.00297 ** pforest -1.70483 0.39992 -4.263 2.02e-05 *** julian -0.29677 0.07795 -3.807 0.00014 *** Coefficients for zero inflation (logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) -0.6765 0.1410 -4.798 1.6e-06 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Log-likelihood: -465.9 on 9 Df AIC = 949.9 svabu> moven3 <- update(moven2, . ~ . | ., zeroinfl=FALSE) svabu> summary(moven3) Call: svabu(formula = count ~ pforest | observ + pforest + julian, data = ovenc, zeroinfl = FALSE) Single visit Binomial - Poisson model Conditional Maximum Likelihood estimates Coefficients for abundance (log link): Estimate Std. Error z value Pr(>|z|) (Intercept) 1.8999 0.2288 8.305 < 2e-16 *** pforest -0.5142 0.1515 -3.394 0.000688 *** Coefficients for detection (logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) -3.79495 0.34949 -10.859 < 2e-16 *** observDW 0.58863 0.51755 1.137 0.255392 observRDW 0.82356 0.29249 2.816 0.004868 ** observSVW 1.08290 0.30262 3.578 0.000346 *** pforest 1.75392 0.25170 6.968 3.21e-12 *** julian -0.33340 0.06735 -4.951 7.40e-07 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Log-likelihood: -730.6 on 8 Df AIC = 1477 svabu> BIC(moven, moven2, moven3) df BIC moven 10 1091.0196 moven2 9 992.9991 moven3 8 1515.6219 ---------- detect example: svocc ---------- svocc> data(datocc) svocc> ## MLE svocc> m00 <- svocc(W ~ x1 | x1 + x3, datocc) svocc> ## PMLE svocc> m01 <- svocc(W ~ x1 | x1 + x3, datocc, penalized=TRUE) svocc> ## print svocc> m00 Call: svocc(formula = W ~ x1 | x1 + x3, data = datocc) Single visit site-occupancy model Maximum Likelihood estimates (optim method) Coefficients for occurrence (cloglog link): (Intercept) x1 0.2216 0.4496 Coefficients for detection (logit link): (Intercept) x1 x3 0.8216 -0.7650 0.3741 svocc> ## summary svocc> summary(m00) Call: svocc(formula = W ~ x1 | x1 + x3, data = datocc) Single visit site-occupancy model Maximum Likelihood estimates (optim method) Occupancy model coefficients with cloglog link: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.2216 0.6695 0.331 0.741 x1 0.4496 0.3896 1.154 0.248 Detection model coefficients with logit link: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.8216 1.1210 0.733 0.464 x1 -0.7650 0.5540 -1.381 0.167 x3 0.3741 0.2848 1.313 0.189 Log-likelihood: -685.6 on 5 Df AIC = 1381 svocc> ## coefficients svocc> coef(m00) sta_(Intercept) sta_x1 det_(Intercept) det_x1 det_x3 0.2215515 0.4496429 0.8216246 -0.7650274 0.3740887 svocc> ## state (occupancy) model estimates svocc> coef(m00, "sta") (Intercept) x1 0.2215515 0.4496429 svocc> ## detection model estimates svocc> coef(m00, "det") (Intercept) x1 x3 0.8216246 -0.7650274 0.3740887 svocc> ## compare estimates svocc> cbind(truth=c(0.6, 0.5, 0.4, -0.5, 0.3), svocc+ mle=coef(m00), pmle=coef(m01)) truth mle pmle sta_(Intercept) 0.6 0.2215515 0.04036723 sta_x1 0.5 0.4496429 0.24861331 det_(Intercept) 0.4 0.8216246 1.12888938 det_x1 -0.5 -0.7650274 -0.61112990 det_x3 0.3 0.3740887 0.47112225 svocc> ## AIC, BIC svocc> AIC(m00) [1] 1381.247 svocc> BIC(m00) [1] 1405.786 svocc> ## log-likelihood svocc> logLik(m00) 'log Lik.' -685.6237 (df=5) svocc> ## variance-covariance matrix svocc> vcov(m00) sta_(Intercept) sta_x1 det_(Intercept) det_x1 sta_(Intercept) 0.4482888 0.21233677 -0.7442351 0.26789589 sta_x1 0.2123368 0.15181417 -0.3340672 0.05230045 det_(Intercept) -0.7442351 -0.33406715 1.2567311 -0.48353333 det_x1 0.2678959 0.05230045 -0.4835333 0.30688349 det_x3 -0.1754498 -0.08306497 0.2943395 -0.10712543 det_x3 sta_(Intercept) -0.17544983 sta_x1 -0.08306497 det_(Intercept) 0.29433945 det_x1 -0.10712543 det_x3 0.08111848 svocc> vcov(m00, model="sta") (Intercept) x1 (Intercept) 0.4482888 0.2123368 x1 0.2123368 0.1518142 svocc> vcov(m00, model="det") (Intercept) x1 x3 (Intercept) 1.2567311 -0.4835333 0.29433945 x1 -0.4835333 0.3068835 -0.10712543 x3 0.2943395 -0.1071254 0.08111848 svocc> ## confidence intervals svocc> confint(m00) 2.5% 97.5% sta_(Intercept) -1.0907301 1.5338330 sta_x1 -0.3140244 1.2133103 det_(Intercept) -1.3755738 3.0188230 det_x1 -1.8507899 0.3207351 det_x3 -0.1841346 0.9323121 svocc> confint(m00, model="sta") 2.5% 97.5% (Intercept) -1.0907301 1.533833 x1 -0.3140244 1.213310 svocc> confint(m00, model="det") 2.5% 97.5% (Intercept) -1.3755738 3.0188230 x1 -1.8507899 0.3207351 x3 -0.1841346 0.9323121 svocc> ## fitted values svocc> ## (conditional probability of occurrence given detection history: svocc> ## if W=1, fitted=1, svocc> ## if W=0, fitted=(phi*(1-delta)) / ((1-delta) + phi * (1-delta)) svocc> summary(fitted(m00)) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.09244 0.42641 0.74961 0.71202 1.00000 1.00000 svocc> ## estimated probabilities: (phi*(1-delta)) / ((1-delta) + phi * (1-delta)) svocc> summary(m00$estimated.probabilities) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.09244 0.28628 0.43414 0.44293 0.59901 0.80221 svocc> ## probability of occurrence (phi) svocc> summary(m00$occurrence.probabilities) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.5490 0.6336 0.7162 0.7121 0.7929 0.8585 svocc> ## probability of detection (delta) svocc> summary(m00$detection.probabilities) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.2660 0.5961 0.6966 0.6811 0.7778 0.9290 svocc> ## model selection svocc> m02 <- svocc(W ~ x1 | x3 + x4, datocc) svocc> m03 <- drop1(m02, model="det") svocc> ## dropping one term at a time, resulting change in AIC svocc> m03 Single visit site-occupancy model Single term deletion on the detection side Model: W ~ x1 | x3 + x4 Dispersion parameter is taken to be 1 Df AIC Delta.AIC 1381.5 0.0000 x3 1 1392.1 10.6055 x4 1 1379.5 -1.9972 svocc> ## updating the model svocc> m04 <- update(m02, . ~ . | . - x4) svocc> m04 Call: svocc(formula = W ~ x1 | x3, data = datocc) Single visit site-occupancy model Maximum Likelihood estimates (optim method) Coefficients for occurrence (cloglog link): (Intercept) x1 -0.24163 -0.01527 Coefficients for detection (logit link): (Intercept) x3 2.2604 0.9657 svocc> ## automatic model selection svocc> ## part of the model (sta/det) must be specified svocc> m05 <- svocc.step(m02, model="det") Start: AIC=1381.49 W ~ x1 | x3 + x4 Df AIC Delta.AIC - x4 1 1379.5 -1.9972 1381.5 0.0000 - x3 1 1392.1 10.6055 Step: AIC=1379.5 W ~ x1 | x3 Df AIC Delta.AIC 1379.5 0 - x3 svocc> summary(m05) Call: svocc(formula = W ~ x1 | x3, data = datocc) Single visit site-occupancy model Maximum Likelihood estimates (optim method) Occupancy model coefficients with cloglog link: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.24163 0.13843 -1.746 0.0809 . x1 -0.01527 0.08583 -0.178 0.8588 Detection model coefficients with logit link: Estimate Std. Error z value Pr(>|z|) (Intercept) 2.2604 1.0523 2.148 0.0317 * x3 0.9657 0.5239 1.843 0.0653 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Log-likelihood: -685.7 on 4 Df AIC = 1379 svocc> ## nonparametric bootstrap svocc> m06 <- bootstrap(m01, B=25) svocc> attr(m06, "bootstrap") cfs sta_(Intercept) 0.04036723 0.2402178 0.7753206 -0.1792962 0.5373761 sta_x1 0.24861331 0.3752813 0.5004229 -0.3248569 1.0280805 det_(Intercept) 1.12888938 0.6532832 0.1655229 1.6509545 0.5609020 det_x1 -0.61112990 -0.4731551 -0.3056861 -0.4016467 -1.0049216 det_x3 0.47112225 0.4259106 0.2704690 1.0503999 0.3881481 sta_(Intercept) 0.3943155 0.61255138 0.7965053 -0.2379533 0.4036633 sta_x1 0.5861091 0.11615250 0.9542351 0.1992000 0.5172542 det_(Intercept) 0.6167747 0.09711886 0.5122373 2.4863923 0.5516786 det_x1 -1.0711661 0.10256412 -0.6153502 -0.5379005 -0.5032798 det_x3 0.2617254 0.14751757 0.2149977 1.4004248 0.2975085 sta_(Intercept) 0.5566774 0.0233945 0.68889332 -0.02553195 0.08126412 sta_x1 0.9830122 0.3231839 -0.05293293 0.22043090 0.41451040 det_(Intercept) 0.5609432 1.6650879 0.12619883 1.41141437 0.98669175 det_x1 -1.2297401 -1.1012976 -0.27290788 -0.45300722 -0.65655924 det_x3 0.3596743 0.7862023 0.19617040 0.80277609 0.30632249 sta_(Intercept) -0.1224226 0.01483467 0.3237822 0.77616868 0.02992533 sta_x1 0.3204974 0.21184657 0.6788141 0.05378419 0.14895352 det_(Intercept) 1.7541092 1.55730461 0.5773976 0.17117805 1.34650883 det_x1 -1.0315339 -1.07778121 -0.7705330 0.07426017 -0.84876741 det_x3 0.8895609 0.41821723 0.4185265 0.18902731 0.74165453 sta_(Intercept) -0.1767476 -0.24604128 -0.02930077 -0.15731837 0.72768951 sta_x1 -0.1280581 0.07074902 0.61840882 -0.09904495 0.09736739 det_(Intercept) 2.4741503 2.32745538 1.53832192 2.09014609 0.22686429 det_x1 -0.2014874 -0.95243720 -1.35205421 -0.16633287 0.01296555 det_x3 1.3320994 1.10931429 0.68754939 0.94951427 0.21394511 sta_(Intercept) 0.07321718 sta_x1 0.46657749 det_(Intercept) 1.53674659 det_x1 -1.45668175 det_x3 0.49290294 attr(,"type") [1] "nonpar" attr(,"ini") sta_(Intercept) sta_x1 det_(Intercept) det_x1 det_x3 0.04036723 0.24861331 1.12888938 -0.61112990 0.47112225 svocc> extractBOOT(m06) $coefficients sta_(Intercept) sta_x1 det_(Intercept) det_x1 det_x3 0.2277520 0.3280228 1.1067028 -0.6502141 0.5700647 $std.error sta_(Intercept) sta_x1 det_(Intercept) det_x1 det_x3 0.3569504 0.3441538 0.7614797 0.4407125 0.3680701 $vcov sta_(Intercept) sta_x1 det_(Intercept) det_x1 sta_(Intercept) 0.12741361 0.05056662 -0.25323548 0.05094126 sta_x1 0.05056662 0.11844182 -0.10077588 -0.08306637 det_(Intercept) -0.25323548 -0.10077588 0.57985135 -0.09737934 det_x1 0.05094126 -0.08306637 -0.09737934 0.19422747 det_x3 -0.11195893 -0.05654749 0.26067137 -0.01631319 det_x3 sta_(Intercept) -0.11195893 sta_x1 -0.05654749 det_(Intercept) 0.26067137 det_x1 -0.01631319 det_x3 0.13547558 $cor sta_(Intercept) sta_x1 det_(Intercept) det_x1 sta_(Intercept) 1.0000000 0.4116267 -0.9316618 0.3238221 sta_x1 0.4116267 1.0000000 -0.3845437 -0.5476680 det_(Intercept) -0.9316618 -0.3845437 1.0000000 -0.2901704 det_x1 0.3238221 -0.5476680 -0.2901704 1.0000000 det_x3 -0.8521584 -0.4464063 0.9300462 -0.1005664 det_x3 sta_(Intercept) -0.8521584 sta_x1 -0.4464063 det_(Intercept) 0.9300462 det_x1 -0.1005664 det_x3 1.0000000 $B [1] 25 $type [1] "nonpar" svocc> summary(m06, type="mle") Call: svocc(formula = W ~ x1 | x1 + x3, data = datocc, penalized = TRUE) Single visit site-occupancy model Maximum Likelihood estimates (optim method) Occupancy model coefficients with cloglog link: Estimate Std. Error z value Pr(>|z|) sta_(Intercept) 0.2216 0.6695 0.331 0.741 sta_x1 0.4496 0.3896 1.154 0.248 Detection model coefficients with logit link: Estimate Std. Error z value Pr(>|z|) det_(Intercept) 0.8216 1.1210 0.733 0.464 det_x1 -0.7650 0.5540 -1.381 0.167 det_x3 0.3741 0.2848 1.313 0.189 Log-likelihood: -686.7 on 5 Df AIC = 1383 svocc> summary(m06, type="pmle") ## no SEs! PMLE!!! Call: svocc(formula = W ~ x1 | x1 + x3, data = datocc, penalized = TRUE) Single visit site-occupancy model Penalized Maximum Likelihood estimates (optim method) Occupancy model coefficients with cloglog link: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.04037 NA NA NA x1 0.24861 NA NA NA Detection model coefficients with logit link: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.1289 NA NA NA x1 -0.6111 NA NA NA x3 0.4711 NA NA NA Log-likelihood: -686.7 on 5 Df AIC = 1383 svocc> summary(m06, type="boot") Call: svocc(formula = W ~ x1 | x1 + x3, data = datocc, penalized = TRUE) Single visit site-occupancy model Maximum Likelihood estimates (optim method) with Nonparametric Bootstrap standard errors (B = 25) Occupancy model coefficients with cloglog link: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.04037 0.35695 0.113 0.91 x1 0.24861 0.34415 0.722 0.47 Detection model coefficients with logit link: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.1289 0.7615 1.482 0.138 x1 -0.6111 0.4407 -1.387 0.166 x3 0.4711 0.3681 1.280 0.201 Log-likelihood: -686.7 on 5 Df AIC = 1383 svocc> ## vcov svocc> #vcov(m06, type="mle") ## this does not work with PMLE svocc> vcov(m06, type="boot") ## this works sta_(Intercept) sta_x1 det_(Intercept) det_x1 sta_(Intercept) 0.12741361 0.05056662 -0.25323548 0.05094126 sta_x1 0.05056662 0.11844182 -0.10077588 -0.08306637 det_(Intercept) -0.25323548 -0.10077588 0.57985135 -0.09737934 det_x1 0.05094126 -0.08306637 -0.09737934 0.19422747 det_x3 -0.11195893 -0.05654749 0.26067137 -0.01631319 det_x3 sta_(Intercept) -0.11195893 sta_x1 -0.05654749 det_(Intercept) 0.26067137 det_x1 -0.01631319 det_x3 0.13547558 svocc> ## confint svocc> confint(m06, type="boot") ## quantile based 2.5% 97.5% sta_(Intercept) -0.2409863 0.78379492 sta_x1 -0.2018576 0.99991282 det_(Intercept) 0.1152938 2.47874102 det_x1 -1.3912895 0.08487415 det_x3 0.1734612 1.35772142 svocc> ## parametric bootstrap svocc> ## sthis is how observations are simulated svocc> head(simulate(m01, 5)) sim_1 sim_2 sim_3 sim_4 sim_5 1 1 1 0 1 0 2 1 1 1 1 1 3 1 0 0 1 0 4 1 0 1 1 0 5 1 0 1 1 1 6 0 0 0 0 1 svocc> m07 <- bootstrap(m01, B=25, type="param") svocc> extractBOOT(m07) $coefficients sta_(Intercept) sta_x1 det_(Intercept) det_x1 det_x3 0.4076715 0.2235850 0.7051518 -0.2149875 0.4232817 $std.error sta_(Intercept) sta_x1 det_(Intercept) det_x1 det_x3 0.4624367 0.3698566 0.5688444 0.4142387 0.1994603 $vcov sta_(Intercept) sta_x1 det_(Intercept) det_x1 sta_(Intercept) 0.21384769 0.07901743 -0.23524337 -0.02588020 sta_x1 0.07901743 0.13679393 -0.07956099 -0.09299702 det_(Intercept) -0.23524337 -0.07956099 0.32358390 0.04731344 det_x1 -0.02588020 -0.09299702 0.04731344 0.17159372 det_x3 -0.06810256 -0.02283834 0.10268871 0.02098074 det_x3 sta_(Intercept) -0.06810256 sta_x1 -0.02283834 det_(Intercept) 0.10268871 det_x1 0.02098074 det_x3 0.03978441 $cor sta_(Intercept) sta_x1 det_(Intercept) det_x1 sta_(Intercept) 1.0000000 0.4619949 -0.8942762 -0.1351029 sta_x1 0.4619949 1.0000000 -0.3781580 -0.6069947 det_(Intercept) -0.8942762 -0.3781580 1.0000000 0.2007892 det_x1 -0.1351029 -0.6069947 0.2007892 1.0000000 det_x3 -0.7383371 -0.3095813 0.9050504 0.2539298 det_x3 sta_(Intercept) -0.7383371 sta_x1 -0.3095813 det_(Intercept) 0.9050504 det_x1 0.2539298 det_x3 1.0000000 $B [1] 25 $type [1] "param" svocc> summary(m07) Call: svocc(formula = W ~ x1 | x1 + x3, data = datocc, penalized = TRUE) Single visit site-occupancy model Maximum Likelihood estimates (optim method) with Parametric Bootstrap standard errors (B = 25) Occupancy model coefficients with cloglog link: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.04037 0.46244 0.087 0.930 x1 0.24861 0.36986 0.672 0.501 Detection model coefficients with logit link: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.1289 0.5688 1.985 0.0472 * x1 -0.6111 0.4142 -1.475 0.1401 x3 0.4711 0.1995 2.362 0.0182 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Log-likelihood: -686.7 on 5 Df AIC = 1383 svocc> data(oven) svocc> ovenc <- oven svocc> ovenc[, c(4:8,10:11)][] <- lapply(ovenc[, c(4:8,10:11)], scale) svocc> ovenc$count01 <- ifelse(ovenc$count > 0, 1, 0) svocc> moven <- svocc(count01 ~ pforest | julian + timeday, ovenc) svocc> summary(moven) Call: svocc(formula = count01 ~ pforest | julian + timeday, data = ovenc) Single visit site-occupancy model Maximum Likelihood estimates (optim method) Occupancy model coefficients with cloglog link: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.2734 0.3047 0.897 0.369 pforest 2.6014 0.4425 5.879 4.13e-09 *** Detection model coefficients with logit link: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.24243 0.11921 2.034 0.0420 * julian -0.21494 0.09397 -2.287 0.0222 * timeday 0.08500 0.09678 0.878 0.3798 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Log-likelihood: -420.5 on 5 Df AIC = 851 svocc> drop1(moven, model="det") Single visit site-occupancy model Single term deletion on the detection side Model: count01 ~ pforest | julian + timeday Dispersion parameter is taken to be 1 Df AIC Delta.AIC 850.95 0.0000 julian 1 854.37 3.4225 timeday 1 849.73 -1.2192 svocc> moven2 <- update(moven, . ~ . | . - timeday) svocc> summary(moven) Call: svocc(formula = count01 ~ pforest | julian + timeday, data = ovenc) Single visit site-occupancy model Maximum Likelihood estimates (optim method) Occupancy model coefficients with cloglog link: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.2734 0.3047 0.897 0.369 pforest 2.6014 0.4425 5.879 4.13e-09 *** Detection model coefficients with logit link: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.24243 0.11921 2.034 0.0420 * julian -0.21494 0.09397 -2.287 0.0222 * timeday 0.08500 0.09678 0.878 0.3798 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Log-likelihood: -420.5 on 5 Df AIC = 851 svocc> BIC(moven, moven2) df BIC moven 5 874.9123 moven2 4 868.9007 svocc> AUC(moven, moven2) AUC moven 0.8076323 moven2 0.8041871 svocc> rocplot(moven) svocc> rocplot(moven2, col=2, add=TRUE) Warning messages: 1: In svabu.fit(y, xother, x[, jj, drop = FALSE], N.max = object$N.max, : negative variance values in optim, NAs produced 2: In svabu.fit(y, xother, x[, jj, drop = FALSE], N.max = object$N.max, : negative variance values in optim, NAs produced 3: In svocc.fit(y, xother, x[, jj, drop = FALSE], link.sta = object$link$sta, : negative variance values in optim, NAs produced 4: In svocc.fit(y, xother, x[, jj, drop = FALSE], link.sta = object$link$sta, : negative variance values in optim, NAs produced > > if (FALSE) { + + library(detect) + data(datocc) + m0 <- svocc(W ~ x1 | x1 + x3, datocc, + method="dc",n.adapt=100,n.update=100,n.iter=100) + extractMLE(m0) + + } > > proc.time() user system elapsed 345.26 4.90 350.15