Compare the log-likelihood values produced by Apollo, mixl, and logitr using the same set of parameter values.

[xindou-cyber]

Hi @jhelvy ,

First of all, thank you very much for the logitr package. The functionality for exploring multiple starting values in WTP space is especially insightful and very nicely implemented.

While trying to better understand how logitr works, I ran into an interesting issue. Using the first example under “WTP space convergence issues in other packages”, the final log-likelihood reported by logitr is -732.2132421.

I then took the corresponding parameter estimates from this model and plugged them into Apollo and the mixl , without re-estimating the model, only evaluating the log-likelihood at those parameter values. In both Apollo and mixl, the log-likelihood at that parameter vector is -1009.453, which is different from the value reported by logitr.

Could you please help me understand what might be causing this difference?

For reference, I have attached the Apollo and mixl code used to compute the log-likelihood at the same parameter values.

Thank you very much for your time and for developing such a helpful package.

Best regards,
Xin

`library(logitr)
library(mlogit)
library(gmnl)
library(apollo)
library(mixl)
library(dplyr)
library(tidyr)

set.seed(1234)

numDraws_wtp <- 50
start_wtp <- c(
scalePar = 0.388098,
feat = 1.738832,
brandhiland = -13.536714,
brandweight = -4.510068,
brandyoplait = 6.663869,
sd_feat = 0.491225,
sd_brandhiland = 5.730571,
sd_brandweight = 11.420500,
sd_brandyoplait = 8.872470
)

yogurt <- subset(logitr::yogurt, logitr::yogurt$id <= 50)

yogurt_price <- yogurt %>%
select(id, obsID, price, brand) %>%
mutate(price = -1*price) %>%
pivot_wider(
names_from = 'brand',
values_from = 'price') %>%
rename(
price_dannon = dannon,
price_hiland = hiland,
price_weight = weight,
price_yoplait = yoplait)
yogurt_feat <- yogurt %>%
select(id, obsID, feat, brand) %>%
pivot_wider(
names_from = 'brand',
values_from = 'feat') %>%
rename(
feat_dannon = dannon,
feat_hiland = hiland,
feat_weight = weight,
feat_yoplait = yoplait)
yogurt_choice <- yogurt %>%
filter(choice == 1) %>%
select(id, obsID, choice = alt)
data_apollo <- yogurt_price %>%
left_join(yogurt_feat, by = c('id', 'obsID')) %>%
left_join(yogurt_choice, by = c('id', 'obsID')) %>%
arrange(id, obsID) %>%
mutate(
av_dannon = 1,
av_hiland = 1,
av_weight = 1,
av_yoplait = 1
)

apollo_probabilities_wtp <- function(
apollo_beta, apollo_inputs, functionality = "estimate"
) {

apollo_attach(apollo_beta, apollo_inputs)
on.exit(apollo_detach(apollo_beta, apollo_inputs))

P <- list()

V <- list()
V[["dannon"]] <- scalePar * (b_feat * feat_dannon - price_dannon)
V[["hiland"]] <- scalePar * (b_brandhiland + b_feat * feat_hiland - price_hiland)
V[["weight"]] <- scalePar * (b_brandweight + b_feat * feat_weight - price_weight)
V[["yoplait"]] <- scalePar * (b_brandyoplait + b_feat * feat_yoplait - price_yoplait)

mnl_settings <- list(
alternatives = c(dannon = 1, hiland = 2, weight = 3, yoplait = 4),
avail = list(
dannon = av_dannon,
hiland = av_hiland,
weight = av_weight,
yoplait = av_yoplait),
choiceVar = choice,
utilities = V
)

P[["model"]] <- apollo_mnl(mnl_settings, functionality)

P <- apollo_panelProd(P, apollo_inputs, functionality)

P = apollo_avgInterDraws(P, apollo_inputs, functionality)

P <- apollo_prepareProb(P, apollo_inputs, functionality)

return(P)
}

apollo_randCoeff <- function(apollo_beta, apollo_inputs) {
randcoeff <- list()
randcoeff[['b_feat']] <- feat + d_feat * sd_feat
randcoeff[['b_brandhiland']] <- brandhiland + d_brandhiland * sd_brandhiland
randcoeff[['b_brandweight']] <- brandweight + d_brandweight * sd_brandweight
randcoeff[['b_brandyoplait']] <- brandyoplait + d_brandyoplait * sd_brandyoplait
return(randcoeff)
}

apollo_control_wtp <- list(
modelName = "MXL_WTP_space",
modelDescr = "MXL model on yogurt choice SP data, in WTP space",
indivID = "id",
mixing = TRUE,
analyticGrad = TRUE,
panelData = TRUE,
nCores = 1
)

apollo_draws_n <- list(
interDrawsType = "halton",
interNDraws = numDraws_wtp,
interUnifDraws = c(),
interNormDraws = c(
"d_feat", "d_brandhiland", "d_brandweight", "d_brandyoplait"),
intraDrawsType = "halton",
intraNDraws = 0,
intraUnifDraws = c(),
intraNormDraws = c()
)

apollo_inputs_wtp <- apollo_validateInputs(
apollo_beta = start_wtp,
apollo_fixed = NULL,
database = data_apollo,
apollo_draws = apollo_draws_n,
apollo_randCoeff = apollo_randCoeff,
apollo_control = apollo_control_wtp
)

data_mixl <- data_apollo
data_mixl$ID <- data_mixl$id
data_mixl$CHOICE <- data_mixl$choice

mixl_wtp <- "
feat_RND = @feat + draw_1 * @sd_feat;
brandhiland_RND = @brandhiland + draw_2 * @sd_brandhiland;
brandweight_RND = @brandweight + draw_3 * @sd_brandweight;
brandyoplait_RND = @brandyoplait + draw_4 * @sd_brandyoplait;
U_1 = @scalePar * (feat_RND * $feat_dannon - $price_dannon);
U_2 = @scalePar * (brandhiland_RND + feat_RND * $feat_hiland - $price_hiland);
U_3 = @scalePar * (brandweight_RND + feat_RND * $feat_weight - $price_weight);
U_4 = @scalePar * (brandyoplait_RND + feat_RND * $feat_yoplait - $price_yoplait);
"
mixl_spec_wtp <- specify_model(mixl_wtp, data_mixl)
availabilities <- generate_default_availabilities(data_mixl, 4)

model_mixl <- estimate(
mixl_spec_wtp, start_wtp,
data_mixl, availabilities,
nDraws = numDraws_wtp
)

LL_apollo_start <- apollo_llCalc(start_wtp, apollo_probabilities_wtp, apollo_inputs_wtp, silent = TRUE)$model[1]
LL_mixl_start <- sum(model_mixl$initLL)

`

[jhelvy]

Just noting this is on my radar to address, but I won't probably for another two weeks until the semester is over