I want to use a mixed logit model to do estimation. Below is my code and I got the error:
Testing influence of parameters......Error in apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, :
Parameter mu_Efficiency_L2 does not influence the log-likelihood of your model!
Code: Select all
library(readr)
library(support.CEs)
library(survival)
library(openxlsx)
library(tidyr)
library(tidymodels)
library(apollo)
library(chatgpt)
library(dplyr)
setwd("C:/Users/sujie/OneDrive - The University of Tokyo/WorldFish/Bangladesh characterization study/DCE/DCE/v2/CSVs")
### 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,
nCores = 7
)
des1 <- rotation.design(attribute.names = list(Efficiency = c ("L0", "L1_breeding", "L2_breeding", "L3_breeding"), Productivity = c ("L0", "L1_productivity", "L2_productivity", "L3_productivity"), Productivity_characteristics = c ("L0", "L1_traits", "L2_traits", "L3_traits"), Pond_management = c ("L0", "L1_management", "L2_management", "L3_management"), Fingerling_cost = c ("500", "1000", "1500", "2000")), nalternatives = 2, nblocks = 1, row.renames = FALSE, randomize = TRUE, seed =6965)
des2 <- rotation.design(attribute.names = list(Efficiency = c ("0", "1", "2", "3"), Productivity = c ("0", "1", "2", "3"), Productivity_characteristics = c ("0", "1", "2", "3"), Pond_management = c ("0", "1", "2", "3"), Fingerling_cost = c ("500", "1000", "1500", "2000")), nalternatives = 2, nblocks = 1, row.renames = FALSE, randomize = TRUE, seed =6965)
alt <- des2$alternatives
alt1 <- alt$alt.1
# Add "1" after variables' name
names(alt1)[names(alt1) == "Efficiency"] <- paste0("Efficiency1")
names(alt1)[names(alt1) == "Productivity"] <- paste0("Productivity1")
names(alt1)[names(alt1) == "Productivity_characteristics"] <- paste0("Productivity_characteristics1")
names(alt1)[names(alt1) == "Pond_management"] <- paste0("Pond_management1")
names(alt1)[names(alt1) == "Fingerling_cost"] <- paste0("Fingerling_cost1")
alt2 <- alt$alt.2
names(alt2)[names(alt2) == "Efficiency"] <- paste0("Efficiency2")
names(alt2)[names(alt2) == "Productivity"] <- paste0("Productivity2")
names(alt2)[names(alt2) == "Productivity_characteristics"] <- paste0("Productivity_characteristics2")
names(alt2)[names(alt2) == "Pond_management"] <- paste0("Pond_management2")
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("1.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("Efficiency", "Productivity", "Productivity_characteristics", "Pond_management"), 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_Efficiency_L0 = 0,
sigma_Efficiency_L0 = 0,
mu_Efficiency_L1 = 0,
sigma_Efficiency_L1 = 0,
mu_Efficiency_L2 = 0,
sigma_Efficiency_L2 = 0,
mu_Efficiency_L3 = 0,
sigma_Efficiency_L3 = 0,
mu_Productivity_L0 = 0,
sigma_Productivity_L0 = 0,
mu_Productivity_L1 = 0,
sigma_Productivity_L1 = 0,
mu_Productivity_L2 = 0,
sigma_Productivity_L2 = 0,
mu_Productivity_L3 = 0,
sigma_Productivity_L3 = 0,
mu_Productivity_characteristics_L0 = 0,
sigma_Productivity_characteristics_L0 = 0,
mu_Productivity_characteristics_L1 = 0,
sigma_Productivity_characteristics_L1 = 0,
mu_Productivity_characteristics_L2 = 0,
sigma_Productivity_characteristics_L2 = 0,
mu_Productivity_characteristics_L3 = 0,
sigma_Productivity_characteristics_L3 = 0,
mu_Pond_management_L0 = 0,
sigma_Pond_management_L0 = 0,
mu_Pond_management_L1 = 0,
sigma_Pond_management_L1 = 0,
mu_Pond_management_L2 = 0,
sigma_Pond_management_L2 = 0,
mu_Pond_management_L3 = 0,
sigma_Pond_management_L3 = 0,
mu_Fingerling_cost = 0,
sigma_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()
# ################################################################# #
#### DEFINE RANDOM COMPONENTS ####
# ################################################################# #
### Set parameters for generating draws
apollo_draws = list(
interDrawsType = "mlhs",
interNDraws = 100,
interUnifDraws = c(),
interNormDraws = c("draws_Efficiency_L0","draws_Efficiency_L1","draws_Efficiency_L2","draws_Efficiency_L3","draws_Productivity_L0","draws_Productivity_L1","draws_Productivity_L2","draws_Productivity_L3","draws_Productivity_characteristics_L0","draws_Productivity_characteristics_L1","draws_Productivity_characteristics_L2","draws_Productivity_characteristics_L3","draws_Pond_management_L0","draws_Pond_management_L1","draws_Pond_management_L2","draws_Pond_management_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_Efficiency_L0"]] =mu_Efficiency_L0 + sigma_Efficiency_L0 * draws_Efficiency_L0
randcoeff[["b_Efficiency_L1"]] =mu_Efficiency_L1 + sigma_Efficiency_L1 * draws_Efficiency_L1
randcoeff[["b_Efficiency_L2"]] =mu_Efficiency_L2 + sigma_Efficiency_L2 * draws_Efficiency_L2
randcoeff[["b_Efficiency_L3"]] =mu_Efficiency_L3 + sigma_Efficiency_L3 * draws_Efficiency_L3
randcoeff[["b_Productivity_L0"]] = mu_Productivity_L0 + sigma_Productivity_L0 * draws_Productivity_L0
randcoeff[["b_Productivity_L1"]] = mu_Productivity_L1 + sigma_Productivity_L1 * draws_Productivity_L1
randcoeff[["b_Productivity_L2"]] = mu_Productivity_L2 + sigma_Productivity_L2 * draws_Productivity_L2
randcoeff[["b_Productivity_L3"]] = mu_Productivity_L3 + sigma_Productivity_L3 * draws_Productivity_L3
randcoeff[["b_Pond_management_L0"]] = mu_Pond_management_L0 + sigma_Pond_management_L0 * draws_Pond_management_L0
randcoeff[["b_Pond_management_L1"]] = mu_Pond_management_L1 + sigma_Pond_management_L1 * draws_Pond_management_L1
randcoeff[["b_Pond_management_L2"]] = mu_Pond_management_L2 + sigma_Pond_management_L2 * draws_Pond_management_L2
randcoeff[["b_Pond_management_L3"]] = mu_Pond_management_L3 + sigma_Pond_management_L3 * draws_Pond_management_L3
randcoeff[["b_Productivity_characteristics_L0"]] = mu_Productivity_characteristics_L0 + sigma_Productivity_characteristics_L0 * draws_Productivity_characteristics_L0
randcoeff[["b_Productivity_characteristics_L1"]] = mu_Productivity_characteristics_L1 + sigma_Productivity_characteristics_L1 * draws_Productivity_characteristics_L1
randcoeff[["b_Productivity_characteristics_L2"]] = mu_Productivity_characteristics_L2 + sigma_Productivity_characteristics_L2 * draws_Productivity_characteristics_L2
randcoeff[["b_Productivity_characteristics_L3"]] = mu_Productivity_characteristics_L3 + sigma_Productivity_characteristics_L3 * draws_Productivity_characteristics_L3
randcoeff[["b_Fingerling_cost"]] = mu_Fingerling_cost + sigma_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_Efficiency_L0 * (Efficiency1==0) + b_Efficiency_L1 * (Efficiency1==1) + b_Efficiency_L1 * (Efficiency1==2) +b_Efficiency_L1 * (Efficiency1==3) +
b_Productivity_L0 * (Productivity1==0)+b_Productivity_L1 * (Productivity1==1) + b_Productivity_L2 * (Productivity1==2)+ b_Productivity_L3 * (Productivity1==3) +
b_Pond_management_L0 * (Pond_management1==0) +b_Pond_management_L1 * (Pond_management1==1) + b_Pond_management_L2 * (Pond_management1==2)+b_Pond_management_L3 * (Pond_management1==3) +
b_Productivity_characteristics_L0 * (Productivity_characteristics1==0) +b_Productivity_characteristics_L1 * (Productivity_characteristics1==1) + b_Productivity_characteristics_L2 * (Productivity_characteristics1==2)+b_Productivity_characteristics_L3 * (Productivity_characteristics1==3) + b_Fingerling_cost * Fingerling_cost1
V[['alt2']] = b_Efficiency_L0 * (Efficiency2==0) + b_Efficiency_L1 * (Efficiency2==1) + b_Efficiency_L1 * (Efficiency2==2) +b_Efficiency_L1 * (Efficiency2==3) +
b_Productivity_L0 * (Productivity2==0) + b_Productivity_L1 * (Productivity2==1) + b_Productivity_L2 * (Productivity2==2)+ b_Productivity_L3 * (Productivity2==3) +
b_Pond_management_L0 * (Pond_management2==0) + b_Pond_management_L1 * (Pond_management2==1) + b_Pond_management_L2 * (Pond_management2==2)+b_Pond_management_L3 * (Pond_management2==3) +
b_Productivity_characteristics_L0 * (Productivity_characteristics2==0) +b_Productivity_characteristics_L1 * (Productivity_characteristics2==1) + b_Productivity_characteristics_L2 * (Productivity_characteristics2==2)+b_Productivity_characteristics_L3 * (Productivity_characteristics2==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)
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMAareaED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)
Best regards,
Joe