inst/doc/convergence.R

In logitr: Logit Models w/Preference & WTP Space Utility Parameterizations

## ----include=FALSE------------------------------------------------------------
# Unfortunately, because not all these packages are available on CRAN 
# for all platforms, this vignette causes errors that will lead to {logitr}
# being taken off CRAN. So instead of actually running all of this code, 
# I just hard-code the results that print out here. You can run it yourself
# to verify that the results are accurate.

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  warning = FALSE,
  message = FALSE,
  eval = FALSE,
  fig.width = 7.252,
  fig.height = 4,
  comment = "#>",
  fig.retina = 3
)

# Read in results from already estimated models so that the
# examples aren't actually run when building this page,
# otherwise it'll take much longer to build
# model_logitr <- readRDS(here::here('inst', 'extdata', 'model_logitr.Rds'))
# model_logitr10 <- readRDS(here::here('inst', 'extdata', 'model_logitr10.Rds'))
# model_mixl1 <- readRDS(here::here('inst', 'extdata', 'model_mixl1.Rds'))
# model_mixl2 <- readRDS(here::here('inst', 'extdata', 'model_mixl2.Rds'))
# model_gmnl1 <- readRDS(here::here('inst', 'extdata', 'model_gmnl1.Rds'))
# model_apollo1 <- readRDS(here::here('inst', 'extdata', 'model_apollo1.Rds'))
# model_apollo2 <- readRDS(here::here('inst', 'extdata', 'model_apollo2.Rds'))

## -----------------------------------------------------------------------------
# library(logitr)
# library(mlogit)
# library(gmnl)
# library(apollo)
# library(mixl)
# library(dplyr)
# library(tidyr)
# 
# set.seed(1234)

## -----------------------------------------------------------------------------
# numDraws_wtp <- 50
# start_wtp <- c(
#     scalePar        = 1,
#     feat            = 0,
#     brandhiland     = 0,
#     brandweight     = 0,
#     brandyoplait    = 0,
#     sd_feat         = 0.1,
#     sd_brandhiland  = 0.1,
#     sd_brandweight  = 0.1,
#     sd_brandyoplait = 0.1
# )

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

## -----------------------------------------------------------------------------
# data_gmnl <- mlogit.data(
#     data     = yogurt,
#     shape    = "long",
#     choice   = "choice",
#     id.var   = 'id',
#     alt.var  = 'alt',
#     chid.var = 'obsID'
# )

## -----------------------------------------------------------------------------
# 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"
# ) {
# 
#     ### 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[["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)
# 
#     ### Define settings for MNL model component
#     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
#     )
# 
#     ### 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)
# }

## -----------------------------------------------------------------------------
# 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
# )
# 
# # Set parameters for generating draws
# 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()
# )
# 
# # Set input
# 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 # Uses the same "wide" format as {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_logitr <- logitr(
#     data      = yogurt,
#     outcome   = 'choice',
#     obsID     = 'obsID',
#     panelID   = 'id',
#     pars      = c('feat', 'brand'),
#     scalePar  = 'price',
#     randPars  = c(feat = "n", brand = "n"),
#     startVals = start_wtp,
#     numDraws  = numDraws_wtp
# )
# 
# summary(model_logitr)

## -----------------------------------------------------------------------------
# model_logitr <- logitr(
#     data      = yogurt,
#     outcome   = 'choice',
#     obsID     = 'obsID',
#     panelID   = 'id',
#     pars      = c('feat', 'brand'),
#     scalePar  = 'price',
#     randPars  = c(feat = "n", brand = "n"),
#     startVals = start_wtp,
#     numDraws  = numDraws_wtp,
#     numMultiStarts = 10
# )
# 
# summary(model_logitr)

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

## -----------------------------------------------------------------------------
# c(logLik(model_logitr), logLik(model_mixl))

## -----------------------------------------------------------------------------
# cbind(coef(model_logitr), coef(model_mixl))

## -----------------------------------------------------------------------------
# model_mixl <- estimate(
#     mixl_spec_wtp, coef(model_logitr),
#     data_mixl, availabilities,
#     nDraws = numDraws_wtp
# )

## -----------------------------------------------------------------------------
# c(logLik(model_logitr), logLik(model_mixl))

## -----------------------------------------------------------------------------
# cbind(coef(model_logitr), coef(model_mixl))

## -----------------------------------------------------------------------------
# model_gmnl <- gmnl(
#     data = data_gmnl,
#     formula = choice ~ price + feat + brand | 0 | 0 | 0 | 1,
#     ranp = c(
#         feat = "n", brandhiland = "n", brandweight = "n",
#         brandyoplait = "n"),
#     fixed = c(TRUE, rep(FALSE, 10), TRUE),
#     model = "gmnl",
#     method = "bfgs",
#     haltons = NA,
#     panel = TRUE,
#     start = c(start_wtp, 0.1, 0.1, 0),
#     R = numDraws_wtp
# )

## -----------------------------------------------------------------------------
# c(logLik(model_logitr), logLik(model_gmnl))

## -----------------------------------------------------------------------------
# cbind(coef(model_logitr), coef(model_gmnl))

## -----------------------------------------------------------------------------
# model_gmnl <- gmnl(
#     data = data_gmnl,
#     formula = choice ~ price + feat + brand | 0 | 0 | 0 | 1,
#     ranp = c(
#         feat = "n", brandhiland = "n", brandweight = "n",
#         brandyoplait = "n"),
#     fixed = c(TRUE, rep(FALSE, 10), TRUE),
#     model = "gmnl",
#     method = "bfgs",
#     haltons = NA,
#     panel = TRUE,
#     start = c(coef(model_logitr), 0.1, 0.1, 0),
#     R = numDraws_wtp
# )

## -----------------------------------------------------------------------------
# c(logLik(model_logitr), logLik(model_gmnl))

## -----------------------------------------------------------------------------
# cbind(coef(model_logitr), coef(model_gmnl))

## -----------------------------------------------------------------------------
# model_apollo <- apollo_estimate(
#     apollo_beta          = start_wtp,
#     apollo_fixed         = NULL,
#     apollo_probabilities = apollo_probabilities_wtp,
#     apollo_inputs        = apollo_inputs_wtp,
#     estimate_settings    = list(printLevel = 0)
# )

## -----------------------------------------------------------------------------
# c(logLik(model_logitr), model_apollo$LLout)

## -----------------------------------------------------------------------------
# cbind(coef(model_logitr), coef(model_apollo))

## -----------------------------------------------------------------------------
# model_apollo <- apollo_estimate(
#     apollo_beta          = coef(model_logitr),
#     apollo_fixed         = NULL,
#     apollo_probabilities = apollo_probabilities_wtp,
#     apollo_inputs        = apollo_inputs_wtp,
#     estimate_settings    = list(printLevel = 0)
# )

## -----------------------------------------------------------------------------
# c(logLik(model_logitr), model_apollo$LLout)

## -----------------------------------------------------------------------------
# cbind(coef(model_logitr), model_apollo$betaStop)