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
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
Code: Select all
load(file = "Model/ML_preference_space.RData")
apollo_modelOutput(fit)
And it seems to be nothing wrong with the model object "fit":
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
Thank you in advance for your help!