Important: Read this before posting to this forum

  1. This forum is for questions related to the use of Apollo. We will answer some general choice modelling questions too, where appropriate, and time permitting. We cannot answer questions about how to estimate choice models with other software packages.
  2. There is a very detailed manual for Apollo available at http://www.ApolloChoiceModelling.com/manual.html. This contains detailed descriptions of the various Apollo functions, and numerous examples are available at http://www.ApolloChoiceModelling.com/examples.html. In addition, help files are available for all functions, using e.g. ?apollo_mnl
  3. Before asking a question on the forum, users are kindly requested to follow these steps:
    1. Check that the same issue has not already been addressed in the forum - there is a search tool.
    2. Ensure that the correct syntax has been used. For any function, detailed instructions are available directly in Apollo, e.g. by using ?apollo_mnl for apollo_mnl
    3. Check the frequently asked questions section on the Apollo website, which discusses some common issues/failures. Please see http://www.apollochoicemodelling.com/faq.html
    4. Make sure that R is using the latest official release of Apollo.
  4. If the above steps do not resolve the issue, then users should follow these steps when posting a question:
    1. provide full details on the issue, including the entire code and output, including any error messages
    2. posts will not immediately appear on the forum, but will be checked by a moderator first. This may take a day or two at busy times. There is no need to submit the post multiple times.

Delta method in random parameter models

Ask questions about post-estimation functions (e.g. prediction, conditionals, etc) or other processing of results.
DavidKL
Posts: 15
Joined: 27 Mar 2021, 16:30

Delta method in random parameter models

Post by DavidKL »

Hi David and Stephane!

How could I calculate standard errors for mean WTP-s in case of random parameter models?

I use the following formulation to calculate mean WTP-s:

Code: Select all

unconditionals=apollo_unconditionals(model,apollo_probabilities,apollo_inputs)
attribute=unconditionals[["b_attribute"]]
price=unconditionals[["b_price"]]
wtp=attribute/price
mean(wtp)
Can I use delta method for this in Apollo?

If so, how could I resolve this?

Thank you for your help!

David
dpalma
Posts: 190
Joined: 24 Apr 2020, 17:54

Re: Delta method in random parameter models

Post by dpalma »

Hi David,

No, you cannot use apollo_deltaMethod to calculate WTP when you have random coefficients. apollo_deltaMethod works only for fixed (as in non-random parameters). Instead, you have a couple of alternatives to calculate WTP when using random coefficients.
  • Do simulation: That is what you did in your code. For this to be valid, you need to make sure that the ratio between the coefficients has finite moments (i.e. finite mean and s.d.). This is true if your numerator (b_attribute) follows a normal distribution and your denominator (b_price) is a fixed parameters; or both are log-normals. But it is not defined if both numerator and denominator follow a normal distribution; and in several other cases. I recommend you looking at Daily, Hess and Train (2012) Assuring finite moments for willingness to pay in random coefficient models
Best
David
DavidKL
Posts: 15
Joined: 27 Mar 2021, 16:30

Re: Delta method in random parameter models

Post by DavidKL »

Hi David!

Thank you very much for your quick and helpful answer!

I have estimated a two-class LC model with random parameters (only within class).

I tried to use WTP-space approach, but class allocation probabilites were very different as in preference-space. So I assume that, my model in WTP-space converged for very different solution as in preference-space.

Consequently, I have tried to derive WTP-s by using the earlier presented formulation (for "attribute" I used normal, while for "price" used negative lognormal distribution) to calculate mean and standard deviation values. Then I have tried to calculate s.e.-s from the standard deviations.

Is that correct?

Best regards,
David
dpalma
Posts: 190
Joined: 24 Apr 2020, 17:54

Re: Delta method in random parameter models

Post by dpalma »

Hi David,

According to Daily et al. (2012) you should be ok doing simulation for a normal coefficient (b_attribute) divided by a negative lognormal (b_price). And yes, calculating the s.d. of the simulated values is an estimator of the s.d. of the WTP. Its precision will depend on the number of draws you are using.

However, you cannot calculate the s.e. of the WTP directly from the simulation. You would need to bootstrap your simulation. In other words, you would need to (i) estimate your model many times drawing subsamples (with repetition) from your full sample, (ii) calculate and record the WTP mean (WTP_mu) and s.d. (WTP_sd) for each subsample, and then (iii) approximate the s.e. by calculating the sd of WTP_mu and WTP_sd across subsamples. You could use apollo_bootstrap to estimate the model in different subsamples, but then you would have to manually calculate the WTP for each subsample. This could get a bit tricky, so I suggest trying to estimate the model in WTP-space again.

I am surprised that the preference and WTP space led to significantly different results. What starting values did you use for WTP-space? WTP-space is usually more difficult to estimate, so it is a good idea to first estimate the model in preference space, calculate the implied WTP, and then re-estimate the model in WTP-space using the implied WTP from preference space as starting values. That should make you start close to the same solution, so both solutions shouldn't diverge too much.

Cheers
David
DavidKL
Posts: 15
Joined: 27 Mar 2021, 16:30

Re: Delta method in random parameter models

Post by DavidKL »

Hi David!

I have tried to use several starting values, but neither LC nor RLC model did not converge to the same solution in the WTP space as it was estimated in the preference space (e.g. I got 0.15 - 0.85 class probability values in preference-space, while 0.47-0.53 in WTP-space).

I will try to apply bootstrap.

Thank you very much for your help!

Best wishes,
David
stephanehess
Site Admin
Posts: 974
Joined: 24 Apr 2020, 16:29

Re: Delta method in random parameter models

Post by stephanehess »

Hi

so in the Mixed Logit model, you combine Normal and Lognormal? Then the WTP space model will be different from preference space as the marginal utility for the non-price attribute will follow a distribution that is a product of Normal and Lognormal

Stephane
--------------------------------
Stephane Hess
www.stephanehess.me.uk
DavidKL
Posts: 15
Joined: 27 Mar 2021, 16:30

Re: Delta method in random parameter models

Post by DavidKL »

Hi Stephane,

thank you for the answer!

I get very different solution not only in latent class model with random parameters, but also in latent class with fix parameters.

My code and the results are the following in preference space:

Code: Select all

### Clear memory
rm(list = ls())

### Load Apollo library
library(apollo)

### Initialise code
apollo_initialise()

### Set core controls
apollo_control = list(
  modelName  ="LC",
  modelDescr ="LC_model_preference_space",
  indivID    ="ID",
  nCores     = 2,
  noDiagnostics = TRUE
)

# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS                     ####
# ################################################################# #

database = read.csv("data.csv",header=TRUE)

# ################################################################# #
#### DEFINE MODEL PARAMETERS                                     ####
# ################################################################# #

### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(asc_nc           = 0,
                b_eg_a       = 0,
                b_eg_b       = 0,
                b_sy_a       = 0,
                b_sy_b       = 0,
                b_gn_a       = 0,
                b_gn_b       = 0,
                b_an_a       = 0,
                b_an_b       = 0,
                b_price_a       = 0,
                b_price_b       = 0,
                delta_a         = 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()

# ################################################################# #
#### DEFINE LATENT CLASS COMPONENTS                              ####
# ################################################################# #

apollo_lcPars=function(apollo_beta, apollo_inputs){
  lcpars = list()
  
  lcpars[["b_eg"]] = list(b_eg_a, b_eg_b)
  lcpars[["b_sy"]] = list(b_sy_a, b_sy_b)
  lcpars[["b_gn"]] = list(b_gn_a, b_gn_b)
  lcpars[["b_an"]] = list(b_an_a, b_an_b)
  lcpars[["b_price"]] = list(b_price_a, b_price_b)
  
  V=list()
  V[["class_a"]] = delta_a
  V[["class_b"]] = 0
  
  mnl_settings = list(
    alternatives = c(class_a=1, class_b=2), 
    avail        = 1, 
    choiceVar    = NA, 
    V            = V
  )
  lcpars[["pi_values"]] = apollo_mnl(mnl_settings, functionality="raw")
  
  lcpars[["pi_values"]] = apollo_firstRow(lcpars[["pi_values"]], apollo_inputs)
  
  return(lcpars)
}

# ################################################################# #
#### 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 settings for MNL model component that are generic across classes
  mnl_settings = list(
    alternatives  = c(alt1=1, alt2=2, alt3=3, altnc=4), 
    avail         = list(alt1=av_alt1, alt2=av_alt2, alt3=av_alt3, altnc=av_altnc), 
    choiceVar     = choice
  )
  
  ### Loop over classes
  s=1
  while(s<=2){
    
    ### Compute class-specific utilities
    V=list()
    V[['alt1']]  = b_eg[[s]]  * ( l_alt1 == 2) + b_sy[[s]]  * ( l_alt1 == 3) + b_gn[[s]]  * ( h_alt1 == 2) + b_an[[s]]  * ( h_alt1 == 3) + b_price[[s]] * price_alt1                           
    V[['alt2']]  = b_eg[[s]]  * ( l_alt2 == 2) + b_sy[[s]]  * ( l_alt2 == 3) + b_gn[[s]]  * ( h_alt2 == 2) + b_an[[s]]  * ( h_alt2 == 3) + b_price[[s]] * price_alt2  
    V[['alt3']]  = b_eg[[s]]  * ( l_alt3 == 2) + b_sy[[s]]  * ( l_alt3 == 3) + b_gn[[s]]  * ( h_alt3 == 2) + b_an[[s]]  * ( h_alt3 == 3) + b_price[[s]] * price_alt3
    V[['altnc']] = asc_nc
    
    mnl_settings$V = V
    mnl_settings$componentName = paste0("Class_",s)
    
    ### Compute within-class choice probabilities using MNL model
    P[[paste0("Class_",s)]] = apollo_mnl(mnl_settings, functionality)
    
    ### Take product across observation for same individual
    P[[paste0("Class_",s)]] = apollo_panelProd(P[[paste0("Class_",s)]], apollo_inputs ,functionality)
    
    s=s+1
  }
  
  ### Compute latent class model probabilities
  lc_settings   = list(inClassProb = P, classProb=pi_values)
  P[["model"]] = apollo_lc(lc_settings, 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)

# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN)                               ----
# ----------------------------------------------------------------- #

modelOutput_settings = list()
modelOutput_settings$printPVal=TRUE
apollo_modelOutput(model,modelOutput_settings)

Code: Select all

Model run at                     : 2021-08-13 15:54:49
Estimation method                : bfgs
Model diagnosis                  : successful convergence 
Number of individuals            : 431
Number of rows in database       : 2586
Number of modelled outcomes      : 2586

Number of cores used             :  2 
Model without mixing

LL(start)                        : -3584.957
LL(0, whole model)               : -3584.957
LL(final, whole model)           : -2626.368
Rho-square (0)                   :  0.2674 
Adj.Rho-square (0)               :  0.264 
AIC                              :  5276.74 
BIC                              :  5347.03 

LL(0,Class_1)                    : -3584.957
LL(final,Class_1)                : -3690.434
LL(0,Class_2)                    : -3584.957
LL(final,Class_2)                : -3629.481

Estimated parameters             :  12
Time taken (hh:mm:ss)            :  00:00:26.97 
     pre-estimation              :  00:00:3.21 
     estimation                  :  00:00:14.33 
     post-estimation             :  00:00:9.43 
Iterations                       :  34  
Min abs eigenvalue of Hessian    :  9.871721 

Estimates:
             Estimate        s.e.   t.rat.(0)  p(1-sided)    Rob.s.e. Rob.t.rat.(0)  p(1-sided)
asc_nc       -6.48929     0.30488    -21.2844    0.000000      0.4934      -13.1521     0.00000
b_eg_a        1.00144     0.10087      9.9280    0.000000      0.1698        5.8985   1.834e-09
b_eg_b        0.43117     0.16505      2.6123    0.004497      0.2714        1.5889     0.05604
b_sy_a        1.27113     0.09032     14.0735    0.000000      0.1363        9.3282     0.00000
b_sy_b        0.88604     0.17588      5.0377   2.355e-07      0.4191        2.1141     0.01725
b_gn_a        0.11563     0.08469      1.3654    0.086071      0.1132        1.0215     0.15351
b_gn_b       -0.03723     0.16579     -0.2246    0.411153      0.2479       -0.1502     0.44031
b_an_a        0.19595     0.07816      2.5070    0.006089      0.1106        1.7717     0.03822
b_an_b        0.02492     0.15906      0.1567    0.437758      0.2212        0.1127     0.45515
b_price_a    -0.32459     0.08086     -4.0140   2.985e-05      0.2894       -1.1217     0.13100
b_price_b    -1.92705     0.10299    -18.7111    0.000000      0.1527      -12.6184     0.00000
delta_a       0.45783     0.22406      2.0433    0.020510      0.7618        0.6010     0.27391


Summary of class allocation for LC model component :
         Mean prob.
Class_1      0.6125
Class_2      0.3875
and in WTP-space:

Code: Select all

### Clear memory
rm(list = ls())

### Load Apollo library
library(apollo)

### Initialise code
apollo_initialise()

### Set core controls
apollo_control = list(
  modelName  ="LC",
  modelDescr ="LC_model_WTP_space",
  indivID    ="ID",
  nCores     = 2,
  noDiagnostics = TRUE
)

# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS                     ####
# ################################################################# #

database = read.csv("data.csv",header=TRUE)

# ################################################################# #
#### DEFINE MODEL PARAMETERS                                     ####
# ################################################################# #

### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(asc_nc           = 0,
                b_eg_a       = 0,
                b_eg_b       = 0,
                b_sy_a       = 0,
                b_sy_b       = 0,
                b_gn_a       = 0,
                b_gn_b       = 0,
                b_an_a       = 0,
                b_an_b       = 0,
                b_price_a       = 0,
                b_price_b       = 0,
                delta_a         = 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()

# ################################################################# #
#### DEFINE LATENT CLASS COMPONENTS                              ####
# ################################################################# #

apollo_lcPars=function(apollo_beta, apollo_inputs){
  lcpars = list()
  
  lcpars[["b_eg"]] = list(b_eg_a, b_eg_b)
  lcpars[["b_sy"]] = list(b_sy_a, b_sy_b)
  lcpars[["b_gn"]] = list(b_gn_a, b_gn_b)
  lcpars[["b_an"]] = list(b_an_a, b_an_b)
  lcpars[["b_price"]] = list(b_price_a, b_price_b)
  
  V=list()
  V[["class_a"]] = delta_a
  V[["class_b"]] = 0
  
  mnl_settings = list(
    alternatives = c(class_a=1, class_b=2), 
    avail        = 1, 
    choiceVar    = NA, 
    V            = V
  )
  lcpars[["pi_values"]] = apollo_mnl(mnl_settings, functionality="raw")
  
  lcpars[["pi_values"]] = apollo_firstRow(lcpars[["pi_values"]], apollo_inputs)
  
  return(lcpars)
}

# ################################################################# #
#### 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 settings for MNL model component that are generic across classes
  mnl_settings = list(
    alternatives  = c(alt1=1, alt2=2, alt3=3, altnc=4), 
    avail         = list(alt1=av_alt1, alt2=av_alt2, alt3=av_alt3, altnc=av_altnc), 
    choiceVar     = choice
  )
  
  ### Loop over classes
  s=1
  while(s<=2){
    
    ### Compute class-specific utilities
    V=list()
    V[['alt1']]  = b_price[[s]]*(b_eg[[s]]  * ( l_alt1 == 2) + b_sy[[s]]  * ( l_alt1 == 3) + b_gn[[s]]  * ( h_alt1 == 2) + b_an[[s]]  * ( h_alt1 == 3) + price_alt1)                           
    V[['alt2']]  = b_price[[s]]*(b_eg[[s]]  * ( l_alt2 == 2) + b_sy[[s]]  * ( l_alt2 == 3) + b_gn[[s]]  * ( h_alt2 == 2) + b_an[[s]]  * ( h_alt2 == 3) + price_alt2)  
    V[['alt3']]  = b_price[[s]]*(b_eg[[s]]  * ( l_alt3 == 2) + b_sy[[s]]  * ( l_alt3 == 3) + b_gn[[s]]  * ( h_alt3 == 2) + b_an[[s]]  * ( h_alt3 == 3) + price_alt3)
    V[['altnc']] = asc_nc
    
    mnl_settings$V = V
    mnl_settings$componentName = paste0("Class_",s)
    
    ### Compute within-class choice probabilities using MNL model
    P[[paste0("Class_",s)]] = apollo_mnl(mnl_settings, functionality)
    
    ### Take product across observation for same individual
    P[[paste0("Class_",s)]] = apollo_panelProd(P[[paste0("Class_",s)]], apollo_inputs ,functionality)
    
    s=s+1
  }
  
  ### Compute latent class model probabilities
  lc_settings   = list(inClassProb = P, classProb=pi_values)
  P[["model"]] = apollo_lc(lc_settings, 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)

# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN)                               ----
# ----------------------------------------------------------------- #

modelOutput_settings = list()
modelOutput_settings$printPVal=TRUE
apollo_modelOutput(model,modelOutput_settings)

Code: Select all

Model run at                     : 2021-08-13 15:51:26
Estimation method                : bfgs
Model diagnosis                  : successful convergence 
Number of individuals            : 431
Number of rows in database       : 2586
Number of modelled outcomes      : 2586

Number of cores used             :  2 
Model without mixing

LL(start)                        : -3584.957
LL(0, whole model)               : -3584.957
LL(final, whole model)           : -2595.486
Rho-square (0)                   :  0.276 
Adj.Rho-square (0)               :  0.2727 
AIC                              :  5214.97 
BIC                              :  5285.27 

LL(0,Class_1)                    : -3584.957
LL(final,Class_1)                : -5574.154
LL(0,Class_2)                    : -3584.957
LL(final,Class_2)                : -3042.518

Estimated parameters             :  12
Time taken (hh:mm:ss)            :  00:00:28.86 
     pre-estimation              :  00:00:3.63 
     estimation                  :  00:00:12.7 
     post-estimation             :  00:00:12.53 
Iterations                       :  35  
Min abs eigenvalue of Hessian    :  16.30801 

Estimates:
             Estimate        s.e.   t.rat.(0)  p(1-sided)    Rob.s.e. Rob.t.rat.(0)  p(1-sided)
asc_nc        -4.9942     0.22575     -22.123    0.000000     0.28488       -17.531    0.000000
b_eg_a        -0.4362     0.15954      -2.734    0.003129     0.18201        -2.397    0.008276
b_eg_b        -1.0927     0.10891     -10.033    0.000000     0.15612        -6.999   1.288e-12
b_sy_a        -0.4696     0.15576      -3.015    0.001284     0.15769        -2.978    0.001449
b_sy_b        -1.5501     0.10829     -14.315    0.000000     0.15548        -9.970    0.000000
b_gn_a         0.5337     0.18222       2.929    0.001701     0.20825         2.563    0.005192
b_gn_b        -0.1433     0.09515      -1.506    0.066082     0.11513        -1.244    0.106690
b_an_a         0.1910     0.14985       1.275    0.101241     0.15761         1.212    0.112811
b_an_b        -0.2065     0.08737      -2.364    0.009040     0.10750        -1.921    0.027350
b_price_a     -2.0322     0.11081     -18.339    0.000000     0.12278       -16.552    0.000000
b_price_b     -0.6759     0.03168     -21.335    0.000000     0.05058       -13.364    0.000000
delta_a       -2.1906     0.17188     -12.745    0.000000     0.18212       -12.028    0.000000


Summary of class allocation for LC model component :
         Mean prob.
Class_1      0.1006
Class_2      0.8994
I tried to use different starting values in WTP space model, but it didn't help.

David
stephanehess
Site Admin
Posts: 974
Joined: 24 Apr 2020, 16:29

Re: Delta method in random parameter models

Post by stephanehess »

David

LC is notoriously susceptible to end up in local optima. In this case, your WTP space model actually ends up in a better solution than the preference space model. But if you used the WTP space results in preference space, you would get the same LL

Did you also try estimation using the EM algorithm?

Stephane
--------------------------------
Stephane Hess
www.stephanehess.me.uk
DavidKL
Posts: 15
Joined: 27 Mar 2021, 16:30

Re: Delta method in random parameter models

Post by DavidKL »

Hi Stephane!

Thank you for the answer and help!

I will try to use the EM algorithm.

Best wishes,
David
aellyson
Posts: 5
Joined: 20 Oct 2021, 19:26

Re: Delta method in random parameter models

Post by aellyson »

dpalma wrote: 09 Aug 2021, 18:17 Hi David,

No, you cannot use apollo_deltaMethod to calculate WTP when you have random coefficients. apollo_deltaMethod works only for fixed (as in non-random parameters). Instead, you have a couple of alternatives to calculate WTP when using random coefficients.
  • Do simulation: That is what you did in your code. For this to be valid, you need to make sure that the ratio between the coefficients has finite moments (i.e. finite mean and s.d.). This is true if your numerator (b_attribute) follows a normal distribution and your denominator (b_price) is a fixed parameters; or both are log-normals. But it is not defined if both numerator and denominator follow a normal distribution; and in several other cases. I recommend you looking at Daily, Hess and Train (2012) Assuring finite moments for willingness to pay in random coefficient models
Best
David
I have a related but different question.

I have estimated mixed multinomial logistic regressions on a stated preference DCE with inter-individual draws for all attributes including time, the attribute we planned to use for willingness to spend time. Time was estimated using uniform draws, and the other attributes were estimated using normal draws. Dummy-coded models clearly show a non-linear specification for time use, so we have specified the time attribute as entering the utility function in quadratic form based on the following two publications:

van der Pol M, Ryan M. Specification of the Utility Function in Discrete Choice Experiments. Value Heal. 2014;17(2):297-301.
Torres C, Hanley N, Riera A. How wrong can you be? Implications of incorrect utility function specification for welfare measurement in choice experiments. J Environ Econ Manage. 2011;62:111-121.

I have tried conducting simulation for willingness to spend time as recommended in your post above, but it produces NaN for all simulated attribute values. I cannot tell if 1) I've made an error, 2) it is invalid to estimate WTP using an attribute whose estimated mean is centered around zero, 3) something else about quadratic specification needs to be revised, 4) something else I don't know. Code and output provided below.

I am also wondering if it is even possible to directly estimate a quadratic WTP attribute in WTP space in Apollo because of the plus/minus functionality in finding quadratic solutions. Do you have any examples of this or can you advise?

Code: Select all

#### Initialise code
> apollo_initialise()
Apollo ignition sequence completed
> 
> #### Check cores
> parallel::detectCores() 
[1] 4
> 
> #### Set core controls
> apollo_control = list(
+   modelName  ="MXL_3",
+   modelDescr ="Mixed logit (MXL) on DCE SP data, dummy-coded attributes, 
+                 reporting time continous and quadratic specification, 
+                 main effects only, all attributes + ASCs random, MLHS draws",
+   indivID    ="RID",
+   mixing    = TRUE, 
+   nCores    = 2
+   #analyticGrad = TRUE # Needs to be explicit for models with inter-intra draws (requires extra RAM).
+ )
> 
> # inter-individual draws are used for variation across individuals
> # intra-individual draws for variation across observations for the same individual
> 
> #### Set database
> database = df1
> 
> ###################################################################
> #### DEFINE MODEL PARAMETERS                                   ####
> ###################################################################
> 
> # Utility U_ij = sum[1 to K](beta_jk * X_jk) + epsilon_ij
> # 
> # where  U_ij is the ultity of profile j for individual i 
> #        X_jk where k = 1...K are attributes in the dce
> 
> ### Vector of parameters, including any that are kept fixed in estimation
> apollo_beta = c(mu_asc_a                    = -3,
+                 sigma_asc_a                 = -0.01,
+                 mu_asc_b                    = -3,
+                 sigma_asc_b                 = -0.01,
+                 b_violate_low               = 0,  # Violate_low is kept fixed
+                 mu_b_investigate            = -3,
+                 sigma_b_investigate_inter   = -0.01,
+                 mu_b_time                   = -3,
+                 sigma_b_time_inter          = -0.01,
+                 mu_b_time_sq                = -3,
+                 sigma_b_time_sq_inter       = -0.01,
+                 mu_b_anonymous              = -3,
+                 sigma_b_anonymous_inter     = -0.01,
+                 mu_b_amnesty                = -3,
+                 sigma_b_amnesty_inter       = -0.01,
+                 mu_b_violate_med            = -3,
+                 sigma_b_violate_med_inter   = -0.01,
+                 mu_b_violate_high           = -3,
+                 sigma_b_violate_high_inter  = -0.01)
> 
> ### 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("b_violate_low")
> 
> #####################################################################
> #### DEFINE RANDOM COMPONENTS                                    ####
> #####################################################################
> 
> ### Set parameters for generating draws
> apollo_draws = list(
+   interDrawsType = "MLHS",
+   interNDraws    = 500,
+   interUnifDraws = c("draws_time_inter", "draws_time_sq_inter"),
+   interNormDraws = c("draws_asc_a", "draws_asc_b",
+                      "draws_anonymous_inter", "draws_amnesty_inter",
+                      "draws_investigate_inter","draws_violate_med_inter",
+                      "draws_violate_high_inter"), 
+   intraDrawsType = "",
+   intraNDraws    = 0,
+   intraNormDraws = c(), 
+   intraUnifDraws = c()
+ )
> 
> ### Create random parameters
> apollo_randCoeff = function(apollo_beta, apollo_inputs){
+   randcoeff = list()
+   
+   randcoeff[["asc_a"]] = mu_asc_a +                   sigma_asc_a *                draws_asc_a
+   randcoeff[["asc_b"]] = mu_asc_b +                   sigma_asc_b *                draws_asc_b
+   randcoeff[["b_time"]] = mu_b_time +                 sigma_b_time_inter *         draws_time_inter
+   randcoeff[["b_time_sq"]] = mu_b_time_sq +           sigma_b_time_sq_inter *      draws_time_sq_inter
+   randcoeff[["b_anonymous"]] = mu_b_anonymous +       sigma_b_anonymous_inter *    draws_anonymous_inter
+   randcoeff[["b_amnesty"]] = mu_b_amnesty +           sigma_b_amnesty_inter *      draws_amnesty_inter
+   randcoeff[["b_investigate"]] = mu_b_investigate +   sigma_b_investigate_inter *  draws_investigate_inter
+   randcoeff[["b_violate_med"]] = mu_b_violate_med +   sigma_b_violate_med_inter *  draws_violate_med_inter
+   randcoeff[["b_violate_high"]] = mu_b_violate_high + sigma_b_violate_high_inter * draws_violate_high_inter
+   # randcoeff[[""]] = mu_b_ + sigma_b__inter * draws__inter
+   
+   # randcoeff[["b_investigate"]] = mu_b_investigate
+   #                                     + sigma_b_investigate_inter     * draws_investigate_inter
+   #                                     + sigma_b_investigate_intra     * draws_investigate_intra 
+   
+   return(randcoeff)
+ }
> 
> #####################################################################
> #### GROUP AND VALIDATE INPUTS                                   ####
> #####################################################################
> 
> apollo_inputs = apollo_validateInputs()
Several observations per individual detected based on the value of RID. Setting panelData in apollo_control
  set to TRUE.
All checks on apollo_control completed.
All checks on database completed.
Generating inter-individual draws ......... Done
> 
> #####################################################################
> #### DEFINE MODEL AND LIKELIHOOD FUNCTION                        ####
> #####################################################################
> 
> apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+   
+   ### Function initialisation: do not change the following three commands
+   ### 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()
+   
+   ### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
+   V = list() 
+   
+   V[['a']] = asc_a + b_anonymous * anonymous_a + b_time * time_a + b_time_sq * time_a * time_a + b_investigate * investigate_a + b_violate_low * (violate_a == 0) + b_violate_med * (violate_a == 1) + b_violate_high * (violate_a == 2) + b_amnesty * amnesty_a 
+   V[['b']] = asc_b + b_anonymous * anonymous_b + b_time * time_b + b_time_sq * time_b * time_b + b_investigate * investigate_b + b_violate_low * (violate_b == 0) + b_violate_med * (violate_b == 1) + b_violate_high * (violate_b == 2) + b_amnesty * amnesty_b 
+   V[['c']] = 0
+   
+   ### Define settings for MNL model component
+   mnl_settings = list(
+     alternatives  = c(a=1, b=2, c=3), 
+     avail         = list(a=1, b=1, c=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)
+   return(P)
+ }
> 
> #####################################################################
> #### MODEL ESTIMATION                                            ####
> #####################################################################
> 
> # ### Optional speedTest
> # speedTest_settings=list(
> #   nDrawsTry = c(100, 200, 500),
> #   nCoresTry = 1:10,
> #   nRep      = 10
> # )
> # 
> # apollo_speedTest(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs, speedTest_settings)
> 
> MXL_3 = apollo_estimate(apollo_beta, apollo_fixed,
+                         apollo_probabilities, apollo_inputs, 
+                         estimate_settings=list(hessianRoutine="maxLik"))

Testing likelihood function...

Overview of choices for MNL model component :
                                       a       b       c
Times available                  8032.00 8032.00 8032.00
Times chosen                     3412.00 3314.00 1306.00
Percentage chosen overall          42.48   41.26   16.26
Percentage chosen when available   42.48   41.26   16.26

Pre-processing likelihood function...
Preparing workers for multithreading...

Testing influence of parameters..................
Starting main estimation
Initial function value: -294565.4 
Initial gradient value:
                  mu_asc_a                sigma_asc_a                   mu_asc_b                sigma_asc_b 
                3411.98161                 -150.96857                 3313.98355                 -143.72177 
          mu_b_investigate  sigma_b_investigate_inter                  mu_b_time         sigma_b_time_inter 
                3876.99834                 -173.52127                17608.96516                 8514.50448 
              mu_b_time_sq      sigma_b_time_sq_inter             mu_b_anonymous    sigma_b_anonymous_inter 
               56644.96516                25272.75713                 3962.99842                 -178.90226 
              mu_b_amnesty      sigma_b_amnesty_inter           mu_b_violate_med  sigma_b_violate_med_inter 
                4031.99838                 -186.53884                 1717.00000                  -41.88887 
         mu_b_violate_high sigma_b_violate_high_inter 
                3502.99699                 -146.35327 
initial  value 294565.439042 
iter   2 value 45985.610138
iter   3 value 9908.723343
iter   4 value 9124.078980
iter   5 value 7984.555492
iter   6 value 7856.865565
iter   7 value 7736.338888
iter   8 value 7667.604144
iter   9 value 7650.549835
iter  10 value 7637.570218
iter  11 value 7566.365249
iter  12 value 7536.284985
iter  13 value 7458.709208
iter  14 value 7395.036999
iter  15 value 7271.401961
iter  16 value 7229.744813
iter  17 value 7201.647174
iter  18 value 7107.481110
iter  19 value 7077.874576
iter  20 value 6887.624995
iter  21 value 6877.596519
iter  22 value 6753.768691
iter  23 value 6724.593713
iter  24 value 6646.528796
iter  25 value 6541.358854
iter  26 value 6516.179882
iter  27 value 6501.455471
iter  28 value 6472.533695
iter  29 value 6466.063839
iter  30 value 6462.003358
iter  31 value 6458.907129
iter  32 value 6457.497947
iter  33 value 6456.925265
iter  34 value 6456.645028
iter  35 value 6456.498729
iter  36 value 6456.457763
iter  37 value 6456.398473
iter  38 value 6456.380287
iter  39 value 6456.379749
iter  40 value 6456.373976
iter  41 value 6456.370618
iter  42 value 6456.369985
iter  43 value 6456.369099
iter  44 value 6456.368444
iter  45 value 6456.367087
iter  46 value 6456.362775
iter  47 value 6456.361955
iter  48 value 6456.361509
iter  49 value 6456.361195
iter  50 value 6456.356222
iter  51 value 6456.353742
iter  52 value 6456.353581
iter  53 value 6456.353452
iter  53 value 6456.353389
iter  53 value 6456.353366
final  value 6456.353366 
converged
Additional convergence test using scaled estimation. Parameters will be scaled by their current estimates
  and additional iterations will be performed.
initial  value 6456.353366 
iter   2 value 6456.352860
iter   3 value 6456.352588
iter   3 value 6456.352566
iter   3 value 6456.352566
final  value 6456.352566 
converged
Estimated parameters:
                              Estimate
mu_asc_a                      -1.50323
sigma_asc_a                   -0.95509
mu_asc_b                      -1.59154
sigma_asc_b                   -0.97108
b_violate_low                  0.00000
mu_b_investigate               1.05741
sigma_b_investigate_inter     -1.16282
mu_b_time                     -0.17974
sigma_b_time_inter             1.45248
mu_b_time_sq                  -0.12585
sigma_b_time_sq_inter         -0.01072
mu_b_anonymous                 1.16539
sigma_b_anonymous_inter       -1.19584
mu_b_amnesty                   1.29390
sigma_b_amnesty_inter         -1.27754
mu_b_violate_med               1.44376
sigma_b_violate_med_inter     -0.93255
mu_b_violate_high              1.74279
sigma_b_violate_high_inter    -1.16227

Computing covariance matrix using numerical methods (maxLik). This may take a while, no progress bar
  displayed.
WARNING: Some eigenvalues of the Hessian are positive, indicating convergence to a saddle point!
Computing score matrix...
Calculating LL(0) for applicable models...
Calculating LL(c) for applicable models...
Calculating LL of each model component...
Warning message:
In sqrt(diag(varcov)) : NaNs produced
> 
> # ----------------------------------------------------------------- #
> #---- FORMATTED OUTPUT (TO SCREEN)                               ----
> # ----------------------------------------------------------------- #
> 
> apollo_modelOutput(MXL_3)
Model run using Apollo for R, version 0.2.6 on Darwin by aliceellyson 
www.ApolloChoiceModelling.com

Model name                       : MXL_3
Model description                : Mixed logit (MXL) on DCE SP data, dummy-coded attributes, 
                reporting time continous and quadratic specification, 
                main effects only, all attributes + ASCs random, MLHS draws
Model run at                     : 2022-01-12 15:29:50
Estimation method                : bfgs
Model diagnosis                  : successful convergence 
Number of individuals            : 1004
Number of rows in database       : 8032
Number of modelled outcomes      : 8032

Number of cores used             :  2 
Number of inter-individual draws : 500 (MLHS)

LL(start)                        : -294565.4
LL(0)                            : -8824.054
LL(C)                            : -8227.245
LL(final)                        : -6456.353
Rho-square (0)                   :  0.2683 
Adj.Rho-square (0)               :  0.2663 
Rho-square (C)                   :  0.2152 
Adj.Rho-square (C)               :  0.2131 
AIC                              :  12948.71 
BIC                              :  13074.55 

Estimated parameters             :  18
Time taken (hh:mm:ss)            :  00:31:21.28 
     pre-estimation              :  00:03:19.52 
     estimation                  :  00:07:7.04 
     post-estimation             :  00:20:54.71 
Iterations                       :  59  
Min abs eigenvalue of Hessian    :  8455.181 
Some eigenvalues of Hessian are positive, indicating potential problems!

Unconstrained optimisation.

Estimates:
                              Estimate        s.e.   t.rat.(0)    Rob.s.e. Rob.t.rat.(0)
mu_asc_a                      -1.50323     0.26855     -5.5975    0.297467       -5.0534
sigma_asc_a                   -0.95509     0.07988    -11.9573    0.097385       -9.8074
mu_asc_b                      -1.59154     0.27041     -5.8857    0.300556       -5.2953
sigma_asc_b                   -0.97108     0.08062    -12.0449    0.093964      -10.3346
b_violate_low                  0.00000          NA          NA          NA            NA
mu_b_investigate               1.05741     0.06812     15.5221    0.070229       15.0566
sigma_b_investigate_inter     -1.16282     0.08624    -13.4837    0.093665      -12.4147
mu_b_time                     -0.17974     0.26788     -0.6710    0.285514       -0.6295
sigma_b_time_inter             1.45248     0.08103     17.9251    0.084762       17.1361
mu_b_time_sq                  -0.12585     0.04955     -2.5395    0.054294       -2.3178
sigma_b_time_sq_inter         -0.01072         NaN         NaN    0.008842       -1.2122
mu_b_anonymous                 1.16539     0.06967     16.7264    0.075518       15.4320
sigma_b_anonymous_inter       -1.19584     0.08475    -14.1095    0.093386      -12.8053
mu_b_amnesty                   1.29390     0.07235     17.8831    0.078648       16.4518
sigma_b_amnesty_inter         -1.27754     0.08509    -15.0131    0.091237      -14.0024
mu_b_violate_med               1.44376     0.10844     13.3137    0.120448       11.9866
sigma_b_violate_med_inter     -0.93255     0.13866     -6.7253    0.161731       -5.7660
mu_b_violate_high              1.74279     0.08884     19.6165    0.096355       18.0872
sigma_b_violate_high_inter    -1.16227     0.08928    -13.0179    0.095641      -12.1524

> # Simulate values 
> N = 10000
> # Attributes
> b_attr_anon = rnorm(N, mean=MXL_3$estimate["mu_b_anonymous"],    sd=MXL_3$estimate["sigma_b_anonymous_inter"])
Warning message:
In rnorm(N, mean = MXL_3$estimate["mu_b_anonymous"], sd = MXL_3$estimate["sigma_b_anonymous_inter"]) :
  NAs produced
> b_attr_immu = rnorm(N, mean=MXL_3$estimate["mu_b_amnesty"],      sd=MXL_3$estimate["sigma_b_amnesty_inter"])
Warning message:
In rnorm(N, mean = MXL_3$estimate["mu_b_amnesty"], sd = MXL_3$estimate["sigma_b_amnesty_inter"]) :
  NAs produced
> b_attr_inve = rnorm(N, mean=MXL_3$estimate["mu_b_investigate"],  sd=MXL_3$estimate["sigma_b_investigate_inter"])
Warning message:
In rnorm(N, mean = MXL_3$estimate["mu_b_investigate"], sd = MXL_3$estimate["sigma_b_investigate_inter"]) :
  NAs produced
> b_attr_MED  = rnorm(N, mean=MXL_3$estimate["mu_b_violate_med"],  sd=MXL_3$estimate["sigma_b_violate_med_inter"])
Warning message:
In rnorm(N, mean = MXL_3$estimate["mu_b_violate_med"], sd = MXL_3$estimate["sigma_b_violate_med_inter"]) :
  NAs produced
> b_attr_HIGH = rnorm(N, mean=MXL_3$estimate["mu_b_violate_high"], sd=MXL_3$estimate["sigma_b_violate_high_inter"])
Warning message:
In rnorm(N, mean = MXL_3$estimate["mu_b_violate_high"], sd = MXL_3$estimate["sigma_b_violate_high_inter"]) :
  NAs produced
  
Post Reply