Running 'tests.R' [269s]
Running the tests in 'tests/tests.R' failed.
Complete output:
  > #devtools::install_github("psolymos/ResourceSelection")
  > library(ResourceSelection)
  ResourceSelection 0.3-6 	 2023-06-27
  > 
  > ## --- run examples with \dontrun sections ---
  > 
  > help_pages <- c("CAIC", "goats", "hoslem.test", "kdepairs",
  +     "makeUsedAvail", "rsf", "simulateUsedAvail", "mep")
  > 
  > for (i in help_pages) {
  +     cat("\n\n---------- ResourceSelection example:", i, "----------\n\n")
  +     eval(parse(text=paste0("example(", i,
  +         ", package = 'ResourceSelection', run.dontrun = TRUE)")))
  + }
  
  
  ---------- ResourceSelection example: CAIC ----------
  
  
  CAIC> ## compare some random models
  CAIC> y <- rnorm(10)
  
  CAIC> a <- lm(y ~ runif(10))
  
  CAIC> b <- lm(y ~ runif(10))
  
  CAIC> 0.5*(AIC(a) + BIC(a))
  [1] 37.39286
  
  CAIC> CAIC(a)
  [1] 37.39286
  
  CAIC> AIC(a)
  [1] 36.93898
  
  CAIC> CAIC(a, alpha=1)
  [1] 36.93898
  
  CAIC> BIC(a)
  [1] 37.84674
  
  CAIC> CAIC(a, alpha=0)
  [1] 37.84674
  
  CAIC> CAIC(a, b)
    df     CAIC
  a  3 37.39286
  b  3 37.78892
  
  CAIC> CAIC(a, b, alpha=0.2)
    df     CAIC
  a  3 37.66519
  b  3 38.06124
  
  CAIC> CAICtable(a, b, alpha=1)
    df     CAIC     dCAIC     wCAIC
  a  3 36.93898 0.0000000 0.5493459
  b  3 37.33504 0.3960561 0.4506541
  
  CAIC> ## you can use global option
  CAIC> ## useful when inside of xv or bootstrap
  CAIC> ## no need for extra argument
  CAIC> getOption("CAIC_alpha")
  NULL
  
  CAIC> op <- options(CAIC_alpha = 0.2)
  
  CAIC> getOption("CAIC_alpha")
  [1] 0.2
  
  CAIC> CAIC(a,b)
    df     CAIC
  a  3 37.66519
  b  3 38.06124
  
  CAIC> options(op)
  
  CAIC> getOption("CAIC_alpha")
  NULL
  
  
  ---------- ResourceSelection example: goats ----------
  
  
  goats> data(goats)
  
  goats> str(goats)
  'data.frame':	19014 obs. of  8 variables:
   $ STATUS   : int  1 1 1 1 1 1 1 1 1 1 ...
   $ ID       : int  1 1 1 1 1 1 1 1 1 1 ...
   $ ELEVATION: int  651 660 316 334 454 343 429 493 400 442 ...
   $ SLOPE    : num  38.5 39.7 20.5 34.1 41.6 ...
   $ ET       : num  35.4 70.7 50 35.4 25 ...
   $ ASPECT   : num  243 270 279 266 258 ...
   $ HLI      : num  0.918 0.884 0.713 0.864 0.935 ...
   $ TASP     : num  0.947 0.699 0.575 0.745 0.829 ...
  
  goats> summary(goats)
       STATUS             ID           ELEVATION        SLOPE      
   Min.   :0.0000   Min.   : 1.000   Min.   : 248   Min.   : 0.00  
   1st Qu.:0.0000   1st Qu.: 3.000   1st Qu.: 685   1st Qu.:25.41  
   Median :0.0000   Median : 6.000   Median : 985   Median :34.18  
   Mean   :0.3333   Mean   : 5.785   Mean   :1000   Mean   :33.65  
   3rd Qu.:1.0000   3rd Qu.: 8.000   3rd Qu.:1322   3rd Qu.:42.68  
   Max.   :1.0000   Max.   :10.000   Max.   :1966   Max.   :70.80  
         ET            ASPECT           HLI               TASP         
   Min.   :  0.0   Min.   : -1.0   Min.   :0.07749   Min.   :-1.00000  
   1st Qu.: 25.0   1st Qu.:139.5   1st Qu.:0.53882   1st Qu.:-0.04676  
   Median : 75.0   Median :178.2   Median :0.73698   Median : 0.53151  
   Mean   :126.8   Mean   :180.7   Mean   :0.68497   Mean   : 0.33343  
   3rd Qu.:180.3   3rd Qu.:229.6   3rd Qu.:0.86565   3rd Qu.: 0.86441  
   Max.   :995.3   Max.   :359.6   Max.   :0.98321   Max.   : 1.00000  
  
  goats> goats$exp.HLI <- exp(goats$HLI)
  
  goats> goats$sin.SLOPE <- sin(pi * goats$SLOPE / 180)
  
  goats> goats$ELEVATION <- scale(goats$ELEVATION)
  
  goats> goats$ET <- scale(goats$ET)
  
  goats> goats$TASP <- scale(goats$TASP)
  
  goats> m1 <- rspf(STATUS ~ TASP + sin.SLOPE + ELEVATION, goats, m=0, B = 99)
  Loading required namespace: pbapply
  
  goats> m2 <- rspf(STATUS ~ TASP + ELEVATION, goats, m=0, B = 99)
  
  goats> summary(m1)
  
  Call:
  rspf(formula = STATUS ~ TASP + sin.SLOPE + ELEVATION, data = goats, m = 0, 
      B = 99)
  
  Resource Selection Probability Function (Logistic RSPF) model
  Non-matched Used-Available design
  Maximum Likelihood estimates
  with Nonparametric Bootstrap standard errors (B = 99)
  
  Fitted probabilities:
       Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
  1.947e-08 4.280e-07 9.977e-07 1.376e-06 1.924e-06 8.793e-06 
  
  Coefficients (logit link):
                Estimate Std. Error z value Pr(>|z|)    
  (Intercept) -16.894536   0.288868  -58.48   <2e-16 ***
  TASP          0.391159   0.014587   26.82   <2e-16 ***
  sin.SLOPE     5.366405   0.087530   61.31   <2e-16 ***
  ELEVATION     0.098295   0.009628   10.21   <2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
  
  Log-likelihood: -5.729e+04 
  BIC = 1.146e+05 
  
  Hosmer and Lemeshow goodness of fit (GOF) test:
  X-squared = 152.4, df = 8, p-value < 2.2e-16
  
  
  goats> summary(m2)
  
  Call:
  rspf(formula = STATUS ~ TASP + ELEVATION, data = goats, m = 0, B = 99)
  
  Resource Selection Probability Function (Logistic RSPF) model
  Non-matched Used-Available design
  Maximum Likelihood estimates
  with Nonparametric Bootstrap standard errors (B = 99)
  
  Fitted probabilities:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.01194 0.58006 0.86178 0.73656 0.95712 0.99835 
  
  Coefficients (logit link):
              Estimate Std. Error z value Pr(>|z|)    
  (Intercept)  1.62906    0.10367   15.71   <2e-16 ***
  TASP         1.86071    0.07501   24.81   <2e-16 ***
  ELEVATION    1.14338    0.08837   12.94   <2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
  
  Log-likelihood: -5.91e+04 
  BIC = 1.182e+05 
  
  Hosmer and Lemeshow goodness of fit (GOF) test:
  X-squared = 174.3, df = 8, p-value < 2.2e-16
  
  
  goats> AIC(m1, m2)
     df      AIC
  m1  4 114578.2
  m2  3 118215.1
  
  goats> plot(m1)
  
  
  ---------- ResourceSelection example: hoslem.test ----------
  
  
  hslm.t> set.seed(123)
  
  hslm.t> n <- 500
  
  hslm.t> x <- rnorm(n)
  
  hslm.t> y <- rbinom(n, 1, plogis(0.1 + 0.5*x))
  
  hslm.t> m <- glm(y ~ x, family=binomial)
  
  hslm.t> hoslem.test(m$y, fitted(m))
  
  	Hosmer and Lemeshow goodness of fit (GOF) test
  
  data:  m$y, fitted(m)
  X-squared = 4.5227, df = 8, p-value = 0.8072
  
  
  
  ---------- ResourceSelection example: kdepairs ----------
  
  
  kdeprs> kdepairs(iris[1:4])
  
  
  ---------- ResourceSelection example: makeUsedAvail ----------
  
  
  mkUsdA> (x <- data.frame(species=rep(1:0,each=4), var1=1:8, var2=11:18))
    species var1 var2
  1       1    1   11
  2       1    2   12
  3       1    3   13
  4       1    4   14
  5       0    5   15
  6       0    6   16
  7       0    7   17
  8       0    8   18
  
  mkUsdA> makeUsedAvail(species ~ var1 + var2, x)
     species var1 var2
  1        1    1   11
  2        1    2   12
  3        1    3   13
  4        1    4   14
  5        0    1   11
  6        0    2   12
  7        0    3   13
  8        0    4   14
  9        0    5   15
  10       0    6   16
  11       0    7   17
  12       0    8   18
  
  
  ---------- ResourceSelection example: rsf ----------
  
  
  rsf> ## --- Simulated data example ---
  rsf> 
  rsf> ## settings
  rsf> n.used <- 1000
  
  rsf> m <- 10
  
  rsf> n <- n.used * m
  
  rsf> set.seed(1234)
  
  rsf> x <- data.frame(x1=rnorm(n), x2=runif(n))
  
  rsf> cfs <- c(1.5,-1,0.5)
  
  rsf> ## fitting Exponential RSF model
  rsf> dat1 <- simulateUsedAvail(x, cfs, n.used, m, link="log")
  
  rsf> m1 <- rsf(status ~ .-status, dat1, m=0, B=0)
  
  rsf> summary(m1)
  
  Call:
  rsf(formula = status ~ . - status, data = dat1, m = 0, B = 0)
  
  Resource Selection Function (Exponential RSF) model
  Non-matched Used-Available design
  Maximum Likelihood estimates
  
  Fitted values:
      Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
   0.03444  0.67639  1.34572  2.34955  2.75012 40.08155 
  
  Coefficients (log link):
     Estimate Std. Error z value Pr(>|z|)
  x1  -1.0036         NA      NA       NA
  x2   0.4506         NA      NA       NA
  
  Log-likelihood: -8714 
  BIC = 1.744e+04 
  
  Hosmer and Lemeshow goodness of fit (GOF) test:
  X-squared = 6.293, df = 8, p-value 0.6145
  
  
  rsf> ## fitting Logistic RSPF model
  rsf> dat2 <- simulateUsedAvail(x, cfs, n.used, m, link="logit")
  
  rsf> m2 <- rspf(status ~ .-status, dat2, m=0, B=0)
  
  rsf> summary(m2)
  
  Call:
  rspf(formula = status ~ . - status, data = dat2, m = 0, B = 0)
  
  Resource Selection Probability Function (Logistic RSPF) model
  Non-matched Used-Available design
  Maximum Likelihood estimates
  
  Fitted probabilities:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.1158  0.8306  0.9151  0.8712  0.9604  0.9986 
  
  Coefficients (logit link):
              Estimate Std. Error z value Pr(>|z|)  
  (Intercept)   1.7294     0.8221   2.104   0.0354 *
  x1           -1.1667     0.5352  -2.180   0.0293 *
  x2            1.2669     1.0942   1.158   0.2469  
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
  
  Log-likelihood: -9197 
  BIC = 1.841e+04 
  
  Hosmer and Lemeshow goodness of fit (GOF) test:
  X-squared = 4.256, df = 8, p-value 0.8334
  
  
  rsf> ## --- Real data analysis from Lele & Keim 2006 ---
  rsf> 
  rsf> 
  rsf> goats$exp.HLI <- exp(goats$HLI)
  
  rsf> goats$sin.SLOPE <- sin(pi * goats$SLOPE / 180)
  
  rsf> goats$ELEVATION <- scale(goats$ELEVATION)
  
  rsf> goats$ET <- scale(goats$ET)
  
  rsf> goats$TASP <- scale(goats$TASP)
  
  rsf> ## Fit two RSPF models:
  rsf> ## global availability (m=0) and bootstrap (B=99)
  rsf> m1 <- rspf(STATUS ~ TASP + sin.SLOPE + ELEVATION, goats, m=0, B = 99)
  
  rsf> m2 <- rspf(STATUS ~ TASP + ELEVATION, goats, m=0, B = 99)
  
  rsf> ## Inspect the summaries
  rsf> summary(m1)
  
  Call:
  rspf(formula = STATUS ~ TASP + sin.SLOPE + ELEVATION, data = goats, m = 0, 
      B = 99)
  
  Resource Selection Probability Function (Logistic RSPF) model
  Non-matched Used-Available design
  Maximum Likelihood estimates
  with Nonparametric Bootstrap standard errors (B = 99)
  
  Fitted probabilities:
       Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
  1.947e-08 4.280e-07 9.977e-07 1.376e-06 1.924e-06 8.793e-06 
  
  Coefficients (logit link):
               Estimate Std. Error z value Pr(>|z|)    
  (Intercept) -16.89454    0.27100 -62.342  < 2e-16 ***
  TASP          0.39116    0.01290  30.315  < 2e-16 ***
  sin.SLOPE     5.36640    0.09501  56.480  < 2e-16 ***
  ELEVATION     0.09829    0.01354   7.258 3.93e-13 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
  
  Log-likelihood: -5.729e+04 
  BIC = 1.146e+05 
  
  Hosmer and Lemeshow goodness of fit (GOF) test:
  X-squared = 152.4, df = 8, p-value < 2.2e-16
  
  
  rsf> summary(m2)
  
  Call:
  rspf(formula = STATUS ~ TASP + ELEVATION, data = goats, m = 0, B = 99)
  
  Resource Selection Probability Function (Logistic RSPF) model
  Non-matched Used-Available design
  Maximum Likelihood estimates
  with Nonparametric Bootstrap standard errors (B = 99)
  
  Fitted probabilities:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.01194 0.58006 0.86178 0.73656 0.95712 0.99835 
  
  Coefficients (logit link):
              Estimate Std. Error z value Pr(>|z|)    
  (Intercept)  1.62906    0.09417   17.30   <2e-16 ***
  TASP         1.86071    0.06706   27.75   <2e-16 ***
  ELEVATION    1.14338    0.07850   14.56   <2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
  
  Log-likelihood: -5.91e+04 
  BIC = 1.182e+05 
  
  Hosmer and Lemeshow goodness of fit (GOF) test:
  X-squared = 174.3, df = 8, p-value < 2.2e-16
  
  
  rsf> ## Compare models: looks like m1 is better supported
  rsf> CAIC(m1, m2)
     df     CAIC
  m1  4 114591.7
  m2  3 118225.2
  
  rsf> ## Visualize the relationships
  rsf> plot(m1)
  
  rsf> mep(m1) # marginal effects similar to plot but with CIs
  
  rsf> kdepairs(m1) # 2D kernel density estimates
  
  rsf> plot(m2)
  
  rsf> kdepairs(m2)
  
  rsf> mep(m2)
  
  rsf> ## fit and compare to null RSF model (not available for RSPF)
  rsf> m3 <- rsf(STATUS ~ TASP + ELEVATION, goats, m=0, B = 0)
  
  rsf> m4 <- rsf.null(Y=goats$STATUS, m=0)
  
  rsf> CAIC(m3, m4)
     df     CAIC
  m3  2 118604.5
  m4  0 119757.1
  
  
  ---------- ResourceSelection example: simulateUsedAvail ----------
  
  
  smltUA> n.used <- 1000
  
  smltUA> m <- 10
  
  smltUA> n <- n.used * m
  
  smltUA> set.seed(1234)
  
  smltUA> x <- data.frame(x1=rnorm(n), x2=runif(n))
  
  smltUA> cfs <- c(1.5,-1,0.5)
  
  smltUA> dat1 <- simulateUsedAvail(x, cfs, n.used, m, link="log")
  
  smltUA> str(dat1)
  'data.frame':	11000 obs. of  3 variables:
   $ status: num  1 1 1 1 1 1 1 1 1 1 ...
   $ x1    : num  -0.19 -2.023 0.422 -0.982 -0.569 ...
   $ x2    : num  0.995 0.542 0.849 0.521 0.531 ...
  
  smltUA> dat2 <- simulateUsedAvail(x, cfs, n.used, m, link="logit")
  
  smltUA> str(dat2)
  'data.frame':	11000 obs. of  3 variables:
   $ status: num  1 1 1 1 1 1 1 1 1 1 ...
   $ x1    : num  0.1955 0.3459 -0.0139 -0.2431 -0.9787 ...
   $ x2    : num  0.0527 0.4838 0.553 0.2334 0.9344 ...
  
  
  ---------- ResourceSelection example: mep ----------
  
  
  mep> data(goats)
  
  mep> goats$ELEVATION <- goats$ELEVATION/1000
  
  mep> goats$TASPc <- cut(goats$TASP, 3, ordered_result=FALSE)
  
  mep> goats$SLOPEc <- cut(goats$SLOPE, 3, ordered_result=TRUE)
  
  mep> fit <- rspf(STATUS ~ TASPc + SLOPEc + ELEVATION + I(ELEVATION^2), goats, m=0, B=0)
  
  mep> op <- par(mfrow=c(2,2))
  
  mep> mep(fit, which=1:4)#, subset=sample.int(nrow(goats), 10^4))
  
  mep> par(op)
  There were 11 warnings (use warnings() to see them)
  > 
  > ## --- test methods ---
  > 
  > ## settings
  > n.used <- 1000
  > m <- 10
  > n <- n.used * m
  > set.seed(1234)
  > x <- data.frame(x1=rnorm(n), x2=runif(n))
  > cfs <- c(1.5,-1,0.5)
  > dat1 <- simulateUsedAvail(x, cfs, n.used, m, link="log")
  > m0 <- rsf.null(dat1$status, m=0)
  > m1 <- rsf(status ~ .-status, dat1, m=0, B=0)
  > m1b <- rsf(status ~ .-status, dat1, m=0, B=9)
  > dat2 <- simulateUsedAvail(x, cfs, n.used, m, link="logit")
  > m2 <- rspf(status ~ .-status, dat2, m=0, B=0)
  > m2b <- rspf(status ~ .-status, dat2, m=0, B=9)
  > 
  > print(m0)
  
  Call:
  rsf.null(Y = dat1$status, m = 0)
  
  Resource Selection Function (Exponential RSF) model
  Non-matched Used-Available design
  Maximum Likelihood estimates
  
  Coefficients (log link):
  character(0)
  
  > print(m1)
  
  Call:
  rsf(formula = status ~ . - status, data = dat1, m = 0, B = 0)
  
  Resource Selection Function (Exponential RSF) model
  Non-matched Used-Available design
  Maximum Likelihood estimates
  
  Coefficients (log link):
       x1       x2  
  -1.0036   0.4506  
  
  > print(m1b)
  
  Call:
  rsf(formula = status ~ . - status, data = dat1, m = 0, B = 9)
  
  Resource Selection Function (Exponential RSF) model
  Non-matched Used-Available design
  Maximum Likelihood estimates
  
  Coefficients (log link):
       x1       x2  
  -1.0036   0.4506  
  
  > print(m2)
  
  Call:
  rspf(formula = status ~ . - status, data = dat2, m = 0, B = 0)
  
  Resource Selection Probability Function (Logistic RSPF) model
  Non-matched Used-Available design
  Maximum Likelihood estimates
  
  Coefficients (logit link):
  (Intercept)           x1           x2  
       1.0311      -0.8973       0.1406  
  
  > print(m2b)
  
  Call:
  rspf(formula = status ~ . - status, data = dat2, m = 0, B = 9)
  
  Resource Selection Probability Function (Logistic RSPF) model
  Non-matched Used-Available design
  Maximum Likelihood estimates
  
  Coefficients (logit link):
  (Intercept)           x1           x2  
       1.0311      -0.8973       0.1406  
  
  > 
  > summary(m0)
  Error in cut.default(yhat, breaks = qq, include.lowest = TRUE) : 
    invalid number of intervals
  
  Call:
  rsf.null(Y = dat1$status, m = 0)
  
  Resource Selection Function (Exponential RSF) model
  Non-matched Used-Available design
  Maximum Likelihood estimates
  
  Fitted values:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
        1       1       1       1       1       1 
  
  Coefficients (log link):
       Estimate Std. Error z value Pr(>|z|)
  
  Log-likelihood: -9210 
  BIC = 1.842e+04 
  
  Warning message:
  In summary.rsf(m0) : Hosmer and Lemeshow test cannot be calculated
  > summary(m1)
  
  Call:
  rsf(formula = status ~ . - status, data = dat1, m = 0, B = 0)
  
  Resource Selection Function (Exponential RSF) model
  Non-matched Used-Available design
  Maximum Likelihood estimates
  
  Fitted values:
      Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
   0.03444  0.67639  1.34572  2.34955  2.75012 40.08155 
  
  Coefficients (log link):
     Estimate Std. Error z value Pr(>|z|)
  x1  -1.0036         NA      NA       NA
  x2   0.4506         NA      NA       NA
  
  Log-likelihood: -8714 
  BIC = 1.744e+04 
  
  Hosmer and Lemeshow goodness of fit (GOF) test:
  X-squared = 6.293, df = 8, p-value 0.6145
  
  > summary(m1b)
  
  Call:
  rsf(formula = status ~ . - status, data = dat1, m = 0, B = 9)
  
  Resource Selection Function (Exponential RSF) model
  Non-matched Used-Available design
  Maximum Likelihood estimates
  with Nonparametric Bootstrap standard errors (B = 9)
  
  Fitted values:
      Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
   0.03444  0.67639  1.34572  2.34955  2.75012 40.08155 
  
  Coefficients (log link):
     Estimate Std. Error z value Pr(>|z|)    
  x1 -1.00362    0.02521 -39.812  < 2e-16 ***
  x2  0.45056    0.12857   3.504 0.000458 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
  
  Log-likelihood: -8714 
  BIC = 1.744e+04 
  
  Hosmer and Lemeshow goodness of fit (GOF) test:
  X-squared = 6.293, df = 8, p-value 0.6145
  
  > summary(m2)
  
  Call:
  rspf(formula = status ~ . - status, data = dat2, m = 0, B = 0)
  
  Resource Selection Probability Function (Logistic RSPF) model
  Non-matched Used-Available design
  Maximum Likelihood estimates
  
  Fitted probabilities:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.1065  0.6273  0.7538  0.7238  0.8474  0.9842 
  
  Coefficients (logit link):
              Estimate Std. Error z value Pr(>|z|)   
  (Intercept)   1.0311     0.6954   1.483  0.13814   
  x1           -0.8973     0.3095  -2.899  0.00374 **
  x2            0.1406     0.3868   0.363  0.71624   
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
  
  Log-likelihood: -9182 
  BIC = 1.839e+04 
  
  Hosmer and Lemeshow goodness of fit (GOF) test:
  X-squared = 4.857, df = 8, p-value 0.7727
  
  > summary(m2b)
  
  Call:
  rspf(formula = status ~ . - status, data = dat2, m = 0, B = 9)
  
  Resource Selection Probability Function (Logistic RSPF) model
  Non-matched Used-Available design
  Maximum Likelihood estimates
  with Nonparametric Bootstrap standard errors (B = 9)
  
  Fitted probabilities:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.1065  0.6273  0.7538  0.7238  0.8474  0.9842 
  
  Coefficients (logit link):
              Estimate Std. Error z value Pr(>|z|)
  (Intercept)   1.0311     2.3533   0.438    0.661
  x1           -0.8973     0.8673  -1.035    0.301
  x2            0.1406     1.3340   0.105    0.916
  
  Log-likelihood: -9182 
  BIC = 1.839e+04 
  
  Hosmer and Lemeshow goodness of fit (GOF) test:
  X-squared = 4.857, df = 8, p-value 0.7727
  
  > 
  > CAIC(m0, m1, m2)
     df     CAIC
  m0  0 18420.88
  m1  2 17437.37
  m2  3 18378.31
  > 
  > family(m1)
  
  Family: binomial 
  Link function: log 
  
  > family(m2)
  
  Family: binomial 
  Link function: logit 
  
  > 
  > coef(m0)
  numeric(0)
  > coef(m1)
          x1         x2 
  -1.0036174  0.4505615 
  > coef(m2)
  (Intercept)          x1          x2 
    1.0311420  -0.8973184   0.1405936 
  > 
  > vcov(m0)
  <0 x 0 matrix>
  > vcov(m1)
               x1           x2
  x1 1.060573e-03 8.629123e-05
  x2 8.629123e-05 1.237935e-02
  > vcov(m1b)
               x1          x2
  x1 0.0006355037 0.001940086
  x2 0.0019400857 0.016531321
  > vcov(m2)
              (Intercept)           x1           x2
  (Intercept)   0.4836213 -0.194518970 -0.093542902
  x1           -0.1945190  0.095774057  0.009148116
  x2           -0.0935429  0.009148116  0.149605995
  > vcov(m2b)
              (Intercept)         x1         x2
  (Intercept)    5.538131 -2.0225172 -2.7366732
  x1            -2.022517  0.7522853  0.9332862
  x2            -2.736673  0.9332862  1.7795414
  > 
  > str(predict(m1, se.fit = FALSE))
   Named num [1:10000] 0.658 0.615 1.73 0.908 -1.736 ...
   - attr(*, "names")= chr [1:10000] "4869" "4131" "8095" "2697" ...
  > str(predict(m1b, se.fit = FALSE))
   Named num [1:10000] 0.658 0.615 1.73 0.908 -1.736 ...
   - attr(*, "names")= chr [1:10000] "4869" "4131" "8095" "2697" ...
  > str(predict(m2, se.fit = FALSE))
   Named num [1:10000] 0.0252 0.9143 2.5706 1.683 1.2692 ...
   - attr(*, "names")= chr [1:10000] "1243" "4380" "3545" "3119" ...
  > str(predict(m2b, se.fit = FALSE))
   Named num [1:10000] 0.0252 0.9143 2.5706 1.683 1.2692 ...
   - attr(*, "names")= chr [1:10000] "1243" "4380" "3545" "3119" ...
  > 
  > str(predict(m1, se.fit = FALSE, newdata = dat1))
   Named num [1:11000] 0.6394 2.2747 -0.0414 1.2205 0.8107 ...
   - attr(*, "names")= chr [1:11000] "3817" "587" "1373" "6082" ...
  > str(predict(m1b, se.fit = FALSE, newdata = dat1))
   Named num [1:11000] 0.6394 2.2747 -0.0414 1.2205 0.8107 ...
   - attr(*, "names")= chr [1:11000] "3817" "587" "1373" "6082" ...
  > str(predict(m2, se.fit = FALSE, newdata = dat1))
   Named num [1:11000] 1.342 2.923 0.771 1.986 1.617 ...
   - attr(*, "names")= chr [1:11000] "3817" "587" "1373" "6082" ...
  > str(predict(m2b, se.fit = FALSE, newdata = dat1))
   Named num [1:11000] 1.342 2.923 0.771 1.986 1.617 ...
   - attr(*, "names")= chr [1:11000] "3817" "587" "1373" "6082" ...
  > 
  > ## cannot provide prediction Std Error without Bootstrap
  > tmp <- try(predict(m1, se.fit = TRUE))
  Error in predict.rsf(m1, se.fit = TRUE) : 
    cannot provide prediction Std Error without Bootstrap
  > stopifnot(inherits(tmp, "try-error"))
  > tmp <- try(predict(m2, se.fit = TRUE))
  Error in predict.rsf(m2, se.fit = TRUE) : 
    cannot provide prediction Std Error without Bootstrap
  > stopifnot(inherits(tmp, "try-error"))
  > tmp <- try(predict(m1, se.fit = TRUE, newdata = dat1))
  Error in predict.rsf(m1, se.fit = TRUE, newdata = dat1) : 
    cannot provide prediction Std Error without Bootstrap
  > stopifnot(inherits(tmp, "try-error"))
  > tmp <- try(predict(m2, se.fit = TRUE, newdata = dat1))
  Error in predict.rsf(m2, se.fit = TRUE, newdata = dat1) : 
    cannot provide prediction Std Error without Bootstrap
  > stopifnot(inherits(tmp, "try-error"))
  > 
  > str(predict(m1b, se.fit = TRUE))
  List of 2
   $ fit   : Named num [1:10000] 0.658 0.615 1.73 0.908 -1.736 ...
    ..- attr(*, "names")= chr [1:10000] "4869" "4131" "8095" "2697" ...
   $ se.fit: Named num [1:10000] 0.0743 0.0175 0.0394 0.0803 0.0601 ...
    ..- attr(*, "names")= chr [1:10000] "4869" "4131" "8095" "2697" ...
  > str(predict(m2b, se.fit = TRUE))
  List of 2
   $ fit   : Named num [1:10000] 0.0252 0.9143 2.5706 1.683 1.2692 ...
    ..- attr(*, "names")= chr [1:10000] "1243" "4380" "3545" "3119" ...
   $ se.fit: Named num [1:10000] 0.527 1.484 3.425 2.953 1.906 ...
    ..- attr(*, "names")= chr [1:10000] "1243" "4380" "3545" "3119" ...
  > str(predict(m1b, se.fit = TRUE, newdata = dat1))
  List of 2
   $ fit   : Named num [1:11000] 0.6394 2.2747 -0.0414 1.2205 0.8107 ...
    ..- attr(*, "names")= chr [1:11000] "3817" "587" "1373" "6082" ...
   $ se.fit: Named num [1:11000] 0.1251 0.0566 0.1159 0.0559 0.0608 ...
    ..- attr(*, "names")= chr [1:11000] "3817" "587" "1373" "6082" ...
  > str(predict(m2b, se.fit = TRUE, newdata = dat1))
  List of 2
   $ fit   : Named num [1:11000] 1.342 2.923 0.771 1.986 1.617 ...
    ..- attr(*, "names")= chr [1:11000] "3817" "587" "1373" "6082" ...
   $ se.fit: Named num [1:11000] 1.52 3.51 1.13 2.63 2.26 ...
    ..- attr(*, "names")= chr [1:11000] "3817" "587" "1373" "6082" ...
  > 
  > confint(m0)
       2.5% 97.5%
  > confint(m1)
          2.5%      97.5%
  x1 -1.067447 -0.9397884
  x2  0.232491  0.6686321
  > confint(m1b)
          2.5%      97.5%
  x1 -1.040089 -0.9630021
  x2  0.184628  0.5796997
  > confint(m2)
                    2.5%      97.5%
  (Intercept) -0.3318736  2.3941576
  x1          -1.5038760 -0.2907608
  x2          -0.6174995  0.8986868
  > confint(m2b)
                    2.5%      97.5%
  (Intercept) -0.3046132  6.7669927
  x1          -3.0660382 -0.4803545
  x2          -2.8109821  1.2574959
  > 
  > str(model.matrix(m1))
   num [1:11000, 1:2] -0.19 -2.023 0.422 -0.982 -0.569 ...
   - attr(*, "dimnames")=List of 2
    ..$ : chr [1:11000] "3817" "587" "1373" "6082" ...
    ..$ : chr [1:2] "x1" "x2"
  > str(model.matrix(m2))
   num [1:11000, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
   - attr(*, "dimnames")=List of 2
    ..$ : chr [1:11000] "7102" "629" "3683" "8460" ...
    ..$ : chr [1:3] "(Intercept)" "x1" "x2"
   - attr(*, "assign")= int [1:3] 0 1 2
  > 
  > str(model.frame(m1))
  'data.frame':	11000 obs. of  3 variables:
   $ status: num  1 1 1 1 1 1 1 1 1 1 ...
   $ x1    : num  -0.19 -2.023 0.422 -0.982 -0.569 ...
   $ x2    : num  0.995 0.542 0.849 0.521 0.531 ...
   - attr(*, "terms")=Classes 'terms', 'formula'  language status ~ (x1 + x2) - status
    .. ..- attr(*, "variables")= language list(status, x1, x2)
    .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
    .. .. ..- attr(*, "dimnames")=List of 2
    .. .. .. ..$ : chr [1:3] "status" "x1" "x2"
    .. .. .. ..$ : chr [1:2] "x1" "x2"
    .. ..- attr(*, "term.labels")= chr [1:2] "x1" "x2"
    .. ..- attr(*, "order")= int [1:2] 1 1
    .. ..- attr(*, "intercept")= int 1
    .. ..- attr(*, "response")= int 1
    .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
    .. ..- attr(*, "predvars")= language list(status, x1, x2)
    .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "numeric" "numeric"
    .. .. ..- attr(*, "names")= chr [1:3] "status" "x1" "x2"
  > str(model.frame(m2))
  'data.frame':	11000 obs. of  3 variables:
   $ status: num  1 1 1 1 1 1 1 1 1 1 ...
   $ x1    : num  -1.0106 0.9961 -0.6497 0.0718 -2.3008 ...
   $ x2    : num  0.4206 0.2424 0.0808 0.7848 0.2483 ...
   - attr(*, "terms")=Classes 'terms', 'formula'  language status ~ (x1 + x2) - status
    .. ..- attr(*, "variables")= language list(status, x1, x2)
    .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
    .. .. ..- attr(*, "dimnames")=List of 2
    .. .. .. ..$ : chr [1:3] "status" "x1" "x2"
    .. .. .. ..$ : chr [1:2] "x1" "x2"
    .. ..- attr(*, "term.labels")= chr [1:2] "x1" "x2"
    .. ..- attr(*, "order")= int [1:2] 1 1
    .. ..- attr(*, "intercept")= int 1
    .. ..- attr(*, "response")= int 1
    .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
    .. ..- attr(*, "predvars")= language list(status, x1, x2)
    .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "numeric" "numeric"
    .. .. ..- attr(*, "names")= chr [1:3] "status" "x1" "x2"
  > 
  > logLik(m0)
  'log Lik.' -9210.44 (df=0)
  > logLik(m1)
  'log Lik.' -8714.233 (df=2)
  > logLik(m2)
  'log Lik.' -9182.475 (df=3)
  > 
  > kdepairs(m1)
  Warning messages:
  1: In par(usr) : argument 1 does not name a graphical parameter
  2: In par(usr) : argument 1 does not name a graphical parameter
  3: In par(usr) : argument 1 does not name a graphical parameter
  > kdepairs(m2)
  Warning messages:
  1: In par(usr) : argument 1 does not name a graphical parameter
  2: In par(usr) : argument 1 does not name a graphical parameter
  3: In par(usr) : argument 1 does not name a graphical parameter
  > 
  > str(fitted(m1))
   Named num [1:11000] 1.895 9.725 0.959 3.389 2.249 ...
   - attr(*, "names")= chr [1:11000] "3817" "587" "1373" "6082" ...
  > str(fitted(m2))
   Named num [1:11000] 0.88 0.543 0.836 0.746 0.958 ...
   - attr(*, "names")= chr [1:11000] "7102" "629" "3683" "8460" ...
  > 
  > str(residuals(m1))
   Named num [1:11000] 0.916 0.679 0.956 0.859 0.902 ...
   - attr(*, "names")= chr [1:11000] "3817" "587" "1373" "6082" ...
  > str(residuals(m2))
   Named num [1:11000] 0.891 0.93 0.896 0.906 0.883 ...
   - attr(*, "names")= chr [1:11000] "7102" "629" "3683" "8460" ...
  > 
  > ## test no-selection model
  > f <- function(n=100, m=10) {
  +     y <- c(rep(1, n), rep(0, n*m))
  +     c(data=log(sum(y==1)/sum(y==0)),
  +         rsf=rsf.null(y, m=0)$results$par,
  +         glm=unname(coef(glm(y ~ 1, family=binomial))))
  + 
  + }
  > vals <- expand.grid(n=c(100, 200, 500), m=c(1,10,20))
  > tst <- t(apply(vals, 1, function(z) f(z[1], z[2])))
  > data.frame(vals, round(tst, 4))
      n  m    data     rsf     glm
  1 100  1  0.0000  0.0000  0.0000
  2 200  1  0.0000  0.0000  0.0000
  3 500  1  0.0000  0.0000  0.0000
  4 100 10 -2.3026 -2.3026 -2.3026
  5 200 10 -2.3026 -2.3026 -2.3026
  6 500 10 -2.3026 -2.3026 -2.3026
  7 100 20 -2.9957 -2.9957 -2.9957
  8 200 20 -2.9957 -2.9957 -2.9957
  9 500 20 -2.9957 -2.6962 -2.9957
  > 
  > proc.time()
     user  system elapsed 
    90.10    8.28   98.46