R Under development (unstable) (2024-10-16 r87241 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 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. > library(mlogitBMA) Loading required package: BMA Loading required package: survival Loading required package: leaps Loading required package: robustbase Attaching package: 'robustbase' The following object is masked from 'package:survival': heart Loading required package: inline Loading required package: rrcov Scalable Robust Estimators with High Breakdown Point (version 1.7-6) Loading required package: abind Loading required package: maxLik Loading required package: miscTools Attaching package: 'miscTools' The following objects are masked from 'package:robustbase': colMedians, rowMedians Please cite the 'maxLik' package as: Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1. If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site: https://r-forge.r-project.org/projects/maxlik/ > > test.specification <- function() { + spec <- mnl.spec(choice ~ fuel + price + cost | hsg2, Car, + varying=5:ncol(Car), sep='') + stopifnot(all(dim(spec$variable.used) == c(6,4))) + stopifnot(all(spec$varying.names == c('fuel', 'price', 'cost'))) + stopifnot(sum(spec$same.coefs) == 3 & !spec$same.coefs['hsg2']) + stopifnot(sum(spec$intercepts)==5 & !spec$intercepts[1]) + cat('\nSpecification test OK.\n') + } > > test.estimation.car <- function() { + cat('\nRunning test for MNL estimation on Car dataset ... ') + est <- estimate.mlogit(choice ~ price + cost | coml5, Car, + varying=5:ncol(Car), sep='') + sest <- summary(est) + stopifnot(all(dim(sest$coefs) == c(12, 4))) + stopifnot(sest$lratio > 0.1) + cat(' OK.\n') + } > > test.estimation.catsup <- function() { + cat('\nRunning test for MNL estimation on Catsup dataset ... ') + est <- estimate.mlogit(choice ~ disp + feat + price, Catsup, + varying=2:(ncol(Catsup)-1), sep='.') + sest <- summary(est) + stopifnot(all(dim(sest$coefs) == c(6, 4))) + stopifnot(sest$lratio > 0.3) + stopifnot(sest$bic > 5083) + print(sest) + cat(' OK.\n') + } > > test.bic.mlogit.car <- function() { + cat('\nRunning test for BMA on Car dataset ... ') + res <- bic.mlogit(choice ~ price + cost + speed + acc + size | hsg2, Car, + varying=5:ncol(Car), sep='', include.intercepts = FALSE, + verbose = TRUE) + stopifnot(all(dim(res$bic.glm$which)==c(1,14))) # 1 model selected, 14 variables in total + cat('... BMA test OK.\n') + } > > test.bic.mlogit.catsup <- function() { + cat('\nRunning test for BMA on Catsup dataset ... ') + res <- bic.mlogit(choice ~ 1 | disp + feat + price, Catsup, + varying=2:(ncol(Catsup)-1), sep='.', + base.choice = 4, + include.intercepts = FALSE, + verbose = TRUE) + summary(res) + stopifnot(all(dim(res$bic.glm$which)==c(2,11))) # 2 models selected, 11 variables in total + cat('... BMA test OK.\n') + est.res <- estimate.mlogit(res, Catsup) + stopifnot(length(est.res)==2) + stopifnot(all(dim(est.res[[1]]$coefs) == c(12, 4))) + stopifnot(all(dim(est.res[[2]]$coefs) == c(11, 4))) + cat('Estimation on the BMA object OK.\n') + } > > # load data > data('Car', package='mlogit') > # convert the choice column into a numerical code, > # since it is the way the alternative-specific variables are constructed > Car[,'choice'] <- as.integer(gsub('^choice', '', Car[,'choice'])) > data('Catsup', package='mlogit') > > test.specification() Specification test OK. > test.estimation.car() Running test for MNL estimation on Car dataset ... OK. > test.estimation.catsup() Running test for MNL estimation on Catsup dataset ... successive function values within relative tolerance limit (reltol) Log-Likelihood: -2518 Null Log-Likelihood: -3879 Likelihood ratio index: 0.3509 AIC: 5047.755 BIC: 5083.374 Sample size: 2798 Iterations: 7 Suggested |t-value| > 2.817208 Convergence statistics: 1.619402e-06 Estimated using BHHH in 0.7920599 s. Coefficients : Estimate Std. Error t-value Pr(>|t|) asc.heinz41 -1.0723 0.08669 -12.368 0 asc.heinz32 -0.9248 0.07622 -12.133 0 asc.hunts32 -2.4261 0.10312 -23.527 0 disp 0.8756 0.09249 9.466 0 feat 0.9086 0.10867 8.361 0 price -1.4024 0.06044 -23.205 0 Response variable: choice Base choice name: heinz28 Base choice index: 1 Frequency of alternatives: heinz28 heinz41 heinz32 hunts32 851 182 1458 307 Equations: ---------- alternative intercept 1 2 3 1 heinz28: 2 heinz41: asc.heinz41 + disp + feat + price 3 heinz32: asc.heinz32 + disp + feat + price 4 hunts32: asc.hunts32 + disp + feat + price OK. > test.bic.mlogit.car() Running test for BMA on Car dataset ... Begg & Gray approximation started. 14 variables considered. 19 models initially selected. Final number of models: 1 ... BMA test OK. > # turn off to speed the tests up > #test.bic.mlogit.catsup() > > proc.time() user system elapsed 10.01 0.95 10.98