Page 1 of 1

Error in Apollo ModelOutput

Posted: 21 Dec 2022, 04:03
by Finley
Hello,

I am having an issue with displaying the output of the mixed logit model after running the model in a remote server, since my personal PC does not have enough memory for the data. When I run the model, the output is not displayed as I expect it to be. Instead, I am seeing an error message which says

Code: Select all

Error in round(model$rho2_0, 4) : 
  non-numeric argument to mathematical function
Strangely, it had worked out well before I updated the apollo package to 0.2.8 version or I directly used my computer to run a small fraction of data.

The model code is:

Code: Select all

# use apollo package to analyze
# ####################################################### #
#### 1. Definition of core settings
# ####################################################### #
### Clear memory
rm(list = ls())

### Install or Load library
if (!require('apollo')) install.packages('apollo')

library(apollo)

### Initialise code
apollo_initialise()

### Set core controls
apollo_control <- list(
  modelName = "Mixed logit model: Preference sapce",
  modelDescr = "Mixed logit model to determine Adults' preferences for a vaccine product",
  indivID = "ID",
  mixing = TRUE, # mixed logit model for random coefficients, FALSE for MNL model
  nCores = 6
  #,weights = "weight"
)
# ####################################################### #
#### 2. Data loading                                   ####
# ####################################################### #
database <- read.csv(file = "database.csv")

# small sample for testing 
# database <- database[1:1440, ]

database <- subset(database, select = c("ID","price.1","risk2.1","risk3.1","duration2.1",
                                        "duration3.1","efficacy2.1","efficacy3.1",
                                        "admin.1","doses2.1","doses3.1","origin.1",
                                        "price.2","risk2.2","risk3.2",
                                        "duration2.2","duration3.2","efficacy2.2",
                                        "efficacy3.2","admin.2","doses2.2",
                                        "doses3.2","origin.2","choice","weight"))

database <- as.data.frame(database)

# ####################################################### #
#### 3. Parameter definition                           ####
# ####################################################### #

### Vector of parameters, including any that are kept fixed
### during estimation

apollo_beta <- c(
  asc = 0,
  b_price_mu = 0,
  b_risk2_mu = 0,
  b_risk3_mu = 0,
  b_duration2_mu = 0,
  b_duration3_mu = 0,
  b_efficacy2_mu = 0,
  b_efficacy3_mu = 0,
  b_oral_mu = 0,
  b_dose2_mu = 0,
  b_dose3_mu = 0,
  b_imported_mu = 0,
  b_price_sigma = 0,
  b_risk2_sigma = 0,
  b_risk3_sigma = 0,
  b_duration2_sigma = 0,
  b_duration3_sigma = 0,
  b_efficacy2_sigma = 0,
  b_efficacy3_sigma = 0,
  b_oral_sigma = 0,
  b_dose2_sigma = 0,
  b_dose3_sigma = 0,
  b_imported_sigma = 0
)

### Vector with names (in quotes) of parameters to be kept fixed at their starting value in apollo_beta, use apollo_beta_fixed = c() if none
apollo_fixed <- c()

# ################################################################# #
#### 4. Define Random Components                                    ####
# ################################################################# #

### Set parameters for generating draws
apollo_draws <- list(
  interDrawsType = "mlhs",
  interNDraws = 500,
  interUnifDraws = c(),
  interNormDraws = paste0("draws_", c(
    "price", "risk2",
    "risk3", "duration2",
    "duration3", "efficacy2",
    "efficacy3", "oral",
    "dose2", "dose3",
    "imported"
  )),
  intraDrawsType = "",
  intraNDraws = 0,
  intraUnifDraws = c(),
  intraNormDraws = c()
)

### Create random parameters
apollo_randCoeff <- function(apollo_beta, apollo_inputs) {
  randcoeff <- list()
  
  randcoeff[["b_price"]] <- b_price_mu + b_price_sigma * draws_price
  randcoeff[["b_risk2"]] <- b_risk2_mu + b_risk2_sigma * draws_risk2
  randcoeff[["b_risk3"]] <- b_risk3_mu + b_risk3_sigma * draws_risk3
  randcoeff[["b_duration2"]] <- b_duration2_mu + b_duration2_sigma * draws_duration2
  randcoeff[["b_duration3"]] <- b_duration3_mu + b_duration3_sigma * draws_duration3
  randcoeff[["b_efficacy2"]] <- b_efficacy2_mu + b_efficacy2_sigma * draws_efficacy2
  randcoeff[["b_efficacy3"]] <- b_efficacy3_mu + b_efficacy3_sigma * draws_efficacy3
  randcoeff[["b_oral"]] <- b_oral_mu + b_oral_sigma * draws_oral
  randcoeff[["b_dose2"]] <- b_dose2_mu + b_dose2_sigma * draws_dose2
  randcoeff[["b_dose3"]] <- b_dose3_mu + b_dose3_sigma * draws_dose3
  randcoeff[["b_imported"]] <- b_imported_mu + b_imported_sigma * draws_imported
  
  return(randcoeff)
}

# ####################################################### #
#### 5. Input validation                               ####
# ####################################################### #

apollo_inputs <- apollo_validateInputs()

# ####################################################### #
#### 6. Likelihood definition                          ####
# ####################################################### #

apollo_probabilities <- function(apollo_beta, apollo_inputs, functionality = "estimate") {
  ### Attach inputs and detach after function exit
  apollo_attach(apollo_beta, apollo_inputs)
  on.exit(apollo_detach(apollo_beta, apollo_inputs))
  
  ### Create list of probabilities P
  P <- list()
  
  V <- list()
  
  ### List of utilities: these must use the same names as
  ### in mnl_settings, order is irrelevant.
  V[["alt1"]] <- asc + b_price * price.1 + b_risk2 * risk2.1 + b_risk3 * risk3.1 +
    b_duration2 * duration2.1 + b_duration3 * duration3.1 + b_efficacy2 * efficacy2.1 +
    b_efficacy3 * efficacy3.1 + b_oral * admin.1 + b_dose2 * doses2.1 +
    b_dose3 * doses3.1 + b_imported * origin.1
  V[["alt2"]] <- asc + b_price * price.2 + b_risk2 * risk2.2 + b_risk3 * risk3.2 +
    b_duration2 * duration2.2 + b_duration3 * duration3.2 + b_efficacy2 * efficacy2.2 +
    b_efficacy3 * efficacy3.2 + b_oral * admin.2 + b_dose2 * doses2.2 +
    b_dose3 * doses3.2 + b_imported * origin.2
  V[["alt3"]] <- 0
  
  ### Define settings for MNL model component
  mnl_settings <- list(
    alternatives = c(alt1 = 1, alt2 = 2, alt3 = 3),
    avail        = list(alt1 = 1, alt2 = 1, alt3 = 1),
    choiceVar    = choice,
    V            = V
  )
  
  ### Compute probabilities using MNL model
  P[["model"]] <- apollo_mnl(mnl_settings, functionality)
  
  ### Take product across observation for same individual
  P = apollo_panelProd(P, apollo_inputs, functionality)
  
  ### Average across inter-individual draws
  P = apollo_avgInterDraws(P, apollo_inputs, functionality)
  
  ### Prepare and return outputs of function
  P = apollo_prepareProb(P, apollo_inputs, functionality)
  
  ### Using sampling weights in estimation and prediction
  # P = apollo_weighting(P, apollo_inputs, functionality)
  
  return(P)
}

# ####################################################### #
#### 7. Model estimation                ####
# ####################################################### #
fit <- apollo_estimate(
  apollo_beta, apollo_fixed,
  apollo_probabilities, apollo_inputs
)

# save(fit, file = "ML_preference_space.RData") # If the file runs in the remote server
And the analysis code is:

Code: Select all

load(file = "Model/ML_preference_space.RData")

apollo_modelOutput(fit)
The full output is:

Image

And it seems to be nothing wrong with the model object "fit":

Image

R session info:

Code: Select all

> sessionInfo()
R version 4.2.0 (2022-04-22 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale:
[1] LC_COLLATE=Chinese (Simplified)_China.utf8  LC_CTYPE=Chinese (Simplified)_China.utf8   
[3] LC_MONETARY=Chinese (Simplified)_China.utf8 LC_NUMERIC=C                               
[5] LC_TIME=Chinese (Simplified)_China.utf8    

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] gpclib_1.5-6      aod_1.3.2         pscl_1.5.5        jtools_2.2.0      reshape2_1.4.4    rgdal_1.5-32     
 [7] mapproj_1.2.8     maptools_1.1-4    sp_1.4-7          mapdata_2.3.0     maps_3.4.0        apollo_0.2.8     
[13] car_3.0-12        carData_3.0-5     forcats_0.5.1     stringr_1.4.0     purrr_0.3.4       readr_2.1.2      
[19] tidyr_1.2.0       tibble_3.1.6      tidyverse_1.3.2   rio_0.5.29        janitor_2.1.0     mlogit_1.1-1     
[25] dfidx_0.0-4       haven_2.5.0       gtsummary_1.5.2   survey_4.1-1      weights_1.0.4     Hmisc_4.7-0      
[31] Formula_1.2-4     survival_3.3-1    fastR2_1.2.2      mosaic_1.8.3      ggridges_0.5.3    mosaicData_0.20.2
[37] ggformula_0.10.1  ggstance_0.3.5    dplyr_1.0.8       Matrix_1.4-1      ggplot2_3.3.5     lattice_0.20-45  
[43] readxl_1.4.0     

loaded via a namespace (and not attached):
  [1] backports_1.4.1     miscTools_0.6-26    plyr_1.8.7          splines_4.2.0       crosstalk_1.2.0    
  [6] leaflet_2.1.1       digest_0.6.29       htmltools_0.5.2     gdata_2.18.0.1      fansi_1.0.3        
 [11] magrittr_2.0.3      checkmate_2.0.0     googlesheets4_1.0.0 cluster_2.1.3       openxlsx_4.2.5     
 [16] tzdb_0.3.0          mosaicCore_0.9.0    modelr_0.1.8        matrixStats_0.62.0  MCMCpack_1.6-3     
 [21] sandwich_3.0-1      jpeg_0.1-9          colorspace_2.0-3    rvest_1.0.2         ggrepel_0.9.1      
 [26] mitools_2.4         rbibutils_2.2.8     xfun_0.30           jsonlite_1.8.0      crayon_1.5.1       
 [31] lme4_1.1-29         zoo_1.8-10          glue_1.6.2          polyclip_1.10-0     gargle_1.2.0       
 [36] gtable_0.3.0        MatrixModels_0.5-0  maxLik_1.5-2        abind_1.4-5         SparseM_1.81       
 [41] scales_1.2.0        mvtnorm_1.1-3       DBI_1.1.2           Rcpp_1.0.8.3        htmlTable_2.4.0    
 [46] tmvnsim_1.0-2       foreign_0.8-82      randtoolbox_1.31.1  httr_1.4.2          htmlwidgets_1.5.4  
 [51] RColorBrewer_1.1-3  ellipsis_0.3.2      mice_3.14.0         pkgconfig_2.0.3     farver_2.1.0       
 [56] dbplyr_2.1.1        nnet_7.3-17         utf8_1.2.2          tidyselect_1.2.0    rlang_1.0.6        
 [61] munsell_0.5.0       cellranger_1.1.0    tools_4.2.0         cli_3.3.0           generics_0.1.2     
 [66] pacman_0.5.1        broom_0.8.0         evaluate_0.15       fastmap_1.1.0       ggdendro_0.1.23    
 [71] yaml_2.3.5          mcmc_0.9-7          fs_1.5.2            knitr_1.38          pander_0.6.5       
 [76] zip_2.2.0           nlme_3.1-157        quantreg_5.88       rngWELL_0.10-7      xml2_1.3.3         
 [81] compiler_4.2.0      rstudioapi_0.13     curl_4.3.2          png_0.1-7           gt_0.5.0           
 [86] reprex_2.0.1        statmod_1.4.36      broom.helpers_1.7.0 tweenr_1.0.2        stringi_1.7.6      
 [91] nloptr_2.0.0        vctrs_0.4.1         pillar_1.7.0        lifecycle_1.0.3     Rdpack_2.3         
 [96] lmtest_0.9-40       data.table_1.14.2   R6_2.5.1            latticeExtra_0.6-29 gridExtra_2.3      
[101] RSGHB_1.2.2         boot_1.3-28         MASS_7.3-56         gtools_3.9.2        assertthat_0.2.1   
[106] withr_2.5.0         mnormt_2.0.2        parallel_4.2.0      hms_1.1.1           rpart_4.1.16       
[111] labelled_2.9.0      coda_0.19-4         minqa_1.2.4         rmarkdown_2.13      snakecase_0.11.0   
[116] googledrive_2.0.0   ggforce_0.3.3       numDeriv_2016.8-1.1 lubridate_1.8.0     base64enc_0.1-3  
I am wondering if anyone has experienced similar issues and knows how to fix this problem. I would greatly appreciate any advice or suggestions you may have.

Thank you in advance for your help!

Re: Error in Apollo ModelOutput

Posted: 28 Feb 2023, 22:44
by dpalma
Hi,

Sorry for taking so long to reply.

I suspect there was a problem calculating the likelihood of the equi-probable model (LL0), or the constants-only model (LLC). We will look into avoiding these issues in v0.2.9.

For now, I would recommend that before you call apollo_modelOutput, you do as follows:

Code: Select all

model$LL0 <- 1
model$LLC <- 1
model$rho2_0 <- 1
model$rho2_C <- 1
apollo_modelOutput(model)
The LLC, LL0, rho2 and adjusted rho2 reported will be wrong, but at least you will be able to see the rest of the output.

Cheers
David