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))
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