inst/doc/v04_modeling_heterogeneity.R

In RprobitB: Bayesian Probit Choice Modeling

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "img/",
  fig.align = "center",
  fig.dim = c(8, 6), 
  out.width = "75%"
)
library("RprobitB")
options("RprobitB_progress" = FALSE)

## ---- message = FALSE---------------------------------------------------------
data("Electricity", package = "mlogit")
Electricity <- as_cov_names(Electricity, c("pf","cl","loc","wk","tod","seas"), 1:4)
data <- prepare_data(
  form = choice ~ pf + cl + loc + wk + tod + seas | 0,
  choice_data = Electricity,
  re = c("cl","loc","wk","tod","seas")
)
model_elec <- fit_model(data, R = 1000, scale = "pf := -1")

## ---- coef-model-elec---------------------------------------------------------
coef(model_elec)

## ---- compute-mixdistr-share--------------------------------------------------
cl_mu <- coef(model_elec)["cl","mean"]
cl_sd <- sqrt(coef(model_elec)["cl","var"])
pnorm(cl_mu / cl_sd)

## ---- cov-mixdistr------------------------------------------------------------
cov_mix(model_elec, cor = TRUE)

## ---- plot-mixture-model-elec, fig.dim = c(8,8)-------------------------------
plot(model_elec, type = "mixture")

## ---- sim-dirichlet-----------------------------------------------------------
set.seed(1)
P_r <- 2
C_true <- 3
N <- c(100,70,30)
(b_true <- matrix(replicate(C_true, rnorm(P_r)), nrow = P_r, ncol = C_true))
(Omega_true <- matrix(replicate(C_true, rwishart(P_r + 1, 0.1*diag(P_r))$W, simplify = TRUE), 
                      nrow = P_r*P_r, ncol = C_true))
beta <- c()
for(c in 1:C_true) for(n in 1:N[c])
  beta <- cbind(beta, rmvnorm(mu = b_true[,c,drop=F], Sigma = matrix(Omega_true[,c,drop=F], ncol = P_r)))
z_true <- rep(1:3, times = N)

## ---- dirichlet-prior---------------------------------------------------------
delta <- 0.1
xi <- numeric(P_r)
D <- diag(P_r)
nu <- P_r + 2
Theta <- diag(P_r)

## ---- dirichlet-inits---------------------------------------------------------
z <- rep(1, ncol(beta))
C <- 1
b <- matrix(0, nrow = P_r, ncol = C)
Omega <- matrix(rep(diag(P_r), C), nrow = P_r*P_r, ncol = C)

## ---- dirichlet-process-app---------------------------------------------------
for(r in 1:100){
  dp <- RprobitB:::update_classes_dp(
    Cmax = 10, beta, z, b, Omega, delta, xi, D, nu, Theta, s_desc = TRUE
    )
  z <- dp$z
  b <- dp$b
  Omega <- dp$Omega
}

## ---- dirichlet-example-plot, fig.dim = c(16,8)-------------------------------
par(mfrow = c(1,2))
plot(t(beta), xlab = bquote(beta[1]), ylab = bquote(beta[2]), pch = 19)
RprobitB:::plot_class_allocation(beta, z, b, Omega, r = 100, perc = 0.95)