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) > m1 <- multicut(Y ~ x2, strata=x0, 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 <- multicut(Y ~ x2, strata=x0, dist=fun) > ocoptions(try_error=FALSE) > > subset(m1, c(3,1)) Multivariate multicut results, dist = poisson Call: multicut.formula(formula = Y ~ x2, strata = x0, dist = "poisson") 2 species > subset(m1, c(TRUE, FALSE, TRUE)) Multivariate multicut results, dist = poisson Call: multicut.formula(formula = Y ~ x2, strata = x0, dist = "poisson") 2 species > subset(m1, c("Spp1", "Spp3")) Multivariate multicut results, dist = poisson Call: multicut.formula(formula = Y ~ x2, strata = x0, dist = "poisson") 2 species > > str(m1$strata) Factor w/ 4 levels "1","2","3","4": 3 2 3 3 2 3 4 3 3 2 ... > > opticut:::.extractOpticut(m1, Z=model.matrix(~m1$strata)[,-1]) $Spp1 $Spp1$coef `(Intercept)` `m1$strata2` `m1$strata3` `m1$strata4` x2 1.44494130 0.05330905 -0.90728235 -0.58512568 -0.20224463 $Spp1$logLik [1] -90.00833 $Spp1$linkinv NULL $Spp2 $Spp2$coef `(Intercept)` `m1$strata2` `m1$strata3` `m1$strata4` x2 1.21028464 0.25189802 0.30228354 0.23281645 -0.02784977 $Spp2$logLik [1] -98.85274 $Spp2$linkinv NULL $Spp3 $Spp3$coef `(Intercept)` `m1$strata2` `m1$strata3` `m1$strata4` x2 -0.5782195 0.1282529 0.4076079 -0.4856141 0.1082685 $Spp3$logLik [1] -49.334 $Spp3$linkinv NULL > opticut:::.extractOpticut(m4, Z=model.matrix(~m1$strata)[,-1]) $Spp1 $Spp1$coef `(Intercept)` `m1$strata2` `m1$strata3` `m1$strata4` x2 1.44494130 0.05330905 -0.90728235 -0.58512568 -0.20224463 $Spp1$logLik 'log Lik.' -90.00833 (df=5) $Spp1$linkinv function (eta) pmax(exp(eta), .Machine$double.eps) $Spp2 $Spp2$coef `(Intercept)` `m1$strata2` `m1$strata3` `m1$strata4` x2 1.21028464 0.25189802 0.30228354 0.23281645 -0.02784977 $Spp2$logLik 'log Lik.' -98.85274 (df=5) $Spp2$linkinv function (eta) pmax(exp(eta), .Machine$double.eps) $Spp3 $Spp3$coef `(Intercept)` `m1$strata2` `m1$strata3` `m1$strata4` x2 -0.5782195 0.1282529 0.4076079 -0.4856141 0.1082685 $Spp3$logLik 'log Lik.' -49.334 (df=5) $Spp3$linkinv function (eta) pmax(exp(eta), .Machine$double.eps) > > u1a <- uncertainty(m1, type="asymp", B=99) > u4a <- try(uncertainty(m4, type="asymp", B=999), silent=TRUE) > stopifnot(inherits(u4a, "try-error")) > > u1b <- uncertainty(m1, type="boot", B=9) > u4b <- uncertainty(m4, type="boot", B=9) > > summary(subset(u1a, c(3,1))) Multivariate multicut uncertainty results type = asymp, B = 99, level = 0.95 split R I Lower Upper Spp1 1+2 0.89 0.4556 0.3049 0.5969 Spp3 2+3 0.34 0.5096 0.2887 0.7382 > summary(subset(u1b, 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 0.4944 0.3594 0.6965 Spp3 3 0.3 0.3928 0.3049 0.4566 > summary(subset(u1b, c("Spp1", "Spp3"))) Multivariate multicut uncertainty results type = boot, B = 9, level = 0.95 split R I Lower Upper Spp1 1+2 1.0 0.4944 0.3594 0.6965 Spp3 3 0.3 0.3928 0.3049 0.4566 > > 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(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) [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(u4b) [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 > > > bestmodel(m1) $Spp1 Call: stats::glm(formula = Y ~ . - 1, family = Family, data = XX) Coefficients: `(Intercept)` Z2 Z3 Z4 x2 1.44494 0.05331 -0.90728 -0.58513 -0.20224 Degrees of Freedom: 50 Total (i.e. Null); 45 Residual Null Deviance: 181.3 Residual Deviance: 52.42 AIC: 190 $Spp2 Call: stats::glm(formula = Y ~ . - 1, family = Family, data = XX) Coefficients: `(Intercept)` Z2 Z3 Z4 x2 1.21028 0.25190 0.30228 0.23282 -0.02785 Degrees of Freedom: 50 Total (i.e. Null); 45 Residual Null Deviance: 311.9 Residual Deviance: 38.11 AIC: 207.7 $Spp3 Call: stats::glm(formula = Y ~ . - 1, family = Family, data = XX) Coefficients: `(Intercept)` Z2 Z3 Z4 x2 -0.5782 0.1283 0.4076 -0.4856 0.1083 Degrees of Freedom: 50 Total (i.e. Null); 45 Residual Null Deviance: 50.64 Residual Deviance: 40.99 AIC: 108.7 > ## 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)` Z2 Z3 Z4 x2 1.44494130 0.05330905 -0.90728235 -0.58512568 -0.20224463 $vcov NULL $dist [1] "poisson" > getMLE(m4, 1, vcov=FALSE) $coef `(Intercept)` Z2 Z3 Z4 x2 1.44494130 0.05330905 -0.90728235 -0.58512568 -0.20224463 $vcov NULL $dist [1] "fun" > getMLE(m1, 1, vcov=TRUE) $coef `(Intercept)` Z2 Z3 Z4 x2 1.44494130 0.05330905 -0.90728235 -0.58512568 -0.20224463 $vcov `(Intercept)` Z2 Z3 Z4 x2 `(Intercept)` 0.029337380 -0.025789183 -0.028086747 -0.027278218 -0.006124488 Z2 -0.025789183 0.046015410 0.026007145 0.026148056 -0.001067382 Z3 -0.028086747 0.026007145 0.065815280 0.026879870 0.003589569 Z4 -0.027278218 0.026148056 0.026879870 0.068288987 0.001950756 x2 -0.006124488 -0.001067382 0.003589569 0.001950756 0.012413778 $dist [1] "poisson" > mle4 <- try(getMLE(m4, 1, vcov=TRUE), silent=TRUE) > stopifnot(inherits(mle4, "try-error")) > > summary(m1) Multivariate multticut results, dist = poisson Call: multicut.formula(formula = Y ~ x2, strata = x0, dist = "poisson") Species model with logLR >= 2: split assoc I logLR Spp1 1+2 +++ 0.4465 10.9 2 species not shown > summary(m4) Multivariate multticut results, dist = fun Call: multicut.formula(formula = Y ~ x2, strata = x0, dist = fun) Species model with logLR >= 2: split assoc I logLR Spp1 1+2 +++ 0.4465 10.9 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 0.89 0.4556 0.30490 0.5969 Spp3 2+3 0.34 0.5096 0.28869 0.7382 Spp2 2+3 0.27 0.2204 0.06132 0.3607 > summary(u1b) Multivariate multicut uncertainty results type = boot, B = 9, level = 0.95 split R I Lower Upper Spp1 1+2 1.0 0.4944 0.35942 0.6965 Spp2 2+3 0.4 0.1687 0.07778 0.2798 Spp3 3 0.3 0.3928 0.30492 0.4566 > summary(u4b) Multivariate multicut uncertainty results type = boot, B = 9, level = 0.95 split R I Lower Upper Spp1 1+2 0.8 0.4746 0.3847 0.5732 Spp3 2+3 0.4 0.5091 0.3992 0.7056 Spp2 2+3+4 0.4 0.2189 0.1545 0.2728 > > print(m1) Multivariate multicut results, dist = poisson Call: multicut.formula(formula = Y ~ x2, strata = x0, dist = "poisson") 3 species > print(m4) Multivariate multicut results, dist = fun Call: multicut.formula(formula = Y ~ x2, strata = x0, dist = fun) 3 species > print(u1a) Multivariate multicut uncertainty results, type = asymp, B = 99 > print(u1b) Multivariate multicut uncertainty results, type = boot, B = 9 > print(u4a) [1] "Error in .opticut1(Y = object$Y[j, i], X = object$X[j, , drop = FALSE], : \n Unable to return full model for custom distribution function.\n" attr(,"class") [1] "try-error" attr(,"condition") > print(u4b) Multivariate multicut uncertainty results, type = boot, B = 9 > > summary(fitted(m1)) Spp1 Spp2 Spp3 Min. :1.215 Min. :3.198 Min. :0.3200 1st Qu.:1.835 1st Qu.:4.104 1st Qu.:0.5528 Median :2.643 Median :4.244 Median :0.6718 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.824 3rd Qu.:4.433 3rd Qu.:0.8069 Max. :5.084 Max. :4.650 Max. :1.0131 > f4 <- try(fitted(m4), silent=TRUE) > stopifnot(inherits(f4, "try-error")) > > summary(predict(m1)) Spp1 Spp2 Spp3 Min. :1.215 Min. :3.198 Min. :0.3200 1st Qu.:1.835 1st Qu.:4.104 1st Qu.:0.5528 Median :2.643 Median :4.244 Median :0.6718 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.824 3rd Qu.:4.433 3rd Qu.:0.8069 Max. :5.084 Max. :4.650 Max. :1.0131 > 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.215 Min. :3.198 Min. :0.3200 1st Qu.:1.835 1st Qu.:4.104 1st Qu.:0.5528 Median :2.643 Median :4.244 Median :0.6718 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.824 3rd Qu.:4.433 3rd Qu.:0.8069 Max. :5.084 Max. :4.650 Max. :1.0131 > pr4 <- try(predict(m4, gnew=x0, xnew=data.frame(x2=x2)), silent=TRUE) > stopifnot(inherits(pr4, "try-error")) > > as.data.frame(m1) data frame with 0 columns and 0 rows > as.data.frame(summary(m1)) data frame with 0 columns and 0 rows > > as.data.frame(u1a) split R I Lower Upper Spp1 1+2 0.89 0.4556349 0.30489703 0.5968533 Spp3 2+3 0.34 0.5096154 0.28868843 0.7381955 Spp2 2+3 0.27 0.2204106 0.06132147 0.3606612 > as.data.frame(summary(u1a)) split R I Lower Upper Spp1 1+2 0.89 0.4556349 0.30489703 0.5968533 Spp3 2+3 0.34 0.5096154 0.28868843 0.7381955 Spp2 2+3 0.27 0.2204106 0.06132147 0.3606612 > > ## --- testing distributions --- > > ocoptions(cut=-Inf) > > ## gaussian > > summary(o <- multicut(Y ~ x2, strata=x0, dist="gaussian")) Multivariate multticut results, dist = gaussian Call: multicut.formula(formula = Y ~ x2, strata = x0, dist = "gaussian") Species models with logLR >= -Inf: split assoc I logLR Spp2 2+3+4 + 0.5254 1.430 Spp1 1+2 +++ 0.8523 10.486 Spp3 2+3 ++ 0.2549 2.046 > summary(fitted(o)) Spp1 Spp2 Spp3 Min. :0.843 Min. :3.171 Min. :0.2807 1st Qu.:1.900 1st Qu.:4.106 1st Qu.:0.5474 Median :2.750 Median :4.245 Median :0.6761 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.892 3rd Qu.:4.439 3rd Qu.:0.8184 Max. :4.688 Max. :4.637 Max. :0.9747 > summary(predict(o)) Spp1 Spp2 Spp3 Min. :0.843 Min. :3.171 Min. :0.2807 1st Qu.:1.900 1st Qu.:4.106 1st Qu.:0.5474 Median :2.750 Median :4.245 Median :0.6761 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.892 3rd Qu.:4.439 3rd Qu.:0.8184 Max. :4.688 Max. :4.637 Max. :0.9747 > summary(predict(o, gnew=x0, xnew=data.frame(x2=x2))) Spp1 Spp2 Spp3 Min. :0.843 Min. :3.171 Min. :0.2807 1st Qu.:1.900 1st Qu.:4.106 1st Qu.:0.5474 Median :2.750 Median :4.245 Median :0.6761 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.892 3rd Qu.:4.439 3rd Qu.:0.8184 Max. :4.688 Max. :4.637 Max. :0.9747 > getMLE(o, 1, vcov=FALSE) $coef `(Intercept)` Z2 Z3 Z4 x2 4.1361005 0.1944123 -2.3344584 -1.6990377 -0.5652299 $vcov NULL $dist [1] "gaussian" > getMLE(o, 1, vcov=TRUE) $coef `(Intercept)` Z2 Z3 Z4 x2 4.1361005 0.1944123 -2.3344584 -1.6990377 -0.5652299 $vcov `(Intercept)` Z2 Z3 Z4 x2 `(Intercept)` 0.28498668 -0.24795983 -0.2689338 -0.26179063 -0.05136681 Z2 -0.24795983 0.45154367 0.2513670 0.25288306 -0.01090231 Z3 -0.26893385 0.25136695 0.4203444 0.25792879 0.02437030 Z4 -0.26179063 0.25288306 0.2579288 0.48752183 0.01235735 x2 -0.05136681 -0.01090231 0.0243703 0.01235735 0.08638505 $dist [1] "gaussian" > ocoptions(fix_fitted=TRUE) > u <- uncertainty(o, type="asymp", B=9) > u <- uncertainty(o, type="boot", B=2) > ocoptions(fix_fitted=FALSE) > > ## poisson > > summary(o <- multicut(Y ~ x2, strata=x0, dist="poisson")) Multivariate multticut results, dist = poisson Call: multicut.formula(formula = Y ~ x2, strata = x0, dist = "poisson") Species models with logLR >= -Inf: split assoc I logLR Spp2 2+3+4 + 0.1500 1.120 Spp1 1+2 +++ 0.4465 10.903 Spp3 2+3 + 0.4191 1.468 > summary(fitted(o)) Spp1 Spp2 Spp3 Min. :1.215 Min. :3.198 Min. :0.3200 1st Qu.:1.835 1st Qu.:4.104 1st Qu.:0.5528 Median :2.643 Median :4.244 Median :0.6718 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.824 3rd Qu.:4.433 3rd Qu.:0.8069 Max. :5.084 Max. :4.650 Max. :1.0131 > summary(predict(o)) Spp1 Spp2 Spp3 Min. :1.215 Min. :3.198 Min. :0.3200 1st Qu.:1.835 1st Qu.:4.104 1st Qu.:0.5528 Median :2.643 Median :4.244 Median :0.6718 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.824 3rd Qu.:4.433 3rd Qu.:0.8069 Max. :5.084 Max. :4.650 Max. :1.0131 > summary(predict(o, gnew=x0, xnew=data.frame(x2=x2))) Spp1 Spp2 Spp3 Min. :1.215 Min. :3.198 Min. :0.3200 1st Qu.:1.835 1st Qu.:4.104 1st Qu.:0.5528 Median :2.643 Median :4.244 Median :0.6718 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.824 3rd Qu.:4.433 3rd Qu.:0.8069 Max. :5.084 Max. :4.650 Max. :1.0131 > getMLE(o, 1, vcov=FALSE) $coef `(Intercept)` Z2 Z3 Z4 x2 1.44494130 0.05330905 -0.90728235 -0.58512568 -0.20224463 $vcov NULL $dist [1] "poisson" > getMLE(o, 1, vcov=TRUE) $coef `(Intercept)` Z2 Z3 Z4 x2 1.44494130 0.05330905 -0.90728235 -0.58512568 -0.20224463 $vcov `(Intercept)` Z2 Z3 Z4 x2 `(Intercept)` 0.029337380 -0.025789183 -0.028086747 -0.027278218 -0.006124488 Z2 -0.025789183 0.046015410 0.026007145 0.026148056 -0.001067382 Z3 -0.028086747 0.026007145 0.065815280 0.026879870 0.003589569 Z4 -0.027278218 0.026148056 0.026879870 0.068288987 0.001950756 x2 -0.006124488 -0.001067382 0.003589569 0.001950756 0.012413778 $dist [1] "poisson" > u <- uncertainty(o, type="asymp", B=9) > u <- uncertainty(o, type="boot", B=2) > > ## binomial > > Y1 <- ifelse(Y > 0, 1, 0) > summary(o <- multicut(Y1 ~ x2, strata=x0, dist="binomial")) Multivariate multticut results, dist = binomial Call: multicut.formula(formula = Y1 ~ x2, strata = x0, dist = "binomial") Species models with logLR >= -Inf: split assoc I logLR Spp2 1+2+3+4 + 2.206e-06 0.000 Spp1 1+2 + 1.000e+00 1.924 Spp3 3 ++ 6.881e-01 2.477 > summary(fitted(o)) Spp1 Spp2 Spp3 Min. :0.7761 Min. :1 Min. :0.3183 1st Qu.:0.8116 1st Qu.:1 1st Qu.:0.4262 Median :0.8426 Median :1 Median :0.4913 Mean :0.8800 Mean :1 Mean :0.5400 3rd Qu.:0.9308 3rd Qu.:1 3rd Qu.:0.7244 Max. :1.0000 Max. :1 Max. :0.7929 > summary(predict(o)) Spp1 Spp2 Spp3 Min. :0.7761 Min. :1 Min. :0.3183 1st Qu.:0.8116 1st Qu.:1 1st Qu.:0.4262 Median :0.8426 Median :1 Median :0.4913 Mean :0.8800 Mean :1 Mean :0.5400 3rd Qu.:0.9308 3rd Qu.:1 3rd Qu.:0.7244 Max. :1.0000 Max. :1 Max. :0.7929 > summary(predict(o, gnew=x0, xnew=data.frame(x2=x2))) Spp1 Spp2 Spp3 Min. :0.7761 Min. :1 Min. :0.3183 1st Qu.:0.8116 1st Qu.:1 1st Qu.:0.4262 Median :0.8426 Median :1 Median :0.4913 Mean :0.8800 Mean :1 Mean :0.5400 3rd Qu.:0.9308 3rd Qu.:1 3rd Qu.:0.7244 Max. :1.0000 Max. :1 Max. :0.7929 > getMLE(o, 1, vcov=FALSE) $coef `(Intercept)` Z2 Z3 Z4 x2 18.6685878 -16.0567498 -17.1461224 -17.0852228 -0.1646947 $vcov NULL $dist [1] "binomial" > getMLE(o, 1, vcov=TRUE) $coef `(Intercept)` Z2 Z3 Z4 x2 18.6685878 -16.0567498 -17.1461224 -17.0852228 -0.1646947 $vcov `(Intercept)` Z2 Z3 Z4 `(Intercept)` 4.244261e+06 -4.244261e+06 -4.244261e+06 -4.244261e+06 Z2 -4.244261e+06 4.244262e+06 4.244261e+06 4.244261e+06 Z3 -4.244261e+06 4.244261e+06 4.244261e+06 4.244261e+06 Z4 -4.244261e+06 4.244261e+06 4.244261e+06 4.244261e+06 x2 -2.304736e-01 -4.775448e-02 1.053559e-01 5.718296e-02 x2 `(Intercept)` -0.23047361 Z2 -0.04775448 Z3 0.10535586 Z4 0.05718296 x2 0.33901875 $dist [1] "binomial" > u <- uncertainty(o, type="asymp", B=9) > u <- uncertainty(o, type="boot", B=2) > > ## negbin > > summary(o <- multicut(Y ~ x2, strata=x0, dist="negbin")) Multivariate multticut results, dist = negbin Call: multicut.formula(formula = Y ~ x2, strata = x0, dist = "negbin") Species models with logLR >= -Inf: split assoc I logLR Spp2 2+3+4 + 0.1500 1.120 Spp1 1+2 +++ 0.4465 10.106 Spp3 2+3 + 0.4191 1.468 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 7: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 8: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 9: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 10: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached > summary(fitted(o)) Spp1 Spp2 Spp3 Min. :1.215 Min. :3.198 Min. :0.3200 1st Qu.:1.835 1st Qu.:4.104 1st Qu.:0.5528 Median :2.643 Median :4.244 Median :0.6718 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.824 3rd Qu.:4.433 3rd Qu.:0.8069 Max. :5.084 Max. :4.650 Max. :1.0131 > summary(predict(o)) Spp1 Spp2 Spp3 Min. :1.215 Min. :3.198 Min. :0.3200 1st Qu.:1.835 1st Qu.:4.104 1st Qu.:0.5528 Median :2.643 Median :4.244 Median :0.6718 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.824 3rd Qu.:4.433 3rd Qu.:0.8069 Max. :5.084 Max. :4.650 Max. :1.0131 > summary(predict(o, gnew=x0, xnew=data.frame(x2=x2))) Spp1 Spp2 Spp3 Min. :1.215 Min. :3.198 Min. :0.3200 1st Qu.:1.835 1st Qu.:4.104 1st Qu.:0.5528 Median :2.643 Median :4.244 Median :0.6718 Mean :2.780 Mean :4.120 Mean :0.6600 3rd Qu.:3.824 3rd Qu.:4.433 3rd Qu.:0.8069 Max. :5.084 Max. :4.650 Max. :1.0131 > getMLE(o, 1, vcov=FALSE) $coef `(Intercept)` Z2 Z3 Z4 x2 1.44494183 0.05330735 -0.90728222 -0.58512617 -0.20224456 $vcov NULL $dist [1] "negbin" > getMLE(o, 1, vcov=TRUE) $coef `(Intercept)` Z2 Z3 Z4 x2 1.44494183 0.05330735 -0.90728222 -0.58512617 -0.20224456 $vcov `(Intercept)` Z2 Z3 Z4 x2 `(Intercept)` 0.029340429 -0.025791781 -0.028089640 -0.027281015 -0.006125119 Z2 -0.025791781 0.046020182 0.026009776 0.026150709 -0.001067526 Z3 -0.028089640 0.026009776 0.065819767 0.026882617 0.003589934 Z4 -0.027281015 0.026150709 0.026882617 0.068294165 0.001950956 x2 -0.006125119 -0.001067526 0.003589934 0.001950956 0.012414812 $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 > 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) > > ## beta > > Y2 <- Y / rowSums(Y) > Y2[Y2 == 0] <- 0.01 > Y2[Y2 == 1] <- 0.99 > summary(o <- multicut(Y2 ~ x2, strata=x0, dist="beta")) Multivariate multticut results, dist = beta Call: multicut.formula(formula = Y2 ~ x2, strata = x0, dist = "beta") Species models with logLR >= -Inf: split assoc I logLR Spp1 1+2 +++ 0.5961 8.986 Spp2 3+4 ++ 0.5036 5.450 Spp3 3 ++ 0.3138 2.202 > summary(fitted(o)) Spp1 Spp2 Spp3 Min. :0.1530 Min. :0.3842 Min. :0.06575 1st Qu.:0.2228 1st Qu.:0.4912 1st Qu.:0.07768 Median :0.2944 Median :0.6153 Median :0.08715 Mean :0.3289 Mean :0.5862 Mean :0.09749 3rd Qu.:0.4420 3rd Qu.:0.6720 3rd Qu.:0.12083 Max. :0.5527 Max. :0.7390 Max. :0.15269 > summary(predict(o)) Spp1 Spp2 Spp3 Min. :0.1530 Min. :0.3842 Min. :0.06575 1st Qu.:0.2228 1st Qu.:0.4912 1st Qu.:0.07768 Median :0.2944 Median :0.6153 Median :0.08715 Mean :0.3289 Mean :0.5862 Mean :0.09749 3rd Qu.:0.4420 3rd Qu.:0.6720 3rd Qu.:0.12083 Max. :0.5527 Max. :0.7390 Max. :0.15269 > summary(predict(o, gnew=x0, xnew=data.frame(x2=x2))) Spp1 Spp2 Spp3 Min. :0.1530 Min. :0.3842 Min. :0.06575 1st Qu.:0.2228 1st Qu.:0.4912 1st Qu.:0.07768 Median :0.2944 Median :0.6153 Median :0.08715 Mean :0.3289 Mean :0.5862 Mean :0.09749 3rd Qu.:0.4420 3rd Qu.:0.6720 3rd Qu.:0.12083 Max. :0.5527 Max. :0.7390 Max. :0.15269 > getMLE(o, 1, vcov=FALSE) $coef `(Intercept)` Z2 Z3 Z4 x2 0.08830764 -0.18840568 -1.37412739 -1.10160840 -0.25095064 (phi) 5.00202161 $vcov NULL $dist [1] "beta" > getMLE(o, 1, vcov=TRUE) $coef `(Intercept)` Z2 Z3 Z4 x2 0.08830764 -0.18840568 -1.37412739 -1.10160840 -0.25095064 (phi) 5.00202161 $vcov `(Intercept)` Z2 Z3 Z4 x2 `(Intercept)` 0.074206595 -0.063932820 -0.070377695 -0.06829718 -0.014616010 Z2 -0.063932820 0.116988890 0.065144796 0.06543289 -0.002737528 Z3 -0.070377695 0.065144796 0.122121101 0.06988850 0.008374297 Z4 -0.068297184 0.065432889 0.069888504 0.13789290 0.004801650 x2 -0.014616010 -0.002737528 0.008374297 0.00480165 0.024749749 (phi) 0.004236217 -0.006876922 -0.058925018 -0.04555634 -0.010922723 (phi) `(Intercept)` 0.004236217 Z2 -0.006876922 Z3 -0.058925018 Z4 -0.045556339 x2 -0.010922723 (phi) 0.895238590 $dist [1] "beta" > u <- uncertainty(o, type="asymp", B=9) > u <- uncertainty(o, type="boot", B=2) > > ## zip > > Yzi <- Y > Yzi[1,] <- 0 > B <- sapply(2:3, function(i) which((1:n) != i)) # jackknife > B[1,] <- 1 > summary(o <- multicut(Yzi ~ x2, strata=x0, dist="zip")) Multivariate multticut results, dist = zip Call: multicut.formula(formula = Yzi ~ x2, strata = x0, dist = "zip") Species models with logLR >= -Inf: split assoc I logLR Spp2 2+3+4 + 0.1317 0.9127 Spp1 1+2 +++ 0.4446 9.3782 Spp3 2+3 + 0.4191 1.4684 > summary(fitted(o)) Spp1 Spp2 Spp3 Min. :1.228 Min. :3.161 Min. :0.3200 1st Qu.:1.827 1st Qu.:4.057 1st Qu.:0.5528 Median :2.669 Median :4.193 Median :0.6718 Mean :2.773 Mean :4.042 Mean :0.6600 3rd Qu.:3.821 3rd Qu.:4.291 3rd Qu.:0.8069 Max. :5.025 Max. :4.478 Max. :1.0131 > summary(predict(o)) Spp1 Spp2 Spp3 Min. :1.228 Min. :3.161 Min. :0.3200 1st Qu.:1.827 1st Qu.:4.057 1st Qu.:0.5528 Median :2.669 Median :4.193 Median :0.6718 Mean :2.773 Mean :4.042 Mean :0.6600 3rd Qu.:3.821 3rd Qu.:4.291 3rd Qu.:0.8069 Max. :5.025 Max. :4.478 Max. :1.0131 > summary(predict(o, gnew=x0, xnew=data.frame(x2=x2))) Spp1 Spp2 Spp3 Min. :1.228 Min. :3.161 Min. :0.3200 1st Qu.:1.827 1st Qu.:4.057 1st Qu.:0.5528 Median :2.669 Median :4.193 Median :0.6718 Mean :2.773 Mean :4.042 Mean :0.6600 3rd Qu.:3.821 3rd Qu.:4.291 3rd Qu.:0.8069 Max. :5.025 Max. :4.478 Max. :1.0131 > getMLE(o, 1, vcov=FALSE) $coef count_`(Intercept)` count_Z2 count_Z3 count_Z4 1.44115996 0.06588016 -0.88994736 -0.55167696 count_x2 zero_(Intercept) -0.19460908 -4.14805421 $vcov NULL $dist [1] "zip" > getMLE(o, 1, vcov=TRUE) $coef count_`(Intercept)` count_Z2 count_Z3 count_Z4 1.44115996 0.06588016 -0.88994736 -0.55167696 count_x2 zero_(Intercept) -0.19460908 -4.14805421 $vcov count_`(Intercept)` count_Z2 count_Z3 count_Z4 count_`(Intercept)` 0.029793908 -0.0266190177 -0.028853282 -0.028093104 count_Z2 -0.026619018 0.0479086757 0.028367776 0.030185899 count_Z3 -0.028853282 0.0283677758 0.070057343 0.032655602 count_Z4 -0.028093104 0.0301858990 0.032655602 0.082206224 count_x2 -0.006996918 0.0006100109 0.005104665 0.003575422 zero_(Intercept) -0.038789458 0.1310254474 0.201569457 0.371437924 count_x2 zero_(Intercept) count_`(Intercept)` -0.0069969182 -0.03878946 count_Z2 0.0006100109 0.13102545 count_Z3 0.0051046652 0.20156946 count_Z4 0.0035754217 0.37143792 count_x2 0.0140756555 0.07803250 zero_(Intercept) 0.0780325040 12.51325999 $dist [1] "zip" > u <- uncertainty(o, type="asymp", B=9) > u <- uncertainty(o, type="boot", B=B) > > ## zinb > > summary(o <- multicut(Yzi ~ x2, strata=x0, dist="zinb")) Multivariate multticut results, dist = zinb Call: multicut.formula(formula = Yzi ~ x2, strata = x0, dist = "zinb") Species models with logLR >= -Inf: split assoc I logLR Spp2 2+3+4 + 0.1317 0.9127 Spp1 1+2 +++ 0.4446 9.3676 Spp3 2+3 + 0.4191 1.4684 > summary(fitted(o)) Spp1 Spp2 Spp3 Min. :1.228 Min. :3.161 Min. :0.3200 1st Qu.:1.827 1st Qu.:4.057 1st Qu.:0.5528 Median :2.669 Median :4.193 Median :0.6718 Mean :2.773 Mean :4.042 Mean :0.6600 3rd Qu.:3.821 3rd Qu.:4.291 3rd Qu.:0.8069 Max. :5.025 Max. :4.478 Max. :1.0131 > summary(predict(o)) Spp1 Spp2 Spp3 Min. :1.228 Min. :3.161 Min. :0.3200 1st Qu.:1.827 1st Qu.:4.057 1st Qu.:0.5528 Median :2.669 Median :4.193 Median :0.6718 Mean :2.773 Mean :4.042 Mean :0.6600 3rd Qu.:3.821 3rd Qu.:4.291 3rd Qu.:0.8069 Max. :5.025 Max. :4.478 Max. :1.0131 > summary(predict(o, gnew=x0, xnew=data.frame(x2=x2))) Spp1 Spp2 Spp3 Min. :1.228 Min. :3.161 Min. :0.3200 1st Qu.:1.827 1st Qu.:4.057 1st Qu.:0.5528 Median :2.669 Median :4.193 Median :0.6718 Mean :2.773 Mean :4.042 Mean :0.6600 3rd Qu.:3.821 3rd Qu.:4.291 3rd Qu.:0.8069 Max. :5.025 Max. :4.478 Max. :1.0131 > getMLE(o, 1, vcov=FALSE) $coef count_`(Intercept)` count_Z2 count_Z3 count_Z4 1.44114046 0.06590496 -0.88994335 -0.55162038 count_x2 zero_(Intercept) -0.19460415 -4.14681253 $vcov NULL $dist [1] "zinb" > getMLE(o, 1, vcov=TRUE) $coef count_`(Intercept)` count_Z2 count_Z3 count_Z4 1.44114046 0.06590496 -0.88994335 -0.55162038 count_x2 zero_(Intercept) -0.19460415 -4.14681253 $vcov count_`(Intercept)` count_Z2 count_Z3 count_Z4 count_`(Intercept)` 0.029795698 -0.02662001 -0.028853791 -0.028091917 count_Z2 -0.026620008 0.04790858 0.028365353 0.030179561 count_Z3 -0.028853791 0.02836535 0.070056227 0.032646105 count_Z4 -0.028091917 0.03017956 0.032646105 0.082190717 count_x2 -0.006997295 0.00060887 0.005102488 0.003569842 zero_(Intercept) -0.038670691 0.13063442 0.201014898 0.370390990 count_x2 zero_(Intercept) count_`(Intercept)` -0.006997295 -0.03867069 count_Z2 0.000608870 0.13063442 count_Z3 0.005102488 0.20101490 count_Z4 0.003569842 0.37039099 count_x2 0.014076195 0.07779180 zero_(Intercept) 0.077791800 12.46631491 $dist [1] "zinb" > 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 > > ## zip2 > > summary(o <- multicut(Yzi ~ x2, strata=x0, dist="zip2")) Multivariate multticut results, dist = zip2 Call: multicut.formula(formula = Yzi ~ x2, strata = x0, dist = "zip2") Species models with logLR >= -Inf: split assoc I logLR Spp2 1+2+4 + 1.0000 0.5788 Spp1 1+2 + 1.0000 1.7272 Spp3 1+3 + 0.9999 0.4652 > summary(fitted(o)) Spp1 Spp2 Spp3 Min. :2.068 Min. :3.727 Min. :0.4573 1st Qu.:2.492 1st Qu.:3.913 1st Qu.:0.5971 Median :2.765 Median :4.056 Median :0.6699 Mean :2.789 Mean :4.040 Mean :0.6603 3rd Qu.:3.052 3rd Qu.:4.135 3rd Qu.:0.7297 Max. :3.593 Max. :4.301 Max. :0.8303 > summary(predict(o)) Spp1 Spp2 Spp3 Min. :2.068 Min. :3.727 Min. :0.4573 1st Qu.:2.492 1st Qu.:3.913 1st Qu.:0.5971 Median :2.765 Median :4.056 Median :0.6699 Mean :2.789 Mean :4.040 Mean :0.6603 3rd Qu.:3.052 3rd Qu.:4.135 3rd Qu.:0.7297 Max. :3.593 Max. :4.301 Max. :0.8303 > summary(predict(o, gnew=x0, xnew=data.frame(x2=x2))) Spp1 Spp2 Spp3 Min. :2.068 Min. :3.727 Min. :0.4573 1st Qu.:2.492 1st Qu.:3.913 1st Qu.:0.5971 Median :2.765 Median :4.056 Median :0.6699 Mean :2.789 Mean :4.040 Mean :0.6603 3rd Qu.:3.052 3rd Qu.:4.135 3rd Qu.:0.7297 Max. :3.593 Max. :4.301 Max. :0.8303 > getMLE(o, 1, vcov=FALSE) $coef ZZ(Intercept) ZZZ2 ZZZ3 ZZZ4 X(Intercept) 21.1660786 -8.7973889 -19.4350532 -19.4105296 1.1817714 Xx2 -0.1535313 $vcov NULL $dist [1] "zip2" > getMLE(o, 1, vcov=TRUE) $coef zero_ZZ(Intercept) zero_ZZZ2 zero_ZZZ3 zero_ZZZ4 -21.1660786 8.7973889 19.4350532 19.4105296 count_X(Intercept) count_Xx2 1.1817714 -0.1535313 $vcov zero_ZZ(Intercept) zero_ZZZ2 zero_ZZZ3 zero_ZZZ4 zero_ZZ(Intercept) 1.557080e+08 -1.557080e+08 -1.557080e+08 -1.557080e+08 zero_ZZZ2 -1.557080e+08 1.557936e+08 1.557080e+08 1.557080e+08 zero_ZZZ3 -1.557080e+08 1.557080e+08 1.557080e+08 1.557080e+08 zero_ZZZ4 -1.557080e+08 1.557080e+08 1.557080e+08 1.557080e+08 count_X(Intercept) -1.097392e-08 -1.933296e-02 8.190644e-03 8.488029e-03 count_Xx2 -3.142452e-08 1.866185e-01 -2.101927e-03 -3.513862e-03 count_X(Intercept) count_Xx2 zero_ZZ(Intercept) -1.094986e-08 -3.151499e-08 zero_ZZZ2 -1.933296e-02 1.866185e-01 zero_ZZZ3 8.190644e-03 -2.101926e-03 zero_ZZZ4 8.488029e-03 -3.513862e-03 count_X(Intercept) 1.041349e-02 -5.758719e-03 count_Xx2 -5.758719e-03 1.242543e-02 $dist [1] "zip2" > u <- uncertainty(o, type="asymp", B=9) > u <- uncertainty(o, type="boot", B=B) > > ## zinb2 > > summary(o <- multicut(Yzi ~ x2, strata=x0, dist="zinb2")) Multivariate multticut results, dist = zinb2 Call: multicut.formula(formula = Yzi ~ x2, strata = x0, dist = "zinb2") Species models with logLR >= -Inf: split assoc I logLR Spp2 1+2+4 + 1.0000 0.5788 Spp1 1+2 + 1.0000 1.7167 Spp3 1+3 + 0.9999 0.4652 Warning message: In sqrt(diag(vc)[np]) : NaNs produced > summary(fitted(o)) Spp1 Spp2 Spp3 Min. :2.068 Min. :3.727 Min. :0.4573 1st Qu.:2.492 1st Qu.:3.913 1st Qu.:0.5971 Median :2.765 Median :4.056 Median :0.6699 Mean :2.789 Mean :4.040 Mean :0.6603 3rd Qu.:3.052 3rd Qu.:4.135 3rd Qu.:0.7297 Max. :3.593 Max. :4.301 Max. :0.8303 > summary(predict(o)) Spp1 Spp2 Spp3 Min. :2.068 Min. :3.727 Min. :0.4573 1st Qu.:2.492 1st Qu.:3.913 1st Qu.:0.5971 Median :2.765 Median :4.056 Median :0.6699 Mean :2.789 Mean :4.040 Mean :0.6603 3rd Qu.:3.052 3rd Qu.:4.135 3rd Qu.:0.7297 Max. :3.593 Max. :4.301 Max. :0.8303 > summary(predict(o, gnew=x0, xnew=data.frame(x2=x2))) Spp1 Spp2 Spp3 Min. :2.068 Min. :3.727 Min. :0.4573 1st Qu.:2.492 1st Qu.:3.913 1st Qu.:0.5971 Median :2.765 Median :4.056 Median :0.6699 Mean :2.789 Mean :4.040 Mean :0.6603 3rd Qu.:3.052 3rd Qu.:4.135 3rd Qu.:0.7297 Max. :3.593 Max. :4.301 Max. :0.8303 > getMLE(o, 1, vcov=FALSE) $coef ZZ(Intercept) ZZZ2 ZZZ3 ZZZ4 X(Intercept) 21.2218745 -8.6301088 -19.4908254 -19.4663004 1.1817691 Xx2 -0.1535308 $vcov NULL $dist [1] "zinb2" > getMLE(o, 1, vcov=TRUE) $coef zero_ZZ(Intercept) zero_ZZZ2 zero_ZZZ3 zero_ZZZ4 -21.2218745 8.6301088 19.4908254 19.4663004 count_X(Intercept) count_Xx2 1.1817691 -0.1535308 $vcov zero_ZZ(Intercept) zero_ZZZ2 zero_ZZZ3 zero_ZZZ4 zero_ZZ(Intercept) 1.646429e+08 -1.646429e+08 -1.646429e+08 -1.646429e+08 zero_ZZZ2 -1.646429e+08 1.647498e+08 1.646429e+08 1.646429e+08 zero_ZZZ3 -1.646429e+08 1.646429e+08 1.646429e+08 1.646429e+08 zero_ZZZ4 -1.646429e+08 1.646429e+08 1.646429e+08 1.646429e+08 count_X(Intercept) -1.956221e-07 -1.931320e-02 8.194367e-03 8.491784e-03 count_Xx2 -8.984079e-08 1.865965e-01 -2.102647e-03 -3.514764e-03 count_X(Intercept) count_Xx2 zero_ZZ(Intercept) -1.964344e-07 -8.958747e-08 zero_ZZZ2 -1.931319e-02 1.865965e-01 zero_ZZZ3 8.194368e-03 -2.102647e-03 zero_ZZZ4 8.491784e-03 -3.514764e-03 count_X(Intercept) 1.041422e-02 -5.759057e-03 count_Xx2 -5.759057e-03 1.242586e-02 $dist [1] "zinb2" > u <- uncertainty(o, type="asymp", B=9) Warning message: In sqrt(diag(vc)[np]) : NaNs produced > u <- uncertainty(o, type="boot", B=B) Warning message: In sqrt(diag(vc)[np]) : NaNs produced > > ## 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 <- multicut(Y ~ x1 + x2, dd, strata=dd$x0, dist="rsf") > o$species $`Sp 1` Univariate multicut results, dist = rsf logLR = 11.52 (logL_null = -6896) Expected values: 1 2 3 1.0000 0.9344 1.3148 > summary(fitted(o)) Sp 1 Min. :0.7265 1st Qu.:1.0471 Median :1.1868 Mean :1.2537 3rd Qu.:1.4362 Max. :2.3192 > summary(predict(o)) Sp 1 Min. :0.7265 1st Qu.:1.0471 Median :1.1868 Mean :1.2537 3rd Qu.:1.4362 Max. :2.3192 > summary(predict(o, gnew=o$strata, xnew=o$X)) Sp 1 Min. :0.7265 1st Qu.:1.0471 Median :1.1868 Mean :1.2537 3rd Qu.:1.4362 Max. :2.3192 > getMLE(o, 1, vcov=FALSE) $coef (Intercept) Z2 Z3 x1 x2 0.00000000 -0.06780372 0.27366076 -0.15097950 0.21535663 $vcov NULL $dist [1] "rsf" > getMLE(o, 1, vcov=TRUE) $coef (Intercept) Z2 Z3 x1 x2 0.00000000 -0.06780372 0.27366076 -0.15097950 0.21535663 $vcov (Intercept) Z2 Z3 x1 x2 (Intercept) NA NA NA NA NA Z2 NA 6.719089e-03 0.0033684705 0.0001488806 2.627768e-05 Z3 NA 3.368470e-03 0.0058581108 0.0001290955 -3.249933e-04 x1 NA 1.488806e-04 0.0001290955 0.0010537759 -2.226946e-04 x2 NA 2.627768e-05 -0.0003249933 -0.0002226946 1.137400e-02 $dist [1] "rsf" > u <- uncertainty(o, type="asymp", B=9) > u <- uncertainty(o, type="boot", B=2) > > ## intercept + partition > o <- multicut(Y ~ 1, dd, strata=x0, dist="rsf") > o$species $`Sp 1` Univariate multicut results, dist = rsf logLR = 12.14 (logL_null = -6909) Expected values: 1 2 3 1.0000 0.9538 1.3441 > summary(fitted(o)) Sp 1 Min. :0.9538 1st Qu.:0.9538 Median :1.0000 Mean :1.1111 3rd Qu.:1.3441 Max. :1.3441 > summary(predict(o)) Sp 1 Min. :0.9538 1st Qu.:0.9538 Median :1.0000 Mean :1.1111 3rd Qu.:1.3441 Max. :1.3441 > summary(predict(o, gnew=o$strata, xnew=NULL)) Sp 1 Min. :0.9538 1st Qu.:0.9538 Median :1.0000 Mean :1.1111 3rd Qu.:1.3441 Max. :1.3441 > getMLE(o, 1, vcov=FALSE) $coef (Intercept) Z2 Z3 0.0000000 -0.0473314 0.2957455 $vcov NULL $dist [1] "rsf" > getMLE(o, 1, vcov=TRUE) $coef (Intercept) Z2 Z3 0.0000000 -0.0473314 0.2957455 $vcov (Intercept) Z2 Z3 (Intercept) NA NA NA Z2 NA 0.006701094 0.003356089 Z3 NA 0.003356089 0.005836988 $dist [1] "rsf" > u <- uncertainty(o, type="asymp", B=9) > u <- uncertainty(o, type="boot", B=2) > > ## rspf > > o <- multicut(Y ~ x1 + x2, dd, strata=x0, dist="rspf") > o$species $`Sp 1` Univariate multicut results, dist = rspf logLR = 4.984 (logL_null = -6889) Expected values: 1 2 3 0.3585 0.3332 0.4808 > summary(fitted(o)) Sp 1 Min. :0.2310 1st Qu.:0.3783 Median :0.4357 Mean :0.4507 3rd Qu.:0.5245 Max. :0.7376 > summary(predict(o)) Sp 1 Min. :0.2310 1st Qu.:0.3783 Median :0.4357 Mean :0.4507 3rd Qu.:0.5245 Max. :0.7376 > summary(predict(o, gnew=o$strata, xnew=o$X)) Sp 1 Min. :0.2310 1st Qu.:0.3783 Median :0.4357 Mean :0.4507 3rd Qu.:0.5245 Max. :0.7376 > getMLE(o, 1, vcov=FALSE) $coef (Intercept) Z2 Z3 x1 x2 -0.5816923 -0.1118811 0.5048065 -0.3034710 0.3949233 $vcov NULL $dist [1] "rspf" > getMLE(o, 1, vcov=TRUE) $coef (Intercept) Z2 Z3 x1 x2 -0.5816923 -0.1118811 0.5048065 -0.3034710 0.3949233 $vcov (Intercept) Z2 Z3 x1 x2 (Intercept) 1.38134824 -0.074928463 0.331959685 -0.24508261 0.24782213 Z2 -0.07492846 0.021925177 -0.006565589 0.01226921 -0.01188509 Z3 0.33195968 -0.006565589 0.107004907 -0.06109673 0.06723105 x1 -0.24508261 0.012269212 -0.061096734 0.04844722 -0.04903877 x2 0.24782213 -0.011885085 0.067231053 -0.04903877 0.09285592 $dist [1] "rspf" > u <- uncertainty(o, type="asymp", B=9) > u <- uncertainty(o, type="boot", B=2) > > ## --- testing ... passing --- > > 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 <- multicut(Y ~ x2, strata=x0, dist="poisson") > ## with offsets: log Area > wo <- multicut(Y ~ x2, strata=x0, dist="poisson", + offset=log(A), weights=rep(1,n)) > agg <- aggregate(data.frame(lam=lam), list(x0=x0), mean) > agg$wout_off <- no$species[[1]]$mu > agg$with_off <- wo$species[[1]]$mu > agg x0 lam wout_off with_off 1 1 4.359026 4.116151 4.116151 2 2 4.413488 48.743325 4.874332 3 3 1.541015 1.970638 1.970638 4 4 1.603424 15.813358 1.581336 > > > proc.time() user system elapsed 7.82 0.40 8.23