Hi
I had a look at the code.
First, there was a simple mistake in the model. The parameters were not multiplying any attributes, e.g. beta_age_wa should be multiplying age, etc. This was the cause for the identification issue.
Code: Select all
V[["walk_alone" ]] = asc_wa + beta_age_wa+ beta_female_wa + beta_distance_wa + beta_density_wa + beta_school_wa + beta_income_wa
Second, I optimised the code so that you can still produce the utilities (and also the betas) iteratively, and it all runs fine, including with analytical gradients. Please see below:
Code: Select all
rm(list = ls())
library(apollo)
library(readr)
library(here)
# Initialisation ===========================
#
#
### Initialise code
apollo_initialise()
### Set core controls
apollo_control = list(
modelName = "mnl_ind3_01",
modelDescr = "MNL on kid trip data - simplified independence with first covariates",
indivID = "id"
)
## class membership
# Data =====================================
# Read in data and do basic manipulations for coherence
#
data <- readr::read_rds(("usa-2017-all.rds"))
database <- data |>
dplyr::transmute(
id = 1:dplyr::n(),
# set up alternatives as named vector
choice = dplyr::case_when(
mode_ind_3 == 720 ~ "drive_parent",
mode_ind_3 == 730 ~ "drive_others",
mode_ind_3 == 810 ~ "walk_alone",
mode_ind_3 == 820 ~ "walk_parent",
mode_ind_3 == 830 ~ "walk_others",
mode_ind_3 == 910 ~ "bike_alone",
mode_ind_3 == 920 ~ "bike_parent",
mode_ind_3 == 930 ~ "bike_others",
TRUE ~ as.character(NA)
),
income = income_k, veh_per_driver, n_adults, non_work_mom, non_work_dad,
age, female, distance, school, density
) |>
dplyr::filter(!is.na(choice))
# Name settings
J = 8
alts = c("walk_alone","walk_others","walk_parent","bike_alone","bike_others","bike_parent",
"drive_others","drive_parent")
alt_short = c("wa","wo","wp","ba","bo","bp","do","dp")
attributes = c("age","female","distance","density","school","income")
# Parameters ===============================
# Create a list of all the parameters we need to estimate
#
apollo_beta=setNames(rep(0,J*(length(attributes)+1)),
c(paste0(rep("asc_",J),alt_short),paste0(rep(paste0("beta_",attributes),each=J),"_",rep(alt_short,length(attributes)))))
apollo_fixed = names(apollo_beta)[grepl("wa", names(apollo_beta))]
apollo_fixed=unique(apollo_fixed)
# Model Definition =========================
# verify all the parameters and data are correct and expected
apollo_inputs = apollo_validateInputs()
# include additional elements in apollo_inputs
apollo_inputs$J = 8
apollo_inputs$alts = alts
apollo_inputs$alt_short = alt_short
apollo_inputs$attributes = attributes
# Define model probabilities and likelihood
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))
### Define settings for MNL
mnl_settings = list(
alternatives = setNames(alts,alts),
choiceVar = choice
)
### Create list of probabilities P
P = list()
V = list()
### compute class-specific utilities
for(j in 1:apollo_inputs$J){
V[[apollo_inputs$alts[j]]] = get(paste0("asc_",apollo_inputs$alt_short[[j]]))
for(a in apollo_inputs$attributes) V[[apollo_inputs$alts[j]]]=V[[apollo_inputs$alts[j]]]+ get(paste0("beta_",a,"_",apollo_inputs$alt_short[[j]])) * get(a)
}
mnl_settings$utilities <- V
### Compute within-class choice probabilities using MNL model
P[["model"]] = apollo_mnl(mnl_settings, functionality)
### Prepare and return outputs of function
P = apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}
# ################################################################# #
#### MODEL ESTIMATION ####
# ################################################################# #
### Estimate model
model = apollo_estimate(apollo_beta, apollo_fixed,
apollo_probabilities, apollo_inputs)
### Show output in screen
apollo_modelOutput(model)
### Save output to file(s)
apollo_saveOutput(model)