ApolloChoiceModelling forum

Forum for users of the Apollo choice modelling software
http://www.apollochoicemodelling.com/forum/

Unreasonable estimate for cost

http://www.apollochoicemodelling.com/forum/viewtopic.php?t=699

Page 1 of 2

Unreasonable estimate for cost

Posted: 09 Jun 2023, 07:15
by Nishihara
Dear Sir,

I am running a mixed logit model and getting a very unreasonable estimate for the log-normalized cost: -1066.11898

Please see the code and results. I cannot figure out what is the issue here because I already run the similar codes with different variable names in different datasets (same project).

It would be highly appreciated if you could let me know the reason behind it and how to fix it.

Thank you very much in advance.

Hara

Code:

Code: Select all

library(readr)
library(support.CEs)
library(survival)
library(openxlsx)
library(tidyr)
library(tidymodels)
library(apollo)
library(chatgpt)
library(dplyr)


setwd("~")

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

### Load Apollo library
library(apollo)

### Initialise code
apollo_initialise()

### Set core controls
apollo_control = list(
  modelName       = "Bangladesh CE",
  modelDescr      = "Mixed logit model on Bangladesh CE data",
  indivID         = "ID",  
  mixing          = TRUE,
  # workInLogs = TRUE,
  nCores          = 15,
  outputDirectory = "~"
)


des1 <- rotation.design(attribute.names = list(Labor_intensity = c ("L0_labor", "L1_labor", "L2_labor", "L3_labor"), Ecosystem_integration = c ("L0_ecosystem", "L1_ecosystem", "L2_ecosystem", "L3_ecosystem"), Government_Support = c ("L0_support", "L1_support", "L2_support", "L3_support"), Management_system = c ("L0_system", "L1_system", "L2_system", "L3_system"), Fingerling_cost = c ("500", "1000", "1500", "2000")), nalternatives = 2, nblocks = 1, row.renames = FALSE, randomize = TRUE, seed = 4592)



des2 <- rotation.design(attribute.names = list(Labor_intensity = c ("0", "1", "2", "3"), Ecosystem_integration = c ("0", "1", "2", "3"), Government_Support = c ("0", "1", "2", "3"), Management_system = c ("0", "1", "2", "3"), Fingerling_cost = c ("500", "1000", "1500", "2000")), nalternatives = 2, nblocks = 1, row.renames = FALSE, randomize = TRUE, seed =4592)

alt <- des2$alternatives  
alt1 <- alt$alt.1 
# Add "1" after variables' name
names(alt1)[names(alt1) == "Labor_intensity"] <- paste0("Labor_intensity1")
names(alt1)[names(alt1) == "Ecosystem_integration"] <- paste0("Ecosystem_integration1")
names(alt1)[names(alt1) == "Government_Support"] <- paste0("Government_Support1")
names(alt1)[names(alt1) == "Management_system"] <- paste0("Management_system1")
names(alt1)[names(alt1) == "Fingerling_cost"] <- paste0("Fingerling_cost1")

alt2 <- alt$alt.2
names(alt2)[names(alt2) == "Labor_intensity"] <- paste0("Labor_intensity2")
names(alt2)[names(alt2) == "Ecosystem_integration"] <- paste0("Ecosystem_integration2")
names(alt2)[names(alt2) == "Government_Support"] <- paste0("Government_Support2")
names(alt2)[names(alt2) == "Management_system"] <- paste0("Management_system2")
names(alt2)[names(alt2) == "Fingerling_cost"] <- paste0("Fingerling_cost2")

alt.all <- cbind(alt1,alt2)

alt.all <- alt.all[,-c(1,3,9:11)]


input1=read.csv("559.csv")

num <- length(unique(input1$ID))
alt.all.expand <- alt.all[rep(seq_len(nrow(alt.all)), times = num), ]

desmat1 <- make.design.matrix(choice.experiment.design = des1, optout = TRUE, categorical.attributes = c("Labor_intensity", "Ecosystem_integration", "Government_Support", "Management_system"), continuous.attributes = c("Fingerling_cost"), unlabeled = TRUE)

dataset <- make.dataset(respondent.dataset = input1, choice.indicators = c("q1", "q2", "q3", "q4", "q5", "q6", "q7", "q8", "q9", "q10", "q11", "q12", "q13", "q14", "q15", "q16"), design.matrix = desmat1)    
choice <- dataset[,c(1,3:5)]
choice <- subset(choice,RES=="TRUE")

database <- cbind(choice,alt.all.expand)
database$Fingerling_cost1 <- as.numeric(as.character(database$Fingerling_cost1))
database$Fingerling_cost2 <- as.numeric(as.character(database$Fingerling_cost2))


### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(asc = 0,
                sigma_asc = 0,
                #mu_Labor_intensity_L0    = 0,
                #sigma_Labor_intensity_L0 = 0,
                mu_Labor_intensity_L1    = 0,
                sigma_Labor_intensity_L1 = 0,
                mu_Labor_intensity_L2    = 0,
                sigma_Labor_intensity_L2 = 0,
                mu_Labor_intensity_L3    = 0,
                sigma_Labor_intensity_L3 = 0,
                #mu_Ecosystem_integration_L0    = 0,
                #sigma_Ecosystem_integration_L0 = 0,
                mu_Ecosystem_integration_L1    = 0,
                sigma_Ecosystem_integration_L1 = 0, 
                mu_Ecosystem_integration_L2    = 0,
                sigma_Ecosystem_integration_L2 = 0, 
                mu_Ecosystem_integration_L3    = 0,
                sigma_Ecosystem_integration_L3 = 0, 
                #mu_Government_Support_L0    = 0,
                #sigma_Government_Support_L0 = 0,
                mu_Government_Support_L1    = 0,
                sigma_Government_Support_L1 = 0,
                mu_Government_Support_L2    = 0,
                sigma_Government_Support_L2 = 0,
                mu_Government_Support_L3    = 0,
                sigma_Government_Support_L3 = 0,
                #mu_Management_system_L0    = 0,
                #sigma_Management_system_L0 = 0,
                mu_Management_system_L1    = 0,
                sigma_Management_system_L1 = 0,
                mu_Management_system_L2    = 0,
                sigma_Management_system_L2 = 0,
                mu_Management_system_L3    = 0,
                sigma_Management_system_L3 = 0,
                mu_log_Fingerling_cost    = -10,
                sigma_log_Fingerling_cost = 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()

# ################################################################# #
#### DEFINE RANDOM COMPONENTS                                    ####
# ################################################################# #

### Set parameters for generating draws
apollo_draws = list(
  interDrawsType = "mlhs",
  interNDraws    = 500,
  interUnifDraws = c(),
  interNormDraws = c("draws_Labor_intensity_L0","draws_Labor_intensity_L1","draws_Labor_intensity_L2","draws_Labor_intensity_L3","draws_Ecosystem_integration_L0","draws_Ecosystem_integration_L1","draws_Ecosystem_integration_L2","draws_Ecosystem_integration_L3","draws_Government_Support_L0","draws_Government_Support_L1","draws_Government_Support_L2","draws_Government_Support_L3","draws_Management_system_L0","draws_Management_system_L1","draws_Management_system_L2","draws_Management_system_L3","draws_Fingerling_cost","draws_asc"),
  intraDrawsType = "mlhs",
  intraNDraws    = 0,
  intraUnifDraws = c(),
  intraNormDraws = c()
)

### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
  randcoeff = list()
  # randcoeff[["b_Labor_intensity_L0"]] =mu_Labor_intensity_L0 + sigma_Labor_intensity_L0 * draws_Labor_intensity_L0 
  randcoeff[["b_Labor_intensity_L1"]] =mu_Labor_intensity_L1 + sigma_Labor_intensity_L1 * draws_Labor_intensity_L1 
  randcoeff[["b_Labor_intensity_L2"]] =mu_Labor_intensity_L2 + sigma_Labor_intensity_L2 * draws_Labor_intensity_L2 
  randcoeff[["b_Labor_intensity_L3"]] =mu_Labor_intensity_L3 + sigma_Labor_intensity_L3 * draws_Labor_intensity_L3 
  #  randcoeff[["b_Ecosystem_integration_L0"]] =  mu_Ecosystem_integration_L0 + sigma_Ecosystem_integration_L0 * draws_Ecosystem_integration_L0
  randcoeff[["b_Ecosystem_integration_L1"]] =  mu_Ecosystem_integration_L1 + sigma_Ecosystem_integration_L1 * draws_Ecosystem_integration_L1
  randcoeff[["b_Ecosystem_integration_L2"]] =  mu_Ecosystem_integration_L2 + sigma_Ecosystem_integration_L2 * draws_Ecosystem_integration_L2
  randcoeff[["b_Ecosystem_integration_L3"]] =  mu_Ecosystem_integration_L3 + sigma_Ecosystem_integration_L3 * draws_Ecosystem_integration_L3
  #  randcoeff[["b_Management_system_L0"]] =  mu_Management_system_L0 + sigma_Management_system_L0 * draws_Management_system_L0 
  randcoeff[["b_Management_system_L1"]] =  mu_Management_system_L1 + sigma_Management_system_L1 * draws_Management_system_L1 
  randcoeff[["b_Management_system_L2"]] =  mu_Management_system_L2 + sigma_Management_system_L2 * draws_Management_system_L2 
  randcoeff[["b_Management_system_L3"]] =  mu_Management_system_L3 + sigma_Management_system_L3 * draws_Management_system_L3 
  #  randcoeff[["b_Government_Support_L0"]] =  mu_Government_Support_L0 + sigma_Government_Support_L0 * draws_Government_Support_L0 
  randcoeff[["b_Government_Support_L1"]] =  mu_Government_Support_L1 + sigma_Government_Support_L1 * draws_Government_Support_L1 
  randcoeff[["b_Government_Support_L2"]] =  mu_Government_Support_L2 + sigma_Government_Support_L2 * draws_Government_Support_L2 
  randcoeff[["b_Government_Support_L3"]] =  mu_Government_Support_L3 + sigma_Government_Support_L3 * draws_Government_Support_L3 
  randcoeff[["b_Fingerling_cost"]] =  -exp(mu_log_Fingerling_cost + sigma_log_Fingerling_cost * draws_Fingerling_cost) 
  randcoeff[["ec"]] =  sigma_asc * draws_asc
  
  return(randcoeff)
}

# ################################################################# #
#### GROUP AND VALIDATE INPUTS                                   ####
# ################################################################# #

apollo_inputs = apollo_validateInputs()

# ################################################################# #
#### DEFINE MODEL AND LIKELIHOOD FUNCTION                        ####
# ################################################################# #

apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
  
  ### Function initialisation: do not ange the following three commands
  ### AareaaFingerling_cost inputs and detaFingerling_cost 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_seareaings, order is irrelevant
  V = list()
  V[['NB']]  = asc + ec
  V[['alt1']]  = #b_Labor_intensity_L0 * (Labor_intensity1==0) + 
    b_Labor_intensity_L1 * (Labor_intensity1==1) + b_Labor_intensity_L2 * (Labor_intensity1==2) +b_Labor_intensity_L3 * (Labor_intensity1==3) + 
    #b_Ecosystem_integration_L0 * (Ecosystem_integration1==0)+
    b_Ecosystem_integration_L1 * (Ecosystem_integration1==1) +  b_Ecosystem_integration_L2 * (Ecosystem_integration1==2)+ b_Ecosystem_integration_L3 * (Ecosystem_integration1==3) + 
    #b_Management_system_L0 * (Management_system1==0) +
    b_Management_system_L1 * (Management_system1==1) + b_Management_system_L2 * (Management_system1==2)+b_Management_system_L3 * (Management_system1==3) + 
    #b_Government_Support_L0 * 
    (Government_Support1==0) +b_Government_Support_L1 * (Government_Support1==1) + b_Government_Support_L2 * (Government_Support1==2)+b_Government_Support_L3 * (Government_Support1==3) + 
    b_Fingerling_cost * Fingerling_cost1 
  V[['alt2']]  = #b_Labor_intensity_L0 * (Labor_intensity2==0) + 
    b_Labor_intensity_L1 * (Labor_intensity2==1) + b_Labor_intensity_L2 * (Labor_intensity2==2) +b_Labor_intensity_L3 * (Labor_intensity2==3) + 
    #b_Ecosystem_integration_L0 * (Ecosystem_integration2==0) + 
    b_Ecosystem_integration_L1 * (Ecosystem_integration2==1) +  b_Ecosystem_integration_L2 * (Ecosystem_integration2==2)+ b_Ecosystem_integration_L3 * (Ecosystem_integration2==3) + 
    #b_Management_system_L0 * (Management_system2==0) + 
    b_Management_system_L1 * (Management_system2==1) + b_Management_system_L2 * (Management_system2==2)+b_Management_system_L3 * (Management_system2==3) + 
    #b_Government_Support_L0 * (Government_Support2==0) +
    b_Government_Support_L1 * (Government_Support2==1) + b_Government_Support_L2 * (Government_Support2==2)+b_Government_Support_L3 * (Government_Support2==3) + 
    b_Fingerling_cost * Fingerling_cost2
  
  
  ### Define seareaings for MNL model component
  mnl_seareaings = list(
    alternatives  = c(alt1="1", alt2="2", NB="3"),
    avail         = list(alt1=1, alt2=1, NB=1),
    choiceVar     = ALT,
    utilities     = V
  )
  
  ### Compute probabilities using MNL model
  P[["model"]] = apollo_mnl(mnl_seareaings, 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                                            ####
# ################################################################# #


model = apollo_estimate(apollo_beta, apollo_fixed,apollo_probabilities, apollo_inputs, estimate_settings = list(maxIterations=500))#,bootstrapSE=30))

# ################################################################# #
#### MODEL OUTPUTS                                               ####
# ################################################################# #

# ----------------------------------------------------------------- #
#---- FORMAareaED OUTPUT (TO SCREEN)                               ----
# ----------------------------------------------------------------- #

#apollo_loadModel("Bangladesh CE")

apollo_modelOutput(model, list(printPVal=2))
Results:

Code: Select all

stimation method                           : bfgs
Model diagnosis                             : successful convergence
Optimisation diagnosis                      : Maximum found
     hessian properties                     : Negative definitive
     maximum eigenvalue                     : 0
Number of individuals                       : 252
Number of rows in database                  : 4032
Number of modelled outcomes                 : 4032

Number of cores used                        :  15 
Number of inter-individual draws            : 500 (mlhs)

LL(start)                                   : -4464.3
LL at equal shares, LL(0)                   : -4429.6
LL at observed shares, LL(C)                : -3790.01
LL(final)                                   : -3195.03
Rho-squared vs equal shares                  :  0.2787 
Adj.Rho-squared vs equal shares              :  0.2724 
Rho-squared vs observed shares               :  0.157 
Adj.Rho-squared vs observed shares           :  0.1501 
AIC                                         :  6446.06 
BIC                                         :  6622.52 

Estimated parameters                        : 28
Time taken (hh:mm:ss)                       :  00:10:25.61 
     pre-estimation                         :  00:02:22.44 
     estimation                             :  00:03:22.77 
          initial estimation                :  00:03:15.28 
          estimation after rescaling        :  00:00:7.5 
     post-estimation                        :  00:04:40.4 
Iterations                                  :  102  
     initial estimation                     :  101 
     estimation after rescaling             :  1 

Unconstrained optimisation.

Estimates:
                                  Estimate        s.e.   t.rat.(0)  p(2-sided)    Rob.s.e. Rob.t.rat.(0)  p(2-sided)
asc                               -3.75892     0.37678    -9.97637    0.000000     0.44821       -8.3865    0.000000
sigma_asc                         -3.66013     0.37002    -9.89167    0.000000     0.43789       -8.3587    0.000000
mu_Labor_intensity_L1             -0.02455     0.12378    -0.19836    0.842764     0.11631       -0.2111    0.832808
sigma_Labor_intensity_L1          -1.10012     0.12803    -8.59265    0.000000     0.13879       -7.9267   2.220e-15
mu_Labor_intensity_L2             -0.28015     0.12110    -2.31336    0.020703     0.13049       -2.1468    0.031807
sigma_Labor_intensity_L2          -0.85888     0.15869    -5.41230   6.222e-08     0.17201       -4.9931   5.941e-07
mu_Labor_intensity_L3             -0.27549     0.07694    -3.58033  3.4316e-04     0.06771       -4.0685   4.732e-05
sigma_Labor_intensity_L3           0.07842     0.12872     0.60924    0.542365     0.07750        1.0119    0.311601
mu_Ecosystem_integration_L1        0.44106     0.12424     3.55012  3.8506e-04     0.13163        3.3507  8.0610e-04
sigma_Ecosystem_integration_L1     0.97570     0.10982     8.88421    0.000000     0.12245        7.9680   1.554e-15
mu_Ecosystem_integration_L2        0.11479     0.09728     1.17992    0.238031     0.11329        1.0132    0.310965
sigma_Ecosystem_integration_L2    -0.07574     0.13440    -0.56357    0.573050     0.09016       -0.8401    0.400842
mu_Ecosystem_integration_L3        0.12817     0.09487     1.35098    0.176702     0.10943        1.1713    0.241488
sigma_Ecosystem_integration_L3     0.51177     0.12595     4.06317   4.841e-05     0.12536        4.0825   4.455e-05
mu_Government_Support_L1           0.48127     0.09919     4.85209   1.222e-06     0.09484        5.0744   3.887e-07
sigma_Government_Support_L1       -0.16938     0.17465    -0.96985    0.332122     0.12701       -1.3336    0.182329
mu_Government_Support_L2           0.68296     0.10267     6.65202   2.891e-11     0.09940        6.8707   6.391e-12
sigma_Government_Support_L2       -0.87376     0.11951    -7.31151   2.642e-13     0.11870       -7.3613   1.823e-13
mu_Government_Support_L3           0.23872     0.10403     2.29486    0.021741     0.10881        2.1940    0.028238
sigma_Government_Support_L3        0.62242     0.14810     4.20272   2.637e-05     0.15362        4.0517   5.085e-05
mu_Management_system_L1           -0.37760     0.12265    -3.07866    0.002079     0.12372       -3.0521    0.002273
sigma_Management_system_L1        -1.26939     0.12358   -10.27198    0.000000     0.11781      -10.7750    0.000000
mu_Management_system_L2           -0.30828     0.13379    -2.30426    0.021208     0.13928       -2.2135    0.026865
sigma_Management_system_L2        -1.16989     0.11492   -10.18017    0.000000     0.11515      -10.1593    0.000000
mu_Management_system_L3            0.09415     0.08764     1.07427    0.282702     0.08782        1.0720    0.283712
sigma_Management_system_L3        -0.07113     0.12247    -0.58076    0.561405     0.09719       -0.7319    0.464248
mu_log_Fingerling_cost         -1066.11898   1.484e+04    -0.07183    0.942736   205.05933       -5.1991   2.003e-07
sigma_log_Fingerling_cost        386.32047  5409.26766     0.07142    0.943065    74.73679        5.1691   2.352e-07

Re: Unreasonable estimate for cost

Posted: 09 Jun 2023, 11:29
by stephanehess
Hi

this would imply that the median is very close to zero. Did you estimate an MNL model - what were the result like for that for this coefficient?

Stephane

Re: Unreasonable estimate for cost

Posted: 13 Jun 2023, 02:25
by Nishihara
Dear Professor,

Thank you very much for your reply.

I tried to do the estimation through MNL but come out this error when I run the apollo_estimate:

"Error in apollo_probabilities(apollo_beta, apollo_inputs, functionality = "validate") :
object '"validate"' not found".

I attached my code below. Thank you so much for your kind help.

Code: Select all

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

### Load Apollo library
library(apollo)

### Initialise code
apollo_initialise()

### Set core controls
apollo_control = list(
  modelName       = "Bangladesh CE",
  modelDescr      = "Mixed logit model on Bangladesh CE data",
  indivID         = "ID",  
  nCores          = 15,
  outputDirectory = "~"
)


des1 <- rotation.design(attribute.names = list(Labor_intensity = c ("L0_labor", "L1_labor", "L2_labor", "L3_labor"), Ecosystem_integration = c ("L0_ecosystem", "L1_ecosystem", "L2_ecosystem", "L3_ecosystem"), Government_Support = c ("L0_support", "L1_support", "L2_support", "L3_support"), Management_system = c ("L0_system", "L1_system", "L2_system", "L3_system"), Fingerling_cost = c ("500", "1000", "1500", "2000")), nalternatives = 2, nblocks = 1, row.renames = FALSE, randomize = TRUE, seed = 4592)

des2 <- rotation.design(attribute.names = list(Labor_intensity = c ("0", "1", "2", "3"), Ecosystem_integration = c ("0", "1", "2", "3"), Government_Support = c ("0", "1", "2", "3"), Management_system = c ("0", "1", "2", "3"), Fingerling_cost = c ("500", "1000", "1500", "2000")), nalternatives = 2, nblocks = 1, row.renames = FALSE, randomize = TRUE, seed =4592)

alt <- des2$alternatives  
alt1 <- alt$alt.1 
# Add "1" after variables' name
names(alt1)[names(alt1) == "Labor_intensity"] <- paste0("Labor_intensity1")
names(alt1)[names(alt1) == "Ecosystem_integration"] <- paste0("Ecosystem_integration1")
names(alt1)[names(alt1) == "Government_Support"] <- paste0("Government_Support1")
names(alt1)[names(alt1) == "Management_system"] <- paste0("Management_system1")
names(alt1)[names(alt1) == "Fingerling_cost"] <- paste0("Fingerling_cost1")

alt2 <- alt$alt.2
names(alt2)[names(alt2) == "Labor_intensity"] <- paste0("Labor_intensity2")
names(alt2)[names(alt2) == "Ecosystem_integration"] <- paste0("Ecosystem_integration2")
names(alt2)[names(alt2) == "Government_Support"] <- paste0("Government_Support2")
names(alt2)[names(alt2) == "Management_system"] <- paste0("Management_system2")
names(alt2)[names(alt2) == "Fingerling_cost"] <- paste0("Fingerling_cost2")

alt.all <- cbind(alt1,alt2)

alt.all <- alt.all[,-c(1,3,9:11)]

input1=read.csv("559.csv")

num <- length(unique(input1$ID))
alt.all.expand <- alt.all[rep(seq_len(nrow(alt.all)), times = num), ]

desmat1 <- make.design.matrix(choice.experiment.design = des1, optout = TRUE, categorical.attributes = c("Labor_intensity", "Ecosystem_integration", "Government_Support", "Management_system"), continuous.attributes = c("Fingerling_cost"), unlabeled = TRUE)

dataset <- make.dataset(respondent.dataset = input1, choice.indicators = c("q1", "q2", "q3", "q4", "q5", "q6", "q7", "q8", "q9", "q10", "q11", "q12", "q13", "q14", "q15", "q16"), design.matrix = desmat1)    
choice <- dataset[,c(1,3:5)]
choice <- subset(choice,RES=="TRUE")

database <- cbind(choice,alt.all.expand)
database$Fingerling_cost1 <- as.numeric(as.character(database$Fingerling_cost1))
database$Fingerling_cost2 <- as.numeric(as.character(database$Fingerling_cost2))


### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(
                asc = 0,
                b_Labor_intensity_L1    = 0,
                b_Labor_intensity_L2    = 0,
                b_Labor_intensity_L3    = 0,
                b_Ecosystem_integration_L1    = 0,
                b_Ecosystem_integration_L2    = 0,
                b_Ecosystem_integration_L3    = 0,
                b_Government_Support_L1    = 0,
                b_Government_Support_L2    = 0,
                b_Government_Support_L3    = 0,
                b_Management_system_L1    = 0,
                b_Management_system_L2    = 0,
                b_Management_system_L3    = 0,
                b_Fingerling_cost    = 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()


# ################################################################# #
#### GROUP AND VALIDATE INPUTS                                   ####
# ################################################################# #

apollo_inputs = apollo_validateInputs()

# ################################################################# #
#### DEFINE MODEL AND LIKELIHOOD FUNCTION                        ####
# ################################################################# #

apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
  
  ### Function initialisation: do not change the following three commands
  ### AareaaFingerling_cost inputs and detaFingerling_cost 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 bst use the same names as in mnl_seareaings, order is irrelevant
  V = list()
  V[['NB']]  = asc 
  V[['alt1']]  =  
    b_Labor_intensity_L1 * (Labor_intensity1==1) + b_Labor_intensity_L2 * (Labor_intensity1==2) +b_Labor_intensity_L3 * (Labor_intensity1==3) + 
    b_Ecosystem_integration_L1 * (Ecosystem_integration1==1) +  b_Ecosystem_integration_L2 * (Ecosystem_integration1==2)+ b_Ecosystem_integration_L3 * (Ecosystem_integration1==3) + 
    b_Management_system_L1 * (Management_system1==1) + b_Management_system_L2 * (Management_system1==2)+b_Management_system_L3 * (Management_system1==3) + 
    (Government_Support1==0) +b_Government_Support_L1 * (Government_Support1==1) + b_Government_Support_L2 * (Government_Support1==2)+b_Government_Support_L3 * (Government_Support1==3) + 
    b_Fingerling_cost * Fingerling_cost1 
  V[['alt2']]  = 
    b_Labor_intensity_L1 * (Labor_intensity2==1) + b_Labor_intensity_L2 * (Labor_intensity2==2) +b_Labor_intensity_L3 * (Labor_intensity2==3) + 
    b_Ecosystem_integration_L1 * (Ecosystem_integration2==1) +  b_Ecosystem_integration_L2 * (Ecosystem_integration2==2)+ b_Ecosystem_integration_L3 * (Ecosystem_integration2==3) + 
    b_Management_system_L1 * (Management_system2==1) + b_Management_system_L2 * (Management_system2==2)+b_Management_system_L3 * (Management_system2==3) + 
    b_Government_Support_L1 * (Government_Support2==1) + b_Government_Support_L2 * (Government_Support2==2)+b_Government_Support_L3 * (Government_Support2==3) + 
    b_Fingerling_cost * Fingerling_cost2
  
  
  ### Define seareaings for MNL model component
  mnl_seareaings = list(
    alternatives  = c(alt1=1, alt2=2, NB=3),
    avail         = list(alt1=1, alt2=1, NB=1),
    choiceVar     = ALT,
    utilities     = 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)
  
  ### 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)#

Re: Unreasonable estimate for cost

Posted: 15 Jun 2023, 13:33
by stephanehess
do you get any error messages before that? Also, what version of Apollo are you using?

Re: Unreasonable estimate for cost

Posted: 16 Jun 2023, 03:12
by Nishihara
Dear Professor,

Thank you very much for your help.

No, I did not get any error before apollo_estimate, and I am using the latest version of 0.2.9.

Re: Unreasonable estimate for cost

Posted: 17 Jun 2023, 08:55
by stephanehess
Hi

are you by any chance running R4.3.0 or another very recent version. The R team made a change that led to this rather unhelpful error message. If you upgrade to R4.3.1, you will get the right messages about your error

Stephane

Re: Unreasonable estimate for cost

Posted: 18 Jun 2023, 08:38
by Nishihara
Dear Professor,

Thank you very much for your help.

It did work after I update to R4.3.1 and no error came out.

However, through MNL, the estimate for the cost is positive but the value is reasonable. I attached the code and results below.

Do you think there is a connection between the unreasonable estimate for cost from the Mixed logit model and the result from MNL?

Thank you very much in advance.

Code: Select all

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

### Load Apollo library
library(apollo)

### Initialise code
apollo_initialise()

### Set core controls
apollo_control = list(
  modelName       = "Bangladesh CE",
  modelDescr      = "MNL on Bangladesh CE data",
  indivID         = "ID",  
  nCores          = 15,
  outputDirectory = "~"
)


des1 <- rotation.design(attribute.names = list(Labor_intensity = c ("L0_labor", "L1_labor", "L2_labor", "L3_labor"), Ecosystem_integration = c ("L0_ecosystem", "L1_ecosystem", "L2_ecosystem", "L3_ecosystem"), Government_Support = c ("L0_support", "L1_support", "L2_support", "L3_support"), Management_system = c ("L0_system", "L1_system", "L2_system", "L3_system"), Fingerling_cost = c ("500", "1000", "1500", "2000")), nalternatives = 2, nblocks = 1, row.renames = FALSE, randomize = TRUE, seed = 4592)



des2 <- rotation.design(attribute.names = list(Labor_intensity = c ("0", "1", "2", "3"), Ecosystem_integration = c ("0", "1", "2", "3"), Government_Support = c ("0", "1", "2", "3"), Management_system = c ("0", "1", "2", "3"), Fingerling_cost = c ("500", "1000", "1500", "2000")), nalternatives = 2, nblocks = 1, row.renames = FALSE, randomize = TRUE, seed =4592)

alt <- des2$alternatives  
alt1 <- alt$alt.1 
# Add "1" after variables' name
names(alt1)[names(alt1) == "Labor_intensity"] <- paste0("Labor_intensity1")
names(alt1)[names(alt1) == "Ecosystem_integration"] <- paste0("Ecosystem_integration1")
names(alt1)[names(alt1) == "Government_Support"] <- paste0("Government_Support1")
names(alt1)[names(alt1) == "Management_system"] <- paste0("Management_system1")
names(alt1)[names(alt1) == "Fingerling_cost"] <- paste0("Fingerling_cost1")

alt2 <- alt$alt.2
names(alt2)[names(alt2) == "Labor_intensity"] <- paste0("Labor_intensity2")
names(alt2)[names(alt2) == "Ecosystem_integration"] <- paste0("Ecosystem_integration2")
names(alt2)[names(alt2) == "Government_Support"] <- paste0("Government_Support2")
names(alt2)[names(alt2) == "Management_system"] <- paste0("Management_system2")
names(alt2)[names(alt2) == "Fingerling_cost"] <- paste0("Fingerling_cost2")

alt.all <- cbind(alt1,alt2)

alt.all <- alt.all[,-c(1,3,9:11)]


input1=read.csv("559.csv")

num <- length(unique(input1$ID))
alt.all.expand <- alt.all[rep(seq_len(nrow(alt.all)), times = num), ]

desmat1 <- make.design.matrix(choice.experiment.design = des1, optout = TRUE, categorical.attributes = c("Labor_intensity", "Ecosystem_integration", "Government_Support", "Management_system"), continuous.attributes = c("Fingerling_cost"), unlabeled = TRUE)

dataset <- make.dataset(respondent.dataset = input1, choice.indicators = c("q1", "q2", "q3", "q4", "q5", "q6", "q7", "q8", "q9", "q10", "q11", "q12", "q13", "q14", "q15", "q16"), design.matrix = desmat1)    
choice <- dataset[,c(1,3:5)]
choice <- subset(choice,RES=="TRUE")

database <- cbind(choice,alt.all.expand)
database$Fingerling_cost1 <- as.numeric(as.character(database$Fingerling_cost1))
database$Fingerling_cost2 <- as.numeric(as.character(database$Fingerling_cost2))


### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(#asc_1 = 0,
                #asc_2 = 0,
                asc = 0,
                b_Labor_intensity_L1    = 0,
                b_Labor_intensity_L2    = 0,
                b_Labor_intensity_L3    = 0,
                b_Ecosystem_integration_L1    = 0,
                b_Ecosystem_integration_L2    = 0,
                b_Ecosystem_integration_L3    = 0,
                b_Government_Support_L1    = 0,
                b_Government_Support_L2    = 0,
                b_Government_Support_L3    = 0,
                b_Management_system_L1    = 0,
                b_Management_system_L2    = 0,
                b_Management_system_L3    = 0,
                b_Fingerling_cost    = 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()


# ################################################################# #
#### GROUP AND VALIDATE INPUTS                                   ####
# ################################################################# #

apollo_inputs = apollo_validateInputs()

# ################################################################# #
#### DEFINE MODEL AND LIKELIHOOD FUNCTION                        ####
# ################################################################# #

apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
  
  ### Function initialisation: do not ange the following three commands
  ### AareaaFingerling_cost inputs and detaFingerling_cost 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 bst use the same names as in mnl_seareaings, order is irrelevant
  V = list()
  V[['NB']]  = asc 
  V[['alt1']]  =  
    b_Labor_intensity_L1 * (Labor_intensity1==1) + b_Labor_intensity_L2 * (Labor_intensity1==2) +b_Labor_intensity_L3 * (Labor_intensity1==3) + 
    b_Ecosystem_integration_L1 * (Ecosystem_integration1==1) +  b_Ecosystem_integration_L2 * (Ecosystem_integration1==2)+ b_Ecosystem_integration_L3 * (Ecosystem_integration1==3) + 
    b_Management_system_L1 * (Management_system1==1) + b_Management_system_L2 * (Management_system1==2)+b_Management_system_L3 * (Management_system1==3) + 
    (Government_Support1==0) +b_Government_Support_L1 * (Government_Support1==1) + b_Government_Support_L2 * (Government_Support1==2)+b_Government_Support_L3 * (Government_Support1==3) + 
    b_Fingerling_cost * Fingerling_cost1 
  V[['alt2']]  = 
    b_Labor_intensity_L1 * (Labor_intensity2==1) + b_Labor_intensity_L2 * (Labor_intensity2==2) +b_Labor_intensity_L3 * (Labor_intensity2==3) + 
    b_Ecosystem_integration_L1 * (Ecosystem_integration2==1) +  b_Ecosystem_integration_L2 * (Ecosystem_integration2==2)+ b_Ecosystem_integration_L3 * (Ecosystem_integration2==3) + 
    b_Management_system_L1 * (Management_system2==1) + b_Management_system_L2 * (Management_system2==2)+b_Management_system_L3 * (Management_system2==3) + 
    b_Government_Support_L1 * (Government_Support2==1) + b_Government_Support_L2 * (Government_Support2==2)+b_Government_Support_L3 * (Government_Support2==3) + 
    b_Fingerling_cost * Fingerling_cost2
  
  
  ### Define seareaings for MNL model component
  mnl_settings = list(
    alternatives  = c(alt1=1, alt2=2, NB=3),
    avail         = list(alt1=1, alt2=1, NB=1),
    choiceVar     = ALT,
    utilities     = 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)
  
  ### 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)#, estimate_settings = list(maxIterations=500))#,bootstrapSE=30))

# ################################################################# #
#### MODEL OUTPUTS                                               ####
# ################################################################# #

# ----------------------------------------------------------------- #
#---- FORMAareaED OUTPUT (TO SCREEN)                               ----
# ----------------------------------------------------------------- #

#apollo_loadModel("Bangladesh CE")

apollo_modelOutput(model, list(printPVal=2))

Results:

Code: Select all

Model run by Jie using Apollo 0.2.9 on R 4.3.1 for Darwin.
www.ApolloChoiceModelling.com

Model name                                  : Bangladesh CE
Model description                           : MNL on Bangladesh CE data
Model run at                                : 2023-06-18 15:30:06.251488
Estimation method                           : bfgs
Model diagnosis                             : successful convergence
Optimisation diagnosis                      : Maximum found
     hessian properties                     : Negative definitive
     maximum eigenvalue                     : -33.620427
Number of individuals                       : 252
Number of rows in database                  : 4032
Number of modelled outcomes                 : 4032

Number of cores used                        :  15 
Model without mixing

LL(start)                                   : -4399.06
LL at equal shares, LL(0)                   : -4429.6
LL at observed shares, LL(C)                : -3790.01
LL(final)                                   : -3746.09
Rho-squared vs equal shares                  :  0.1543 
Adj.Rho-squared vs equal shares              :  0.1511 
Rho-squared vs observed shares               :  0.0116 
Adj.Rho-squared vs observed shares           :  0.0084 
AIC                                         :  7520.18 
BIC                                         :  7608.4 

Estimated parameters                        : 14
Time taken (hh:mm:ss)                       :  00:00:33.63 
     pre-estimation                         :  00:00:30.38 
     estimation                             :  00:00:1.57 
          initial estimation                :  00:00:1.48 
          estimation after rescaling        :  00:00:0.1 
     post-estimation                        :  00:00:1.67 
Iterations                                  :  24  
     initial estimation                     :  23 
     estimation after rescaling             :  1 

Unconstrained optimisation.

Estimates:
                              Estimate        s.e.   t.rat.(0)  p(2-sided)    Rob.s.e. Rob.t.rat.(0)  p(2-sided)
asc                           -0.87549     0.13369   -6.548407   5.815e-11     0.17704     -4.945173   7.608e-07
b_Labor_intensity_L1           0.02344     0.08044    0.291366    0.770772     0.07137      0.328369    0.742632
b_Labor_intensity_L2          -0.24491     0.08528   -2.871874    0.004080     0.07777     -3.149327    0.001636
b_Labor_intensity_L3          -0.19265     0.06766   -2.847238    0.004410     0.05646     -3.411967  6.4496e-04
b_Ecosystem_integration_L1     0.32801     0.08401    3.904426   9.445e-05     0.08495      3.861238  1.1281e-04
b_Ecosystem_integration_L2  6.9484e-04     0.08003    0.008682    0.993073     0.08431      0.008242    0.993424
b_Ecosystem_integration_L3     0.05349     0.07349    0.727858    0.466700     0.07971      0.671022    0.502206
b_Government_Support_L1        0.55402     0.07804    7.099463   1.252e-12     0.06418      8.631818    0.000000
b_Government_Support_L2        0.73521     0.07047   10.433013    0.000000     0.06987     10.522770    0.000000
b_Government_Support_L3        0.35307     0.07868    4.487211   7.216e-06     0.07989      4.419425   9.896e-06
b_Management_system_L1        -0.24548     0.07947   -3.089073    0.002008     0.08585     -2.859529    0.004243
b_Management_system_L2        -0.23672     0.08649   -2.736802    0.006204     0.08010     -2.955154    0.003125
b_Management_system_L3         0.06865     0.07039    0.975334    0.329395     0.05913      1.160980    0.245650
b_Fingerling_cost           2.1246e-04   5.348e-05    3.972637   7.108e-05   4.337e-05      4.898439   9.660e-07

Re: Unreasonable estimate for cost

Posted: 18 Jun 2023, 08:45
by stephanehess
Hi

you say the estimate is positive but the value is reasonable. How can a positive cost coefficient be reasoable?

If the estimate is positive in MNL, then there is clearly a problem, and if you then use a negative lognormal, it will collapse to close to zero. This is why should always start with MNL and try to understand what is happening with your data.

Also, you have 3 alternatives and should estimat 2 constants, not just 1 - there may be order effects for example

Stephane

Re: Unreasonable estimate for cost

Posted: 18 Jun 2023, 09:03
by Nishihara
Dear Professor,

Thank you very much for your prompt reply.

Now I understand the problem.

May I ask if I have results with positive values of cost estimate from MNL, what should I do with the analysis?

Does it mean that the data is problematic and the analysis is undoable?

Thank you very much for your help.

Re: Unreasonable estimate for cost

Posted: 18 Jun 2023, 09:48
by stephanehess
Difficult to say as it really depends on your data. I suggest you study your data before you do any modelling and understand whether people really choose the options more when they are more expensive as that's what your results seem to tell you
All times are UTC+01:00
Page 1 of 2

Powered by phpBB® Forum Software © phpBB Limited