Error in MDCEV with outside good
Posted: 03 Mar 2024, 16:59
Hello.
I am running an MDCEV model with outside good for traveling time use. I have issues with list of the alternatives and the gamma.
Here is my code.
Here is the results.
I carefully follow the example code of MDCEV_with_outside_good on your website. I set the alpha_base to (-40.27) and put it in the apollo_fixed based on your suggestion here viewtopic.php?t=808.
But then I got an error as above, list of alternatives != list of gamma. If I make the alternative and the gamma the same, doesn't that become an MDCEV model without outside good?
Would you help me with this issue?
Also, would you let me know why did you suggest to set the alpha_base to (-40.27) and put it in the apollo_fixed based as you mentioned here viewtopic.php?t=808?
I ran the example code of MDCEV_with_outside_good on your website, and got an error exactly like here viewtopic.php?t=808. After I set the alpha_base to (-40.27) and put it in the apollo_fixed. I got the warning as this.
Would you explain what is going on with the example code of MDCEV_with_outside_good?
I'd appreciate your help.
Best,
Naomi
I am running an MDCEV model with outside good for traveling time use. I have issues with list of the alternatives and the gamma.
Here is my code.
Code: Select all
# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS ####
# ################################################################# #
### Clear memory
# rm(list = ls())
### Load Apollo library
library(apollo)
### Initialise code
apollo_initialise()
### Set core controls
apollo_control = list(
modelName = "MDCEV_care",
modelDescr = "MDCEV model on time use data, alpha-gamma profile with outside good and socio-demographics",
indivID = "TUCASEID",
outputDirectory = "output"
)
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #
### Loading data from package
### if data is to be loaded from a file (e.g. called data.csv),
### the code would be: database = read.csv("data.csv",header=TRUE)
database = read.csv("atus2022_mdcev.csv",header=TRUE)
### Create consumption variables for combined activities
# database$t_outside = rowSums(database[,c("t_a01", "t_a06", "t_a10", "t_a11", "t_a12")]) # outside good: time spent at home and travelling
# database$t_leisure = rowSums(database[,c("t_a07", "t_a08", "t_a09")])
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(alpha_base = -40.27,
gamma_trv_personal = 1,
gamma_trv_care = 1,
gamma_trv_work = 1,
gamma_trv_school = 1,
gamma_trv_leisure = 1,
gamma_trv_other = 1,
delta_trv_personal = 0,
delta_trv_care = 0,
delta_trv_work = 0,
delta_trv_school = 0,
delta_trv_leisure = 0,
delta_trv_other = 1,
delta_age = 0,
delta_gender = 0,
sig = 1)
### 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("alpha_base", "sig")
# ################################################################# #
#### GROUP AND VALIDATE INPUTS ####
# ################################################################# #
apollo_inputs = apollo_validateInputs()
# ################################################################# #
#### DEFINE MODEL AND LIKELIHOOD FUNCTION ####
# ################################################################# #
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()
### Define individual alternatives
alternatives = c("activities",
"trv_personal",
"trv_care",
"trv_work",
"trv_school",
"trv_leisure",
"trv_other")
### Define availabilities
avail = list(activities = av_act,
trv_personal = av_trv_personal,
trv_care = av_trv_care,
trv_work = av_trv_work,
trv_school = av_trv_school,
trv_leisure = av_trv_leisure,
trv_other = av_trv_other)
### Define continuous consumption for individual alternatives
continuousChoice = list(activities = act/60,
trv_personal = trv_personal/60,
trv_care = trv_care/60,
trv_work = trv_work/60,
trv_school = trv_school/60,
trv_leisure = trv_leisure/60,
trv_other = trv_other/60)
### Define utilities for individual alternatives
V = list()
V[["activities" ]] = 0
V[["trv_personal" ]] = delta_trv_personal + delta_age * TEAGE + delta_gender * TESEX
V[["trv_care" ]] = delta_trv_care + delta_age * TEAGE + delta_gender * TESEX
V[["trv_work" ]] = delta_trv_work + delta_age * TEAGE + delta_gender * TESEX
V[["trv_school" ]] = delta_trv_school + delta_age * TEAGE + delta_gender * TESEX
V[["trv_leisure" ]] = delta_trv_leisure + delta_age * TEAGE + delta_gender * TESEX
V[["trv_other" ]] = delta_trv_other + delta_age * TEAGE + delta_gender * TESEX
# V[["outside"]] = 0
# V[["work"]] = delta_work + delta_work_FT * occ_full_time + delta_work_wknd * weekend
# V[["school"]] = delta_school + delta_school_young * (age<=30)
# V[["shopping"]] = delta_shopping
# V[["private"]] = delta_private
# V[["leisure"]] = delta_leisure + delta_leisure_wknd*weekend
### Define alpha parameters
alpha = list(activities = 1 /(1 + exp(-alpha_base)),
trv_personal = 1 /(1 + exp(-alpha_base)),
trv_care = 1 /(1 + exp(-alpha_base)),
trv_work = 1 /(1 + exp(-alpha_base)),
trv_school = 1 /(1 + exp(-alpha_base)),
trv_leisure = 1 /(1 + exp(-alpha_base)),
trv_other = 1 /(1 + exp(-alpha_base)))
### Define gamma parameters
gamma = list(trv_personal = gamma_trv_personal,
trv_care = gamma_trv_care,
trv_work = gamma_trv_work,
trv_school = gamma_trv_school,
trv_leisure = gamma_trv_leisure,
trv_other = gamma_trv_other)
### Define costs for individual alternatives
cost = list(activities = 1,
trv_personal = 1,
trv_care = 1,
trv_work = 1,
trv_school = 1,
trv_leisure = 1,
trv_other = 1)
### Define settings for MDCEV model
mdcev_settings <- list(alternatives = alternatives,
avail = avail,
continuousChoice = continuousChoice,
utilities = V,
alpha = alpha,
gamma = gamma,
sigma = sig,
cost = cost,
budget = 24)
### Compute probabilities using MDCEV model
P[["model"]] = apollo_mdcev(mdcev_settings, functionality)
### Take product across observation for same individual
# P = apollo_panelProd(P, apollo_inputs, functionality)
### Prepare and return outputs of function
P = apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}
# ################################################################# #
#### MODEL ESTIMATION ####
# ################################################################# #
model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO FILE, using model name) ----
# ----------------------------------------------------------------- #
# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS ####
# ################################################################# #
### Clear memory
# rm(list = ls())
### Load Apollo library
library(apollo)
### Initialise code
apollo_initialise()
### Set core controls
apollo_control = list(
modelName ="MDCEV_care",
modelDescr ="MDCEV model on time use data, alpha-gamma profile with outside good and socio-demographics",
indivID ="TUCASEID",
outputDirectory = "output"
)
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #
### Loading data from package
### if data is to be loaded from a file (e.g. called data.csv),
### the code would be: database = read.csv("data.csv",header=TRUE)
database = read.csv("atus2022_mdcev.csv",header=TRUE)
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(alpha_base = -40.27,
gamma_trv_personal = 1,
gamma_trv_care = 1,
gamma_trv_work = 1,
gamma_trv_school = 1,
gamma_trv_leisure = 1,
gamma_trv_other = 1,
delta_personal = 0,
delta_care = 0,
delta_work = 0,
delta_school = 0,
delta_leisure = 0,
delta_other = 0,
delta_age = 0,
delta_gender = 0,
sig = 1)
### 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("alpha_base","sig")
# ################################################################# #
#### GROUP AND VALIDATE INPUTS ####
# ################################################################# #
apollo_inputs = apollo_validateInputs()
# ################################################################# #
#### DEFINE MODEL AND LIKELIHOOD FUNCTION ####
# ################################################################# #
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()
### Define individual alternatives
alternatives = c("activities",
"trv_personal",
"trv_care",
"trv_work",
"trv_school",
"trv_leisure",
"trv_other")
### Define availabilities
avail = list(activities = av_act,
trv_personal = av_trv_personal,
trv_care = av_trv_care,
trv_work = av_trv_work,
trv_school = av_trv_school,
trv_leisure = av_trv_leisure,
trv_other = av_trv_other)
### Define continuous consumption for individual alternatives
continuousChoice = list(activities = act/60,
trv_personal = trv_personal/60,
trv_care = trv_care/60,
trv_work = trv_work/60,
trv_school = trv_school/60,
trv_leisure = trv_leisure/60,
trv_other = trv_other/60)
### Define utilities for individual alternatives
V = list()
V[["activities"]] = 0
V[["trv_personal" ]] = delta_personal + delta_age * TEAGE + delta_gender * TESEX
V[["trv_care" ]] = delta_care + delta_age * TEAGE + delta_gender * TESEX
V[["trv_work" ]] = delta_work + delta_age * TEAGE + delta_gender * TESEX
V[["trv_school" ]] = delta_school + delta_age * TEAGE + delta_gender * TESEX
V[["trv_leisure" ]] = delta_leisure + delta_age * TEAGE + delta_gender * TESEX
V[["trv_other" ]] = delta_other + delta_age * TEAGE + delta_gender * TESEX
# V[["work"]] = delta_work + delta_work_FT * occ_full_time + delta_work_wknd * weekend
# V[["school"]] = delta_school + delta_school_young * (age<=30)
# V[["shopping"]] = delta_shopping
# V[["private"]] = delta_private
# V[["leisure"]] = delta_leisure + delta_leisure_wknd*weekend
### Define alpha parameters
alpha = list(activities = 1 /(1 + exp(-alpha_base)),
trv_personal = 1 /(1 + exp(-alpha_base)),
trv_care = 1 /(1 + exp(-alpha_base)),
trv_work = 1 /(1 + exp(-alpha_base)),
trv_school = 1 /(1 + exp(-alpha_base)),
trv_leisure = 1 /(1 + exp(-alpha_base)),
trv_other = 1 /(1 + exp(-alpha_base)))
### Define gamma parameters
gamma = list(
trv_personal = gamma_trv_personal,
trv_care = gamma_trv_care,
trv_work = gamma_trv_work,
trv_school = gamma_trv_school,
trv_leisure = gamma_trv_leisure,
trv_other = gamma_trv_other)
### Define costs for individual alternatives
cost = list(activities = 1,
trv_personal = 1,
trv_care = 1,
trv_work = 1,
trv_school = 1,
trv_leisure = 1,
trv_other = 1)
### Define settings for MDCEV model
mdcev_settings <- list(alternatives = alternatives,
avail = avail,
continuousChoice = continuousChoice,
utilities = V,
alpha = alpha,
gamma = gamma,
sigma = sig,
cost = cost,
budget = 24)
### Compute probabilities using MDCEV model
P[["model"]] = apollo_mdcev(mdcev_settings, functionality)
### Take product across observation for same individual
# P = apollo_panelProd(P, apollo_inputs, functionality)
### Prepare and return outputs of function
P = apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}
# ################################################################# #
#### MODEL ESTIMATION ####
# ################################################################# #
model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO FILE, using model name) ----
# ----------------------------------------------------------------- #
# apollo_saveOutput(model)
Here is the results.
Code: Select all
> # ################################################################# #
> #### LOAD LIBRARY AND DEFINE CORE SETTINGS ####
> # ## .... [TRUNCATED]
> ### Initialise code
> apollo_initialise()
Apollo ignition sequence completed
> ### Set core controls
> apollo_control = list(
+ modelName = "MDCEV_care",
+ modelDescr = "MDCEV model on time use data, alpha-gamma profil ..." ... [TRUNCATED]
> # ################################################################# #
> #### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
> # ## .... [TRUNCATED]
> ### Create consumption variables for combined activities
> # database$t_outside = rowSums(database[,c("t_a01", "t_a06", "t_a10", "t_a11", "t_a12")]) .... [TRUNCATED]
> ### 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
> apoll .... [TRUNCATED]
> # ################################################################# #
> #### GROUP AND VALIDATE INPUTS ####
> # ## .... [TRUNCATED]
All checks on apollo_control completed.
All checks on database completed.
> # ################################################################# #
> #### DEFINE MODEL AND LIKELIHOOD FUNCTION ####
> # ## .... [TRUNCATED]
> # ################################################################# #
> #### MODEL ESTIMATION ####
> # ## .... [TRUNCATED]
WARNING: Element alpha_base in 'apollo_fixed' is constrained to a value other than zero or one. This
may be intentional. If not, stop this function by pressing the "Escape" key and adjust the
starting values accordingly.
Current process will resume in 5 seconds unless interrupted by the user.....
Preparing user-defined functions.
Error in if (!test) { : missing value where TRUE/FALSE needed
In addition: Warning messages:
1: In mdcev_settings$alternatives != names(mdcev_settings$gamma) :
longer object length is not a multiple of shorter object length
2: In mdcev_settings$alternatives != names(mdcev_settings$gamma) :
longer object length is not a multiple of shorter object lengthBut then I got an error as above, list of alternatives != list of gamma. If I make the alternative and the gamma the same, doesn't that become an MDCEV model without outside good?
Would you help me with this issue?
Also, would you let me know why did you suggest to set the alpha_base to (-40.27) and put it in the apollo_fixed based as you mentioned here viewtopic.php?t=808?
I ran the example code of MDCEV_with_outside_good on your website, and got an error exactly like here viewtopic.php?t=808. After I set the alpha_base to (-40.27) and put it in the apollo_fixed. I got the warning as this.
Code: Select all
Warning messages:
1: In log((inputs$continuousChoice[[j]]/inputs$gamma[[j]]) + 1) :
NaNs produced
2: In log(inputs$continuousChoice[[j]] + inputs$gamma[[j]]) :
NaNs produced
3: In log((inputs$continuousChoice[[j]]/inputs$gamma[[j]]) + 1) :
NaNs produced
4: In log(inputs$continuousChoice[[j]] + inputs$gamma[[j]]) :
NaNs producedI'd appreciate your help.
Best,
Naomi