R Under development (unstable) (2025-07-11 r88405 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/opticut") > library(opticut) Loading required package: pbapply opticut 0.1-4 2025-07-12 > > ## --- testing methods --- > > ocoptions(try_error=TRUE) > set.seed(2345) > n <- 50 > x0 <- sample(1:4, n, TRUE) > x1 <- ifelse(x0 %in% 1:2, 1, 0) > x2 <- rnorm(n, 0.5, 1) > x3 <- ifelse(x0 %in% 2:4, 1, 0) > mu1 <- 0.5 + 1*x1 + -0.2*x2 > mu2 <- 1 + 0.5*x3 > mu3 <- rep(0, n) > Y1 <- rpois(n, exp(mu1)) > Y2 <- rpois(n, exp(mu2)) > Y3 <- rpois(n, exp(mu3)) > Y <- cbind(Spp1=Y1, Spp2=Y2, Spp3=Y3) > comb <- allComb(x0)[,1:6] > attr(comb, "collapse") <- NULL > attr(comb, "comb") <- NULL > colnames(comb) <- letters[1:6] > m1 <- opticut(Y ~ x2, strata=x0, dist="poisson", comb="rank") > m2 <- opticut(Y ~ x2, strata=x0, dist="poisson", comb="all") > m3 <- opticut(Y ~ x2, strata=comb, dist="poisson") > fun <- function(Y, X, linkinv, ...) { + mod <- stats::glm(Y ~ .-1, data=X, family="poisson", ...) + list(coef=coef(mod), + logLik=logLik(mod), + linkinv=family(mod)$linkinv) + } > m4 <- opticut(Y ~ x2, strata=x0, dist=fun) > ocoptions(try_error=FALSE) > > subset(m1, c(3,1)) Multivariate opticut results, comb = rank, dist = poisson Call: opticut.formula(formula = Y ~ x2, strata = x0, dist = "poisson", comb = "rank") 2 species, 3 binary splits > subset(m2, c(TRUE, FALSE, TRUE)) Multivariate opticut results, comb = all, dist = poisson Call: opticut.formula(formula = Y ~ x2, strata = x0, dist = "poisson", comb = "all") 2 species, 7 binary splits > subset(m3, c("Spp1", "Spp3")) Multivariate opticut results, comb = NA, dist = poisson Call: opticut.formula(formula = Y ~ x2, strata = comb, dist = "poisson") 2 species, 6 binary splits > > str(m1$strata) Factor w/ 4 levels "1","2","3","4": 3 2 3 3 2 3 4 3 3 2 ... > str(m2$strata) int [1:50, 1:7] 0 0 0 0 0 0 0 0 0 0 ... - attr(*, "dimnames")=List of 2 ..$ : chr [1:50] "3" "2" "3" "3" ... ..$ : chr [1:7] "1" "2" "3" "4" ... - attr(*, "collapse")= chr "+" - attr(*, "comb")= chr "all" > str(m3$strata) int [1:50, 1:6] 0 0 0 0 0 0 0 0 0 0 ... - attr(*, "dimnames")=List of 2 ..$ : chr [1:50] "3" "2" "3" "3" ... ..$ : chr [1:6] "a" "b" "c" "d" ... > str(m4$strata) Factor w/ 4 levels "1","2","3","4": 3 2 3 3 2 3 4 3 3 2 ... > > opticut:::.extractOpticut(m1) $Spp1 $Spp1$coef `(Intercept)` Z1 x2 0.6774482 0.7934881 -0.1945135 $Spp1$logLik [1] -90.67904 $Spp1$linkinv NULL $Spp2 $Spp2$coef `(Intercept)` Z1 x2 1.21288873 0.26682754 -0.03234522 $Spp2$logLik [1] -98.92923 $Spp2$linkinv NULL $Spp3 $Spp3$coef `(Intercept)` Z1 x2 -1.04471822 0.71031616 0.07022314 $Spp3$logLik [1] -49.74613 $Spp3$linkinv NULL > opticut:::.extractOpticut(m2) $Spp1 $Spp1$coef `(Intercept)` Z1 x2 0.6774482 0.7934881 -0.1945135 $Spp1$logLik [1] -90.67904 $Spp1$linkinv NULL $Spp2 $Spp2$coef `(Intercept)` Z1 x2 1.47971628 -0.26682754 -0.03234522 $Spp2$logLik [1] -98.92923 $Spp2$linkinv NULL $Spp3 $Spp3$coef `(Intercept)` Z1 x2 -0.33440206 -0.71031616 0.07022314 $Spp3$logLik [1] -49.74613 $Spp3$linkinv NULL > opticut:::.extractOpticut(m3) $Spp1 $Spp1$coef `(Intercept)` Z1 x2 0.6774482 0.7934881 -0.1945135 $Spp1$logLik [1] -90.67904 $Spp1$linkinv NULL $Spp2 $Spp2$coef `(Intercept)` Z1 x2 1.47971628 -0.26682754 -0.03234522 $Spp2$logLik [1] -98.92923 $Spp2$linkinv NULL $Spp3 $Spp3$coef `(Intercept)` Z1 x2 -0.33440206 -0.71031616 0.07022314 $Spp3$logLik [1] -49.74613 $Spp3$linkinv NULL > opticut:::.extractOpticut(m4) $Spp1 $Spp1$coef `(Intercept)` Z1 x2 0.6774482 0.7934881 -0.1945135 $Spp1$logLik 'log Lik.' -90.67904 (df=3) $Spp1$linkinv function (eta) pmax(exp(eta), .Machine$double.eps) $Spp2 $Spp2$coef `(Intercept)` Z1 x2 1.21288873 0.26682754 -0.03234522 $Spp2$logLik 'log Lik.' -98.92923 (df=3) $Spp2$linkinv function (eta) pmax(exp(eta), .Machine$double.eps) $Spp3 $Spp3$coef `(Intercept)` Z1 x2 -1.04471822 0.71031616 0.07022314 $Spp3$logLik 'log Lik.' -49.74613 (df=3) $Spp3$linkinv function (eta) pmax(exp(eta), .Machine$double.eps) > > u1a <- uncertainty(m1, type="asymp", B=99) > u2a <- uncertainty(m2, type="asymp", B=99) > u3a <- uncertainty(m3, type="asymp", B=99) > ## asymp needs Hessian: dist=fun cannot provide that > u4a <- try(uncertainty(m4, type="asymp", B=999), silent=TRUE) > stopifnot(inherits(u4a, "try-error")) > > u1b <- uncertainty(m1, type="boot", B=9) > u2b <- uncertainty(m2, type="boot", B=9) > u3b <- uncertainty(m3, type="boot", B=9) > u4b <- uncertainty(m4, type="boot", B=9) > > summary(subset(u1b, c(3,1))) Multivariate multicut uncertainty results type = boot, B = 9, level = 0.95 split R I Lower Upper Spp1 1+2 1 0.3967 0.3139 0.4764 Spp3 1+2+3 1 0.4992 0.2459 0.8167 > summary(subset(u2b, c(TRUE, FALSE, TRUE))) Multivariate multicut uncertainty results type = boot, B = 9, level = 0.95 split R I Lower Upper Spp1 1+2 1 0.3754 0.1910 0.4945 Spp3 4 1 0.3541 0.1082 0.6908 > summary(subset(u3b, c("Spp1", "Spp3"))) Multivariate multicut uncertainty results type = boot, B = 9, level = 0.95 split R I Lower Upper Spp3 d 1 0.3776 0.2010 0.6349 Spp1 e 1 0.3350 0.2377 0.4593 > > u1c <- uncertainty(m1, type="multi", B=9) > ## type=multi cannot use object with comb=all > u2c <- try(uncertainty(m2, type="multi", B=9), silent=TRUE) > stopifnot(inherits(u2c, "try-error")) > ## type=multi cannot use object with comb=NA (custom partitions) > u3c <- try(uncertainty(m3, type="multi", B=9), silent=TRUE) > stopifnot(inherits(u3c, "try-error")) > u4c <- uncertainty(m4, type="multi", B=9) > > strata(m1) [1] 3 2 3 3 2 3 4 3 3 2 1 4 1 2 2 1 1 1 4 1 3 4 4 1 3 2 4 3 4 4 4 1 3 2 3 2 4 3 [39] 4 3 3 2 3 1 1 2 2 2 3 2 Levels: 1 2 3 4 > strata(m2) [1] 3 2 3 3 2 3 4 3 3 2 1 4 1 2 2 1 1 1 4 1 3 4 4 1 3 2 4 3 4 4 4 1 3 2 3 2 4 3 [39] 4 3 3 2 3 1 1 2 2 2 3 2 Levels: 1 2 3 4 > strata(m3) 3 2 3 3 2 3 4 3 3 2 1 001001 010010 001001 001001 010010 001001 000100 001001 001001 010010 100011 4 1 2 2 1 1 1 4 1 3 4 000100 100011 010010 010010 100011 100011 100011 000100 100011 001001 000100 4 1 3 2 4 3 4 4 4 1 3 000100 100011 001001 010010 000100 001001 000100 000100 000100 100011 001001 2 3 2 4 3 4 3 3 2 3 1 010010 001001 010010 000100 001001 000100 001001 001001 010010 001001 100011 1 2 2 2 3 2 100011 010010 010010 010010 001001 010010 Levels: 000100 001001 010010 100011 > strata(m4) [1] 3 2 3 3 2 3 4 3 3 2 1 4 1 2 2 1 1 1 4 1 3 4 4 1 3 2 4 3 4 4 4 1 3 2 3 2 4 3 [39] 4 3 3 2 3 1 1 2 2 2 3 2 Levels: 1 2 3 4 > strata(u1b) > strata(u2b) > strata(u3b) > strata(u4b) > > bestmodel(m1) $Spp1 Call: stats::glm(formula = Y ~ . - 1, family = Family, data = XX) Coefficients: `(Intercept)` Z1 x2 0.6774 0.7935 -0.1945 Degrees of Freedom: 50 Total (i.e. Null); 47 Residual Null Deviance: 181.3 Residual Deviance: 53.76 AIC: 187.4 $Spp2 Call: stats::glm(formula = Y ~ . - 1, family = Family, data = XX) Coefficients: `(Intercept)` Z1 x2 1.21289 0.26683 -0.03235 Degrees of Freedom: 50 Total (i.e. Null); 47 Residual Null Deviance: 311.9 Residual Deviance: 38.26 AIC: 203.9 $Spp3 Call: stats::glm(formula = Y ~ . - 1, family = Family, data = XX) Coefficients: `(Intercept)` Z1 x2 -1.04472 0.71032 0.07022 Degrees of Freedom: 50 Total (i.e. Null); 47 Residual Null Deviance: 50.64 Residual Deviance: 41.81 AIC: 105.5 > bestmodel(m2) $Spp1 Call: stats::glm(formula = Y ~ . - 1, family = Family, data = XX) Coefficients: `(Intercept)` Z1 x2 0.6774 0.7935 -0.1945 Degrees of Freedom: 50 Total (i.e. Null); 47 Residual Null Deviance: 181.3 Residual Deviance: 53.76 AIC: 187.4 $Spp2 Call: stats::glm(formula = Y ~ . - 1, family = Family, data = XX) Coefficients: `(Intercept)` Z1 x2 1.47972 -0.26683 -0.03235 Degrees of Freedom: 50 Total (i.e. Null); 47 Residual Null Deviance: 311.9 Residual Deviance: 38.26 AIC: 203.9 $Spp3 Call: stats::glm(formula = Y ~ . - 1, family = Family, data = XX) Coefficients: `(Intercept)` Z1 x2 -0.33440 -0.71032 0.07022 Degrees of Freedom: 50 Total (i.e. Null); 47 Residual Null Deviance: 50.64 Residual Deviance: 41.81 AIC: 105.5 > bestmodel(m3) $Spp1 Call: stats::glm(formula = Y ~ . - 1, family = Family, data = XX) Coefficients: `(Intercept)` Z1 x2 0.6774 0.7935 -0.1945 Degrees of Freedom: 50 Total (i.e. Null); 47 Residual Null Deviance: 181.3 Residual Deviance: 53.76 AIC: 187.4 $Spp2 Call: stats::glm(formula = Y ~ . - 1, family = Family, data = XX) Coefficients: `(Intercept)` Z1 x2 1.47972 -0.26683 -0.03235 Degrees of Freedom: 50 Total (i.e. Null); 47 Residual Null Deviance: 311.9 Residual Deviance: 38.26 AIC: 203.9 $Spp3 Call: stats::glm(formula = Y ~ . - 1, family = Family, data = XX) Coefficients: `(Intercept)` Z1 x2 -0.33440 -0.71032 0.07022 Degrees of Freedom: 50 Total (i.e. Null); 47 Residual Null Deviance: 50.64 Residual Deviance: 41.81 AIC: 105.5 > ## dist=fun cannot return the best model (--> uncertainty(type=asymm) fails) > bm4 <- try(bestmodel(m4), silent=TRUE) # dist=fun problem > stopifnot(inherits(bm4, "try-error")) > > getMLE(m1, 1, vcov=FALSE) $coef `(Intercept)` Z1 x2 0.6774482 0.7934881 -0.1945135 $vcov NULL $dist [1] "poisson" > getMLE(m2, 1, vcov=FALSE) $coef `(Intercept)` Z1 x2 0.6774482 0.7934881 -0.1945135 $vcov NULL $dist [1] "poisson" > getMLE(m3, 1, vcov=FALSE) $coef `(Intercept)` Z1 x2 0.6774482 0.7934881 -0.1945135 $vcov NULL $dist [1] "poisson" > getMLE(m4, 1, vcov=FALSE) $coef `(Intercept)` Z1 x2 0.6774482 0.7934881 -0.1945135 $vcov NULL $dist [1] "fun" > getMLE(m1, 1, vcov=TRUE) $coef `(Intercept)` Z1 x2 0.6774482 0.7934881 -0.1945135 $vcov `(Intercept)` Z1 x2 `(Intercept)` 0.020846432 -0.019077412 -0.003239102 Z1 -0.019077412 0.032241518 -0.003530473 x2 -0.003239102 -0.003530473 0.012395193 $dist [1] "poisson" > getMLE(m2, 1, vcov=TRUE) $coef `(Intercept)` Z1 x2 0.6774482 0.7934881 -0.1945135 $vcov `(Intercept)` Z1 x2 `(Intercept)` 0.020846432 -0.019077412 -0.003239102 Z1 -0.019077412 0.032241518 -0.003530473 x2 -0.003239102 -0.003530473 0.012395193 $dist [1] "poisson" > getMLE(m3, 1, vcov=TRUE) $coef `(Intercept)` Z1 x2 0.6774482 0.7934881 -0.1945135 $vcov `(Intercept)` Z1 x2 `(Intercept)` 0.020846432 -0.019077412 -0.003239102 Z1 -0.019077412 0.032241518 -0.003530473 x2 -0.003239102 -0.003530473 0.012395193 $dist [1] "poisson" > mle4 <- try(getMLE(m4, 1, vcov=TRUE), silent=TRUE) > stopifnot(inherits(mle4, "try-error")) > > str(bestpart(m1)) int [1:50, 1:3] 0 1 0 0 1 0 0 0 0 1 ... - attr(*, "dimnames")=List of 2 ..$ : chr [1:50] "3" "2" "3" "3" ... ..$ : chr [1:3] "Spp1" "Spp2" "Spp3" > str(bestpart(m2)) int [1:50, 1:3] 0 1 0 0 1 0 0 0 0 1 ... - attr(*, "dimnames")=List of 2 ..$ : chr [1:50] "3" "2" "3" "3" ... ..$ : chr [1:3] "Spp1" "Spp2" "Spp3" > str(bestpart(m2, pos_only=TRUE)) num [1:50, 1:3] 0 1 0 0 1 0 0 0 0 1 ... - attr(*, "dimnames")=List of 2 ..$ : chr [1:50] "3" "2" "3" "3" ... ..$ : chr [1:3] "Spp1" "Spp2" "Spp3" > str(bestpart(m3)) int [1:50, 1:3] 0 1 0 0 1 0 0 0 0 1 ... - attr(*, "dimnames")=List of 2 ..$ : chr [1:50] "001001" "010010" "001001" "001001" ... ..$ : chr [1:3] "Spp1" "Spp2" "Spp3" > str(bestpart(m4)) int [1:50, 1:3] 0 1 0 0 1 0 0 0 0 1 ... - attr(*, "dimnames")=List of 2 ..$ : chr [1:50] "3" "2" "3" "3" ... ..$ : chr [1:3] "Spp1" "Spp2" "Spp3" > bestpart(u1c) Spp1 Spp2 Spp3 1 0.9 0.1 0.5 2 1.0 0.8 0.6 3 0.0 0.8 1.0 4 0.2 0.5 0.0 > > summary(m1) Multivariate opticut results, comb = rank, dist = poisson Call: opticut.formula(formula = Y ~ x2, strata = x0, dist = "poisson", comb = "rank") Best supported model with logLR >= 2: split assoc I mu0 mu1 logLR w Spp1 1+2 +++ 0.3772 1.969 4.353 10.23 0.9618 3 binary splits 2 species not shown > summary(m2) Multivariate opticut results, comb = all, dist = poisson Call: opticut.formula(formula = Y ~ x2, strata = x0, dist = "poisson", comb = "all") Best supported model with logLR >= 2: split assoc I mu0 mu1 logLR w Spp1 1+2 +++ 0.3772 1.969 4.353 10.23 0.9612 7 binary splits 2 species not shown > summary(m3) Multivariate opticut results, comb = NA, dist = poisson Call: opticut.formula(formula = Y ~ x2, strata = comb, dist = "poisson") Best supported model with logLR >= 2: split assoc I mu0 mu1 logLR w Spp1 e +++ 0.3772 1.969 4.353 10.23 0.9612 6 binary splits 2 species not shown > summary(m4) Multivariate opticut results, comb = rank, dist = fun Call: opticut.formula(formula = Y ~ x2, strata = x0, dist = fun) Best supported model with logLR >= 2: split assoc I mu0 mu1 logLR w Spp1 1+2 +++ 0.3772 1.969 4.353 10.23 0.9618 3 binary splits 2 species not shown > summary(u1a) Multivariate multicut uncertainty results type = asymp, B = 99, level = 0.95 split R I Lower Upper Spp1 1+2 1 0.3840 0.20839 0.5310 Spp3 1+2+3 1 0.3572 0.06558 0.6827 Spp2 2+3+4 1 0.1572 0.02199 0.3037 > summary(u2a) Multivariate multicut uncertainty results type = asymp, B = 99, level = 0.95 split R I Lower Upper Spp2 1 1 0.1449 0.02016 0.3344 Spp1 1+2 1 0.3638 0.21196 0.4776 Spp3 4 1 0.3207 0.02642 0.6977 > summary(u3a) Multivariate multicut uncertainty results type = asymp, B = 99, level = 0.95 split R I Lower Upper Spp2 a 1 0.1359 0.01219 0.2661 Spp3 d 1 0.3631 0.02222 0.7124 Spp1 e 1 0.3808 0.25061 0.5009 > > print(m1) Multivariate opticut results, comb = rank, dist = poisson Call: opticut.formula(formula = Y ~ x2, strata = x0, dist = "poisson", comb = "rank") 3 species, 3 binary splits > print(m2) Multivariate opticut results, comb = all, dist = poisson Call: opticut.formula(formula = Y ~ x2, strata = x0, dist = "poisson", comb = "all") 3 species, 7 binary splits > print(m3) Multivariate opticut results, comb = NA, dist = poisson Call: opticut.formula(formula = Y ~ x2, strata = comb, dist = "poisson") 3 species, 6 binary splits > print(m4) Multivariate opticut results, comb = rank, dist = fun Call: opticut.formula(formula = Y ~ x2, strata = x0, dist = fun) 3 species, 3 binary splits > print(u1a) Multivariate multicut uncertainty results, type = asymp, B = 99 > print(u2a) Multivariate multicut uncertainty results, type = asymp, B = 99 > print(u3a) Multivariate multicut uncertainty results, type = asymp, B = 99 > > summary(fitted(m1)) Spp1 Spp2 Spp3 Min. :1.354 Min. :3.182 Min. :0.3349 1st Qu.:1.753 1st Qu.:4.142 1st Qu.:0.6858 Median :2.242 Median :4.289 Median :0.7232 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.810 3rd Qu.:4.378 3rd Qu.:0.7543 Max. :4.923 Max. :4.518 Max. :0.8351 > summary(fitted(m2)) Spp1 Spp2 Spp3 Min. :1.354 Min. :3.182 Min. :0.3349 1st Qu.:1.753 1st Qu.:4.142 1st Qu.:0.6858 Median :2.242 Median :4.289 Median :0.7232 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.810 3rd Qu.:4.378 3rd Qu.:0.7543 Max. :4.923 Max. :4.518 Max. :0.8351 > summary(fitted(m3)) Spp1 Spp2 Spp3 Min. :1.354 Min. :3.182 Min. :0.3349 1st Qu.:1.753 1st Qu.:4.142 1st Qu.:0.6858 Median :2.242 Median :4.289 Median :0.7232 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.810 3rd Qu.:4.378 3rd Qu.:0.7543 Max. :4.923 Max. :4.518 Max. :0.8351 > f4 <- try(fitted(m4), silent=TRUE) > stopifnot(inherits(f4, "try-error")) > > summary(predict(m1)) Spp1 Spp2 Spp3 Min. :1.354 Min. :3.182 Min. :0.3349 1st Qu.:1.753 1st Qu.:4.142 1st Qu.:0.6858 Median :2.242 Median :4.289 Median :0.7232 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.810 3rd Qu.:4.378 3rd Qu.:0.7543 Max. :4.923 Max. :4.518 Max. :0.8351 > summary(predict(m2)) Spp1 Spp2 Spp3 Min. :1.354 Min. :3.182 Min. :0.3349 1st Qu.:1.753 1st Qu.:4.142 1st Qu.:0.6858 Median :2.242 Median :4.289 Median :0.7232 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.810 3rd Qu.:4.378 3rd Qu.:0.7543 Max. :4.923 Max. :4.518 Max. :0.8351 > summary(predict(m3)) Spp1 Spp2 Spp3 Min. :1.354 Min. :3.182 Min. :0.3349 1st Qu.:1.753 1st Qu.:4.142 1st Qu.:0.6858 Median :2.242 Median :4.289 Median :0.7232 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.810 3rd Qu.:4.378 3rd Qu.:0.7543 Max. :4.923 Max. :4.518 Max. :0.8351 > p4 <- try(predict(m4), silent=TRUE) > stopifnot(inherits(p4, "try-error")) > summary(predict(m1, gnew=x0, xnew=data.frame(x2=x2))) Spp1 Spp2 Spp3 Min. :1.354 Min. :3.182 Min. :0.3349 1st Qu.:1.753 1st Qu.:4.142 1st Qu.:0.6858 Median :2.242 Median :4.289 Median :0.7232 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.810 3rd Qu.:4.378 3rd Qu.:0.7543 Max. :4.923 Max. :4.518 Max. :0.8351 > pr2 <- try(predict(m2, gnew=x0, xnew=data.frame(x2=x2)), silent=TRUE) > stopifnot(inherits(pr2, "try-error")) > pr3 <- try(predict(m3, gnew=x0, xnew=data.frame(x2=x2)), silent=TRUE) > stopifnot(inherits(pr3, "try-error")) > pr4 <- try(predict(m4, gnew=x0, xnew=data.frame(x2=x2)), silent=TRUE) > stopifnot(inherits(pr4, "try-error")) > > as.data.frame(m1) split assoc I mu0 mu1 logLR w Spp1 1+2 +++ 0.3771596 1.968847 4.353309 10.2319 0.9617597 > as.data.frame(summary(m1)) split assoc I mu0 mu1 logLR w Spp1 1+2 +++ 0.3771596 1.968847 4.353309 10.2319 0.9617597 > as.data.frame(u1a) split R I Lower Upper Spp1 1+2 1 0.3839813 0.20839145 0.5309920 Spp3 1+2+3 1 0.3572272 0.06558442 0.6826718 Spp2 2+3+4 1 0.1571915 0.02199145 0.3037229 > as.data.frame(summary(u1a)) split R I Lower Upper Spp1 1+2 1 0.3839813 0.20839145 0.5309920 Spp3 1+2+3 1 0.3572272 0.06558442 0.6826718 Spp2 2+3+4 1 0.1571915 0.02199145 0.3037229 > > ## --- testing distributions --- > > ocoptions(cut=-Inf) > > ## gaussian > > summary(o <- opticut(Y ~ x2, strata=x0, dist="gaussian")) Multivariate opticut results, comb = rank, dist = gaussian Call: opticut.formula(formula = Y ~ x2, strata = x0, dist = "gaussian") Best supported models with logLR >= -Inf: split assoc I mu0 mu1 logLR w Spp2 2+3+4 + 0.4661 3.3797 4.3898 1.318 0.5037 Spp1 1+2 +++ 0.7966 2.0516 4.2299 9.876 0.9822 Spp3 2+3 + 0.1578 0.4464 0.7647 1.370 0.3359 3 binary splits > u <- uncertainty(o, type="asymp", B=9) > u <- uncertainty(o, type="boot", B=2) > u <- uncertainty(o, type="multi", B=2) > summary(fitted(o)) Spp1 Spp2 Spp3 Min. :1.009 Min. :3.150 Min. :0.4063 1st Qu.:1.729 1st Qu.:4.148 1st Qu.:0.4885 Median :2.413 Median :4.292 Median :0.7372 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.859 3rd Qu.:4.376 3rd Qu.:0.8011 Max. :4.572 Max. :4.507 Max. :0.8907 > summary(predict(o)) Spp1 Spp2 Spp3 Min. :1.009 Min. :3.150 Min. :0.4063 1st Qu.:1.729 1st Qu.:4.148 1st Qu.:0.4885 Median :2.413 Median :4.292 Median :0.7372 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.859 3rd Qu.:4.376 3rd Qu.:0.8011 Max. :4.572 Max. :4.507 Max. :0.8907 > summary(predict(o, gnew=x0, xnew=data.frame(x2=x2))) Spp1 Spp2 Spp3 Min. :1.009 Min. :3.150 Min. :0.4063 1st Qu.:1.729 1st Qu.:4.148 1st Qu.:0.4885 Median :2.413 Median :4.292 Median :0.7372 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.859 3rd Qu.:4.376 3rd Qu.:0.8011 Max. :4.572 Max. :4.507 Max. :0.8907 > getMLE(o, 1, vcov=FALSE) $coef `(Intercept)` Z1 x2 2.0516249 2.1783196 -0.5411427 attr(,"sigma") [1] 1.579978 $vcov NULL $dist [1] "gaussian" > getMLE(o, 1, vcov=TRUE) $coef `(Intercept)` Z1 x2 2.0516249 2.1783196 -0.5411427 attr(,"sigma") [1] 1.579978 $vcov `(Intercept)` Z1 x2 `(Intercept)` 0.10392305 -0.08323834 -0.03106000 Z1 -0.08323834 0.20840377 -0.02497049 x2 -0.03106000 -0.02497049 0.08413492 $dist [1] "gaussian" > > ## poisson > > summary(o <- opticut(Y ~ x2, strata=x0, dist="poisson")) Multivariate opticut results, comb = rank, dist = poisson Call: opticut.formula(formula = Y ~ x2, strata = x0, dist = "poisson") Best supported models with logLR >= -Inf: split assoc I mu0 mu1 logLR w Spp3 1+2+3 + 0.3409 0.3518 0.7158 1.056 0.3614 Spp2 2+3+4 + 0.1326 3.3632 4.3917 1.044 0.4724 Spp1 1+2 +++ 0.3772 1.9688 4.3533 10.232 0.9618 3 binary splits > u <- uncertainty(o, type="asymp", B=9) > u <- uncertainty(o, type="boot", B=2) > u <- uncertainty(o, type="multi", B=2) > summary(fitted(o)) Spp1 Spp2 Spp3 Min. :1.354 Min. :3.182 Min. :0.3349 1st Qu.:1.753 1st Qu.:4.142 1st Qu.:0.6858 Median :2.242 Median :4.289 Median :0.7232 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.810 3rd Qu.:4.378 3rd Qu.:0.7543 Max. :4.923 Max. :4.518 Max. :0.8351 > summary(predict(o)) Spp1 Spp2 Spp3 Min. :1.354 Min. :3.182 Min. :0.3349 1st Qu.:1.753 1st Qu.:4.142 1st Qu.:0.6858 Median :2.242 Median :4.289 Median :0.7232 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.810 3rd Qu.:4.378 3rd Qu.:0.7543 Max. :4.923 Max. :4.518 Max. :0.8351 > summary(predict(o, gnew=x0, xnew=data.frame(x2=x2))) Spp1 Spp2 Spp3 Min. :1.354 Min. :3.182 Min. :0.3349 1st Qu.:1.753 1st Qu.:4.142 1st Qu.:0.6858 Median :2.242 Median :4.289 Median :0.7232 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.810 3rd Qu.:4.378 3rd Qu.:0.7543 Max. :4.923 Max. :4.518 Max. :0.8351 > getMLE(o, 1, vcov=FALSE) $coef `(Intercept)` Z1 x2 0.6774482 0.7934881 -0.1945135 $vcov NULL $dist [1] "poisson" > getMLE(o, 1, vcov=TRUE) $coef `(Intercept)` Z1 x2 0.6774482 0.7934881 -0.1945135 $vcov `(Intercept)` Z1 x2 `(Intercept)` 0.020846432 -0.019077412 -0.003239102 Z1 -0.019077412 0.032241518 -0.003530473 x2 -0.003239102 -0.003530473 0.012395193 $dist [1] "poisson" > > ## binomial > > Y1 <- ifelse(Y > 0, 1, 0) > summary(o <- opticut(Y1 ~ x2, strata=x0, dist="binomial")) Multivariate opticut results, comb = rank, dist = binomial Call: opticut.formula(formula = Y1 ~ x2, strata = x0, dist = "binomial") Best supported models with logLR >= -Inf: split assoc I mu0 mu1 logLR w Spp1 1 + 1.000e+00 0.8528 1.0000 1.441 0.4355 Spp3 3 ++ 6.031e-01 0.4128 0.7396 2.286 0.5923 Spp2 4 + 2.407e-06 1.0000 1.0000 0.000 0.3333 3 binary splits > u <- uncertainty(o, type="asymp", B=9) > u <- uncertainty(o, type="boot", B=2) > u <- uncertainty(o, type="multi", B=2) > summary(fitted(o)) Spp1 Spp2 Spp3 Min. :0.8401 Min. :1 Min. :0.3807 1st Qu.:0.8484 1st Qu.:1 1st Qu.:0.4266 Median :0.8517 Median :1 Median :0.4548 Mean :0.8800 Mean :1 Mean :0.5400 3rd Qu.:0.8563 3rd Qu.:1 3rd Qu.:0.7217 Max. :1.0000 Max. :1 Max. :0.7973 > summary(predict(o)) Spp1 Spp2 Spp3 Min. :0.8401 Min. :1 Min. :0.3807 1st Qu.:0.8484 1st Qu.:1 1st Qu.:0.4266 Median :0.8517 Median :1 Median :0.4548 Mean :0.8800 Mean :1 Mean :0.5400 3rd Qu.:0.8563 3rd Qu.:1 3rd Qu.:0.7217 Max. :1.0000 Max. :1 Max. :0.7973 > summary(predict(o, gnew=x0, xnew=data.frame(x2=x2))) Spp1 Spp2 Spp3 Min. :0.8401 Min. :1 Min. :0.3807 1st Qu.:0.8484 1st Qu.:1 1st Qu.:0.4266 Median :0.8517 Median :1 Median :0.4548 Mean :0.8800 Mean :1 Mean :0.5400 3rd Qu.:0.8563 3rd Qu.:1 3rd Qu.:0.7217 Max. :1.0000 Max. :1 Max. :0.7973 > getMLE(o, 1, vcov=FALSE) $coef `(Intercept)` Z1 x2 1.7566135 16.8363082 -0.0446101 $vcov NULL $dist [1] "binomial" > getMLE(o, 1, vcov=TRUE) $coef `(Intercept)` Z1 x2 1.7566135 16.8363082 -0.0446101 $vcov `(Intercept)` Z1 x2 `(Intercept)` 0.2735044 -1.785699e-01 -0.15373362 Z1 -0.1785699 4.253702e+06 -0.03484206 x2 -0.1537336 -3.484206e-02 0.30537278 $dist [1] "binomial" > > ## negbin > > summary(o <- opticut(Y ~ x2, strata=x0, dist="negbin")) Multivariate opticut results, comb = rank, dist = negbin Call: opticut.formula(formula = Y ~ x2, strata = x0, dist = "negbin") Best supported models with logLR >= -Inf: split assoc I mu0 mu1 logLR w Spp3 1+2+3 + 0.3409 0.3518 0.7158 1.056 0.3614 Spp2 2+3+4 + 0.1326 3.3632 4.3917 1.044 0.4724 Spp1 1+2 +++ 0.3772 1.9688 4.3533 9.435 0.9617 3 binary splits There were 26 warnings (use warnings() to see them) > u <- uncertainty(o, type="asymp", B=9) Warning messages: 1: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 2: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 3: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 4: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 5: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 6: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached > u <- uncertainty(o, type="boot", B=2) There were 18 warnings (use warnings() to see them) > u <- uncertainty(o, type="multi", B=2) There were 50 or more warnings (use warnings() to see the first 50) > summary(fitted(o)) Spp1 Spp2 Spp3 Min. :1.354 Min. :3.182 Min. :0.3349 1st Qu.:1.753 1st Qu.:4.142 1st Qu.:0.6858 Median :2.242 Median :4.289 Median :0.7232 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.810 3rd Qu.:4.378 3rd Qu.:0.7543 Max. :4.923 Max. :4.518 Max. :0.8351 Warning messages: 1: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 2: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 3: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 4: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 5: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 6: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached > summary(predict(o)) Spp1 Spp2 Spp3 Min. :1.354 Min. :3.182 Min. :0.3349 1st Qu.:1.753 1st Qu.:4.142 1st Qu.:0.6858 Median :2.242 Median :4.289 Median :0.7232 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.810 3rd Qu.:4.378 3rd Qu.:0.7543 Max. :4.923 Max. :4.518 Max. :0.8351 Warning messages: 1: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 2: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 3: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 4: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 5: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 6: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached > summary(predict(o, gnew=x0, xnew=data.frame(x2=x2))) Spp1 Spp2 Spp3 Min. :1.354 Min. :3.182 Min. :0.3349 1st Qu.:1.753 1st Qu.:4.142 1st Qu.:0.6858 Median :2.242 Median :4.289 Median :0.7232 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.810 3rd Qu.:4.378 3rd Qu.:0.7543 Max. :4.923 Max. :4.518 Max. :0.8351 Warning messages: 1: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 2: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 3: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 4: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 5: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 6: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached > getMLE(o, 1, vcov=FALSE) $coef `(Intercept)` Z1 x2 0.6774487 0.7934871 -0.1945135 attr(,"theta") [1] 35394.97 $vcov NULL $dist [1] "negbin" Warning messages: 1: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 2: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached > getMLE(o, 1, vcov=TRUE) $coef `(Intercept)` Z1 x2 0.6774487 0.7934871 -0.1945135 attr(,"theta") [1] 35394.97 $vcov `(Intercept)` Z1 x2 `(Intercept)` 0.020847618 -0.019078364 -0.003239456 Z1 -0.019078364 0.032243989 -0.003530873 x2 -0.003239456 -0.003530873 0.012396289 $dist [1] "negbin" Warning messages: 1: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 2: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached > > ## beta > > Y2 <- Y / rowSums(Y) > Y2[Y2 == 0] <- 0.01 > Y2[Y2 == 1] <- 0.99 > summary(o <- opticut(Y2 ~ x2, strata=x0, dist="beta")) Multivariate opticut results, comb = rank, dist = beta Call: opticut.formula(formula = Y2 ~ x2, strata = x0, dist = "beta") Best supported models with logLR >= -Inf: split assoc I mu0 mu1 logLR w Spp1 1+2 +++ 0.5203 0.23630 0.4951 8.534 0.9823 Spp2 3+4 ++ 0.4028 0.45155 0.6592 4.933 0.8362 Spp3 3 ++ 0.2814 0.07598 0.1279 2.062 0.5841 3 binary splits > u <- uncertainty(o, type="asymp", B=9) > u <- uncertainty(o, type="boot", B=2) > u <- uncertainty(o, type="multi", B=2) > summary(fitted(o)) Spp1 Spp2 Spp3 Min. :0.1604 Min. :0.4219 Min. :0.07058 1st Qu.:0.2104 1st Qu.:0.4840 1st Qu.:0.07836 Median :0.2678 Median :0.6300 Median :0.08342 Mean :0.3286 Mean :0.5862 Mean :0.09748 3rd Qu.:0.4524 3rd Qu.:0.6842 3rd Qu.:0.12199 Max. :0.5346 Max. :0.7362 Max. :0.15100 > summary(predict(o)) Spp1 Spp2 Spp3 Min. :0.1604 Min. :0.4219 Min. :0.07058 1st Qu.:0.2104 1st Qu.:0.4840 1st Qu.:0.07836 Median :0.2678 Median :0.6300 Median :0.08342 Mean :0.3286 Mean :0.5862 Mean :0.09748 3rd Qu.:0.4524 3rd Qu.:0.6842 3rd Qu.:0.12199 Max. :0.5346 Max. :0.7362 Max. :0.15100 > summary(predict(o, gnew=x0, xnew=data.frame(x2=x2))) Spp1 Spp2 Spp3 Min. :0.1604 Min. :0.4219 Min. :0.07058 1st Qu.:0.2104 1st Qu.:0.4840 1st Qu.:0.07836 Median :0.2678 Median :0.6300 Median :0.08342 Mean :0.3286 Mean :0.5862 Mean :0.09748 3rd Qu.:0.4524 3rd Qu.:0.6842 3rd Qu.:0.12199 Max. :0.5346 Max. :0.7362 Max. :0.15100 > getMLE(o, 1, vcov=FALSE) $coef `(Intercept)` Z1 x2 (phi) -1.1730452 1.1534178 -0.2503281 4.9126628 $vcov NULL $dist [1] "beta" > getMLE(o, 1, vcov=TRUE) $coef `(Intercept)` Z1 x2 (phi) -1.1730452 1.1534178 -0.2503281 4.9126628 $vcov `(Intercept)` Z1 x2 (phi) `(Intercept)` 0.035029668 -0.029567617 -0.007726522 -0.04877094 Z1 -0.029567617 0.063846471 -0.008524017 0.04906146 x2 -0.007726522 -0.008524017 0.024865289 -0.01081438 (phi) -0.048770941 0.049061461 -0.010814377 0.86132474 $dist [1] "beta" > > ## zip > > Yzi <- Y > Yzi[1,] <- 0 > #B <- replicate(2, sample.int(n, replace=TRUE)) > B <- sapply(2:3, function(i) which((1:n) != i)) # jackknife > B[1,] <- 1 > summary(o <- opticut(Yzi ~ x2, strata=x0, dist="zip")) Multivariate opticut results, comb = rank, dist = zip Call: opticut.formula(formula = Yzi ~ x2, strata = x0, dist = "zip") Best supported models with logLR >= -Inf: split assoc I mu0 mu1 logLR w Spp3 1+2+3 + 0.3409 0.3518 0.7158 1.0563 0.3614 Spp2 2+3+4 + 0.1247 3.3647 4.3234 0.8992 0.4831 Spp1 1+2 +++ 0.3743 1.9837 4.3573 8.6714 0.8961 3 binary splits > u <- uncertainty(o, type="asymp", B=9) > u <- uncertainty(o, type="boot", B=B) > u <- uncertainty(o, type="multi", B=B) > summary(fitted(o)) Spp1 Spp2 Spp3 Min. :1.366 Min. :3.160 Min. :0.3349 1st Qu.:1.760 1st Qu.:4.048 1st Qu.:0.6858 Median :2.238 Median :4.195 Median :0.7232 Mean :2.777 Mean :4.041 Mean :0.6600 3rd Qu.:3.800 3rd Qu.:4.284 3rd Qu.:0.7543 Max. :4.882 Max. :4.424 Max. :0.8351 > summary(predict(o)) Spp1 Spp2 Spp3 Min. :1.366 Min. :3.160 Min. :0.3349 1st Qu.:1.760 1st Qu.:4.048 1st Qu.:0.6858 Median :2.238 Median :4.195 Median :0.7232 Mean :2.777 Mean :4.041 Mean :0.6600 3rd Qu.:3.800 3rd Qu.:4.284 3rd Qu.:0.7543 Max. :4.882 Max. :4.424 Max. :0.8351 > summary(predict(o, gnew=x0, xnew=data.frame(x2=x2))) Spp1 Spp2 Spp3 Min. :1.366 Min. :3.160 Min. :0.3349 1st Qu.:1.760 1st Qu.:4.048 1st Qu.:0.6858 Median :2.238 Median :4.195 Median :0.7232 Mean :2.777 Mean :4.041 Mean :0.6600 3rd Qu.:3.800 3rd Qu.:4.284 3rd Qu.:0.7543 Max. :4.882 Max. :4.424 Max. :0.8351 > getMLE(o, 1, vcov=FALSE) $coef count_`(Intercept)` count_Z1 count_x2 zero_(Intercept) 0.6849872 0.7868550 -0.1901101 -5.0354253 $vcov NULL $dist [1] "zip" > getMLE(o, 1, vcov=TRUE) $coef count_`(Intercept)` count_Z1 count_x2 zero_(Intercept) 0.6849872 0.7868550 -0.1901101 -5.0354253 $vcov count_`(Intercept)` count_Z1 count_x2 count_`(Intercept)` 0.025429054 -0.023093033 -0.001224714 count_Z1 -0.023093033 0.035749338 -0.005321943 count_x2 -0.001224714 -0.005321943 0.014137405 zero_(Intercept) 0.561944205 -0.496757613 0.310637380 zero_(Intercept) count_`(Intercept)` 0.5619442 count_Z1 -0.4967576 count_x2 0.3106374 zero_(Intercept) 76.4121113 $dist [1] "zip" > > ## zinb > > summary(o <- opticut(Yzi ~ x2, strata=x0, dist="zinb")) Multivariate opticut results, comb = rank, dist = zinb Call: opticut.formula(formula = Yzi ~ x2, strata = x0, dist = "zinb") Best supported models with logLR >= -Inf: split assoc I mu0 mu1 logLR w Spp3 1+2+3 + 0.3409 0.3518 0.7158 1.0562 0.3614 Spp2 2+3+4 + 0.1247 3.3648 4.3233 0.8992 0.4831 Spp1 1+2 +++ 0.3744 1.9832 4.3572 8.6606 0.8961 3 binary splits Warning message: In sqrt(diag(vc)[np]) : NaNs produced > u <- uncertainty(o, type="asymp", B=9) > u <- uncertainty(o, type="boot", B=B) Warning message: In sqrt(diag(vc)[np]) : NaNs produced > u <- uncertainty(o, type="multi", B=B) Warning messages: 1: In sqrt(diag(vc)[np]) : NaNs produced 2: In sqrt(diag(vc)[np]) : NaNs produced 3: In sqrt(diag(vc)[np]) : NaNs produced 4: In sqrt(diag(vc)[np]) : NaNs produced 5: In sqrt(diag(vc)[np]) : NaNs produced 6: In sqrt(diag(vc)[np]) : NaNs produced 7: In sqrt(diag(vc)[np]) : NaNs produced > summary(fitted(o)) Spp1 Spp2 Spp3 Min. :1.366 Min. :3.160 Min. :0.3349 1st Qu.:1.760 1st Qu.:4.048 1st Qu.:0.6858 Median :2.238 Median :4.195 Median :0.7232 Mean :2.777 Mean :4.041 Mean :0.6600 3rd Qu.:3.801 3rd Qu.:4.284 3rd Qu.:0.7543 Max. :4.884 Max. :4.424 Max. :0.8351 > summary(predict(o)) Spp1 Spp2 Spp3 Min. :1.366 Min. :3.160 Min. :0.3349 1st Qu.:1.760 1st Qu.:4.048 1st Qu.:0.6858 Median :2.238 Median :4.195 Median :0.7232 Mean :2.777 Mean :4.041 Mean :0.6600 3rd Qu.:3.801 3rd Qu.:4.284 3rd Qu.:0.7543 Max. :4.884 Max. :4.424 Max. :0.8351 > summary(predict(o, gnew=x0, xnew=data.frame(x2=x2))) Spp1 Spp2 Spp3 Min. :1.366 Min. :3.160 Min. :0.3349 1st Qu.:1.760 1st Qu.:4.048 1st Qu.:0.6858 Median :2.238 Median :4.195 Median :0.7232 Mean :2.777 Mean :4.041 Mean :0.6600 3rd Qu.:3.801 3rd Qu.:4.284 3rd Qu.:0.7543 Max. :4.884 Max. :4.424 Max. :0.8351 > getMLE(o, 1, vcov=FALSE) $coef count_`(Intercept)` count_Z1 count_x2 zero_(Intercept) 0.6847032 0.7871204 -0.1902707 -5.0775528 $vcov NULL $dist [1] "zinb" > getMLE(o, 1, vcov=TRUE) $coef count_`(Intercept)` count_Z1 count_x2 zero_(Intercept) 0.6847032 0.7871204 -0.1902707 -5.0775528 $vcov count_`(Intercept)` count_Z1 count_x2 count_`(Intercept)` 0.025599425 -0.023243502 -0.001099795 count_Z1 -0.023243502 0.035883636 -0.005431618 count_x2 -0.001099795 -0.005431618 0.014189671 zero_(Intercept) 0.611491196 -0.540367930 0.339598153 zero_(Intercept) count_`(Intercept)` 0.6114912 count_Z1 -0.5403679 count_x2 0.3395982 zero_(Intercept) 86.5527911 $dist [1] "zinb" > > ## zip2 > > summary(o <- opticut(Yzi ~ x2, strata=x0, dist="zip2")) Multivariate opticut results, comb = rank, dist = zip2 Call: opticut.formula(formula = Yzi ~ x2, strata = x0, dist = "zip2") Best supported models with logLR >= -Inf: split assoc I mu0 mu1 logLR w Spp3 1+2+3 + 0.9999 0.7224 1 0.3692 0.3471 Spp2 1+2+4 + 1.0000 0.9533 1 0.5788 0.4359 Spp1 1+2 + 0.9999 0.8508 1 1.7269 0.5670 3 binary splits > u <- uncertainty(o, type="asymp", B=9) > u <- uncertainty(o, type="boot", B=B) > u <- uncertainty(o, type="multi", B=B) > summary(fitted(o)) Spp1 Spp2 Spp3 Min. :2.063 Min. :3.727 Min. :0.4609 1st Qu.:2.493 1st Qu.:3.913 1st Qu.:0.6425 Median :2.770 Median :4.056 Median :0.6820 Mean :2.789 Mean :4.040 Mean :0.6601 3rd Qu.:3.048 3rd Qu.:4.135 3rd Qu.:0.7148 Max. :3.593 Max. :4.301 Max. :0.8012 > summary(predict(o)) Spp1 Spp2 Spp3 Min. :2.063 Min. :3.727 Min. :0.4609 1st Qu.:2.493 1st Qu.:3.913 1st Qu.:0.6425 Median :2.770 Median :4.056 Median :0.6820 Mean :2.789 Mean :4.040 Mean :0.6601 3rd Qu.:3.048 3rd Qu.:4.135 3rd Qu.:0.7148 Max. :3.593 Max. :4.301 Max. :0.8012 > summary(predict(o, gnew=x0, xnew=data.frame(x2=x2))) Spp1 Spp2 Spp3 Min. :2.063 Min. :3.727 Min. :0.4609 1st Qu.:2.493 1st Qu.:3.913 1st Qu.:0.6425 Median :2.770 Median :4.056 Median :0.6820 Mean :2.789 Mean :4.040 Mean :0.6601 3rd Qu.:3.048 3rd Qu.:4.135 3rd Qu.:0.7148 Max. :3.593 Max. :4.301 Max. :0.8012 > getMLE(o, 1, vcov=FALSE) $coef zero_ZZ(Intercept) zero_ZZZ1 count_X(Intercept) count_Xx2 1.7411170 10.2218839 1.1817729 -0.1535508 $vcov NULL $dist [1] "zip2" > getMLE(o, 1, vcov=TRUE) $coef zero_ZZ(Intercept) zero_ZZZ1 count_X(Intercept) count_Xx2 1.7411170 10.2218839 1.1817729 -0.1535508 $vcov zero_ZZ(Intercept) zero_ZZZ1 count_X(Intercept) zero_ZZ(Intercept) 0.383870384 -3.792870e-01 0.008310887 zero_ZZZ1 -0.379287023 1.230008e+04 -0.012476152 count_X(Intercept) 0.008310887 -1.247615e-02 0.010412712 count_Xx2 -0.002669227 4.290973e-02 -0.005757329 count_Xx2 zero_ZZ(Intercept) -0.002669227 zero_ZZZ1 0.042909732 count_X(Intercept) -0.005757329 count_Xx2 0.012423049 $dist [1] "zip2" > > ## zinb2 > > summary(o <- opticut(Yzi ~ x2, strata=x0, dist="zinb2")) Multivariate opticut results, comb = rank, dist = zinb2 Call: opticut.formula(formula = Yzi ~ x2, strata = x0, dist = "zinb2") Best supported models with logLR >= -Inf: split assoc I mu0 mu1 logLR w Spp3 1+2+3 + 1 0.7225 1 0.3693 0.3471 Spp2 1+2+4 + 1 0.9533 1 0.5788 0.4359 Spp1 1+2 + 1 0.8508 1 1.7165 0.5659 3 binary splits Warning messages: 1: In sqrt(diag(vc)[np]) : NaNs produced 2: In sqrt(diag(vc)[np]) : NaNs produced 3: In sqrt(diag(vc)[np]) : NaNs produced > u <- uncertainty(o, type="asymp", B=9) > u <- uncertainty(o, type="boot", B=B) Warning messages: 1: In sqrt(diag(vc)[np]) : NaNs produced 2: In sqrt(diag(vc)[np]) : NaNs produced 3: In sqrt(diag(vc)[np]) : NaNs produced 4: In sqrt(diag(vc)[np]) : NaNs produced > u <- uncertainty(o, type="multi", B=B) Warning messages: 1: In sqrt(diag(vc)[np]) : NaNs produced 2: In sqrt(diag(vc)[np]) : NaNs produced 3: In sqrt(diag(vc)[np]) : NaNs produced 4: In sqrt(diag(vc)[np]) : NaNs produced 5: In sqrt(diag(vc)[np]) : NaNs produced > summary(fitted(o)) Spp1 Spp2 Spp3 Min. :2.063 Min. :3.727 Min. :0.4610 1st Qu.:2.493 1st Qu.:3.913 1st Qu.:0.6425 Median :2.770 Median :4.056 Median :0.6820 Mean :2.789 Mean :4.040 Mean :0.6601 3rd Qu.:3.048 3rd Qu.:4.135 3rd Qu.:0.7148 Max. :3.593 Max. :4.301 Max. :0.8012 > summary(predict(o)) Spp1 Spp2 Spp3 Min. :2.063 Min. :3.727 Min. :0.4610 1st Qu.:2.493 1st Qu.:3.913 1st Qu.:0.6425 Median :2.770 Median :4.056 Median :0.6820 Mean :2.789 Mean :4.040 Mean :0.6601 3rd Qu.:3.048 3rd Qu.:4.135 3rd Qu.:0.7148 Max. :3.593 Max. :4.301 Max. :0.8012 > summary(predict(o, gnew=x0, xnew=data.frame(x2=x2))) Spp1 Spp2 Spp3 Min. :2.063 Min. :3.727 Min. :0.4610 1st Qu.:2.493 1st Qu.:3.913 1st Qu.:0.6425 Median :2.770 Median :4.056 Median :0.6820 Mean :2.789 Mean :4.040 Mean :0.6601 3rd Qu.:3.048 3rd Qu.:4.135 3rd Qu.:0.7148 Max. :3.593 Max. :4.301 Max. :0.8012 > getMLE(o, 1, vcov=FALSE) $coef zero_ZZ(Intercept) zero_ZZZ1 count_X(Intercept) count_Xx2 1.7411203 11.1979779 1.1817712 -0.1535531 $vcov NULL $dist [1] "zinb2" > getMLE(o, 1, vcov=TRUE) $coef zero_ZZ(Intercept) zero_ZZZ1 count_X(Intercept) count_Xx2 1.7411203 11.1979779 1.1817712 -0.1535531 $vcov zero_ZZ(Intercept) zero_ZZZ1 count_X(Intercept) zero_ZZ(Intercept) 0.383925932 -3.792715e-01 0.008315533 zero_ZZZ1 -0.379271486 3.264459e+04 -0.012475106 count_X(Intercept) 0.008315533 -1.247511e-02 0.010413640 count_Xx2 -0.002670217 4.291023e-02 -0.005757730 count_Xx2 zero_ZZ(Intercept) -0.002670217 zero_ZZZ1 0.042910228 count_X(Intercept) -0.005757730 count_Xx2 0.012423524 $dist [1] "zinb2" > > ## rsf > > library(ResourceSelection) ResourceSelection 0.3-6 2023-06-27 > n.used <- 1000 > m <- 1 > n <- n.used * m > set.seed(1234) > x <- data.frame(x0=as.factor(sample(1:3, n, replace=TRUE)), + x1=rnorm(n), x2=runif(n)) > cfs <- c(1, -0.5, 0.1, -1, 0.5) > dd <- simulateUsedAvail(x, cfs, n.used, m, link="logit") > > Y <- dd$status > X <- model.matrix(~ x1 + x2, dd) > > ## intercept + partition + covariates > o <- opticut(Y ~ x1 + x2, dd, strata=dd$x0, dist="rsf") > o$species $`Sp 1` Univariate opticut results, comb = rank, dist = rsf I = 0.1533; w = 0.998; H = 0.996; logL_null = -6896 Best supported models with logLR >= -Inf: assoc I mu0 mu1 logLR w 3 +++ 0.1533 1 1.362 11.192 0.997994 1+3 ++ 0.1071 1 1.240 4.982 0.002006 2 binary splits > u <- uncertainty(o, type="asymp", B=9) > u <- uncertainty(o, type="boot", B=2) > u <- uncertainty(o, type="multi", B=2) > summary(fitted(o)) Sp 1 Min. :0.754 1st Qu.:1.081 Median :1.226 Mean :1.298 3rd Qu.:1.492 Max. :2.396 > summary(predict(o)) Sp 1 Min. :0.754 1st Qu.:1.081 Median :1.226 Mean :1.298 3rd Qu.:1.492 Max. :2.396 > summary(predict(o, gnew=o$strata, xnew=o$X)) Sp 1 Min. :0.754 1st Qu.:1.081 Median :1.226 Mean :1.298 3rd Qu.:1.492 Max. :2.396 > getMLE(o, 1, vcov=FALSE) $coef (Intercept) Z1 x1 x2 0.0000000 0.3090649 -0.1496463 0.2161110 $vcov NULL $dist [1] "rsf" > getMLE(o, 1, vcov=TRUE) $coef (Intercept) Z1 x1 x2 0.0000000 0.3090649 -0.1496463 0.2161110 $vcov (Intercept) Z1 x1 x2 (Intercept) NA NA NA NA Z1 NA 4.168452e-03 5.152092e-05 -0.0003378792 x1 NA 5.152092e-05 1.050105e-03 -0.0002191716 x2 NA -3.378792e-04 -2.191716e-04 0.0113575726 $dist [1] "rsf" > ## intercept + partition > o <- opticut(Y ~ 1, dd, strata=x0, dist="rsf") > o$species $`Sp 1` Univariate opticut results, comb = rank, dist = rsf I = 0.1585; w = 0.9994; H = 0.9988; logL_null = -6909 Best supported models with logLR >= -Inf: assoc I mu0 mu1 logLR w 3 +++ 0.1585 1 1.377 11.975 0.9993852 1+3 ++ 0.1029 1 1.229 4.581 0.0006148 2 binary splits > u <- uncertainty(o, type="asymp", B=9) > u <- uncertainty(o, type="boot", B=2) > u <- uncertainty(o, type="multi", B=2) > summary(fitted(o)) Sp 1 Min. :1.000 1st Qu.:1.000 Median :1.000 Mean :1.138 3rd Qu.:1.377 Max. :1.377 > summary(predict(o)) Sp 1 Min. :1.000 1st Qu.:1.000 Median :1.000 Mean :1.138 3rd Qu.:1.377 Max. :1.377 > summary(predict(o, gnew=o$strata, xnew=NULL)) Sp 1 Min. :1.000 1st Qu.:1.000 Median :1.000 Mean :1.138 3rd Qu.:1.377 Max. :1.377 > getMLE(o, 1, vcov=FALSE) $coef (Intercept) Z1 0.0000000 0.3197308 $vcov NULL $dist [1] "rsf" > getMLE(o, 1, vcov=TRUE) $coef (Intercept) Z1 0.0000000 0.3197308 $vcov (Intercept) Z1 (Intercept) NA NA Z1 NA 0.004156167 $dist [1] "rsf" > > ## rspf > > o <- opticut(Y ~ x1 + x2, dd, strata=x0, dist="rspf") > o$species $`Sp 1` Univariate opticut results, comb = rank, dist = rspf I = 0.2751; w = 0.9762; H = 0.9535; logL_null = -6889 Best supported models with logLR >= -Inf: assoc I mu0 mu1 logLR w 3 ++ 0.2751 0.3493 0.4857 4.6531 0.97617 1+3 + 0.2362 0.9126 0.9441 0.9403 0.02383 2 binary splits > u <- uncertainty(o, type="asymp", B=9) > u <- uncertainty(o, type="boot", B=2) > u <- uncertainty(o, type="multi", B=2) > summary(fitted(o)) Sp 1 Min. :0.2303 1st Qu.:0.3810 Median :0.4424 Mean :0.4562 3rd Qu.:0.5331 Max. :0.7434 > summary(predict(o)) Sp 1 Min. :0.2303 1st Qu.:0.3810 Median :0.4424 Mean :0.4562 3rd Qu.:0.5331 Max. :0.7434 > summary(predict(o, gnew=o$strata, xnew=o$X)) Sp 1 Min. :0.2303 1st Qu.:0.3810 Median :0.4424 Mean :0.4562 3rd Qu.:0.5331 Max. :0.7434 > getMLE(o, 1, vcov=FALSE) $coef (Intercept) Z1 x1 x2 -0.6222194 0.5648166 -0.3042349 0.4060796 $vcov NULL $dist [1] "rspf" > getMLE(o, 1, vcov=TRUE) $coef (Intercept) Z1 x1 x2 -0.6222194 0.5648166 -0.3042349 0.4060796 $vcov (Intercept) Z1 x1 x2 (Intercept) 1.2575477 0.35181230 -0.23394056 0.2524899 Z1 0.3518123 0.11885164 -0.06704858 0.0777048 x1 -0.2339406 -0.06704858 0.04856497 -0.0522358 x2 0.2524899 0.07770480 -0.05223580 0.1007597 $dist [1] "rspf" > > ## --- ... in uncertainty should produce an error --- > > set.seed(1234) > n <- 50 > x0 <- sample(1:4, n, TRUE) > x1 <- ifelse(x0 %in% 1:2, 1, 0) > x2 <- rnorm(n, 0.5, 1) > lam <- exp(0.5 + 1*x1 + -0.2*x2) > A <- ifelse(x0 %in% c(1,3), 1, 10) > Y <- rpois(n, lam*A) > > ## no offset: incorrect > no <- opticut(Y ~ x2, strata=x0, dist="poisson", comb="rank") > ## with offsets: log Area > wo <- opticut(Y ~ x2, strata=x0, dist="poisson", + offset=log(A), weights=rep(1,n), comb="rank") > no$species[[1]] Univariate opticut results, comb = rank, dist = poisson I = 0.6803; w = 1; H = 1; logL_null = -550.9 Best supported models with logLR >= -Inf: assoc I mu0 mu1 logLR w 2 +++ 0.6803 9.255 48.65 334.7 1.000e+00 2+4 +++ 0.8345 2.829 31.36 278.9 6.070e-25 1+2+4 +++ 0.8615 1.971 26.49 168.9 9.954e-73 3 binary splits > wo$species[[1]] Univariate opticut results, comb = rank, dist = poisson I = 0.5023; w = 0.9997; H = 0.9995; logL_null = -272.3 Best supported models with logLR >= -Inf: assoc I mu0 mu1 logLR w 1+2 +++ 0.5023 1.604 4.841 141.7 9.997e-01 2 +++ 0.4834 1.698 4.875 133.3 2.240e-04 1+2+3 +++ 0.4938 1.580 4.663 131.5 3.882e-05 3 binary splits > > nu <- uncertainty(no, type="multi", B=2) > ## passing ... is not enough for resampling, treated as user error > wu <- try(uncertainty(wo, type="multi", B=2), silent=TRUE) > stopifnot(inherits(wu, "try-error")) > > ## --- zip2 & zinb2 coef inversion > > ## implementation: > ## - MLE returns modified coefs (P of 1 in ZI) > ## - .opticut1 returns: > ## -1*coef[1:2] > ## linkinv: binomial(link)$linkinv(eta) > ## - asymp uncertainty uses MLE > ## bestmodel returns unmodified coefs (P of 0 in ZI) > > ## less 0 in g=1 stratum: assoc is 1+ or 2- > yzi <- c(rep(0, 10), rpois(40, 6), rep(0, 30), rpois(20, 4)) > g <- rep(1:2, each=50) > table(yzi, g) g yzi 1 2 0 11 30 1 0 1 2 1 4 3 4 6 4 6 1 5 7 3 6 4 2 7 8 2 8 4 1 9 1 0 10 2 0 11 1 0 12 1 0 > o1 <- opticut(yzi, strata=g, dist="zip2") > o2 <- opticut(yzi, strata=g, dist="zinb2") > ## assoc must be positive for comb=rank > stopifnot(o1$species[[1]]$assoc == 1) > stopifnot(o2$species[[1]]$assoc == 1) > ## MLE is negative (prob of 0) > stopifnot(getMLE(o1, 1)$coef[2] >= 0) > stopifnot(getMLE(o2, 1)$coef[2] >= 0) > stopifnot(coef(bestmodel(o1, 1)[[1]], "zero")[2] < 0) > stopifnot(coef(bestmodel(o2, 1)[[1]], "zero")[2] < 0) > ## uncertainty should show 1+ (mu0 < mu1) > o1$species[[1]] Univariate opticut results, comb = rank, dist = zip2 I = 0.6871; w = 1; H = 1; logL_null = -201.2 Best supported model with logLR >= -Inf: assoc I mu0 mu1 logLR w 1 ++ 0.6871 0.4019 0.7837 7.69 1 1 binary split > o2$species[[1]] Univariate opticut results, comb = rank, dist = zinb2 I = 0.6884; w = 1; H = 1; logL_null = -201 Best supported model with logLR >= -Inf: assoc I mu0 mu1 logLR w 1 ++ 0.6884 0.4026 0.785 7.69 1 1 binary split > u <- uncertainty(o1, type="asymp", B=9)$uncertainty[[1]] > stopifnot(all(u$mu0 < u$mu1)) > u <- uncertainty(o2, type="asymp", B=9)$uncertainty[[1]] > stopifnot(all(u$mu0 < u$mu1)) > B <- sapply(2:10, function(i) which((1:length(g)) != i)) # jackknife > u <- uncertainty(o1, type="boot", B=B)$uncertainty[[1]] > stopifnot(all(u$mu0 < u$mu1)) > u <- uncertainty(o2, type="boot", B=B)$uncertainty[[1]] > stopifnot(all(u$mu0 < u$mu1)) > u <- uncertainty(o1, type="multi", B=B)$uncertainty[[1]] > stopifnot(all(u$mu0 < u$mu1)) > u <- uncertainty(o2, type="multi", B=B)$uncertainty[[1]] > stopifnot(all(u$mu0 < u$mu1)) > > proc.time() user system elapsed 22.26 0.92 23.18