Suppose we estimate the following utility model describing preferences for cars:
\[ \tilde{u}_j = \beta_1 x_j^{\mathrm{price}} + \beta_2 x_j^{\mathrm{mpg}} + \beta_3 x_j^{\mathrm{elec}} + \varepsilon_j \]
Solution:
\[ \tilde{u}_j = \beta_1 x_j^{\mathrm{price}} + \beta_2 x_j^{\mathrm{mpg}} + \beta_3 x_j^{\mathrm{elec}} + \beta_4 x_j^{\mathrm{price}}\delta^\mathrm{B} + \beta_5 x_j^{\mathrm{mpg}}\delta^\mathrm{B} + \beta_6 x_j^{\mathrm{elec}}\delta^\mathrm{B} + \varepsilon_j \]
Solution:
Effect | Group A | Group B |
---|---|---|
\(x_j^\mathrm{price}\) | \(\beta_1\) | \(\beta_1 + \beta_4\) |
\(x_j^\mathrm{mpg}\) | \(\beta_2\) | \(\beta_2 + \beta_5\) |
\(\delta_j^\mathrm{elec}\) | \(\beta_3\) | \(\beta_3 + \beta_6\) |
Suppose we estimate the following utility model describing preferences for chocolate bars between two groups: A & B
\[ \tilde{u}_j = \beta_1 x_j^{\mathrm{price}} + \beta_2 x_j^{\mathrm{caco}} + \beta_3 x_j^{\mathrm{price}}\delta_j^{\mathrm{B}} + \beta_4 x_j^{\mathrm{cacao}}\delta_j^{\mathrm{B}} + \varepsilon_j \]
The estimated model produces the following coefficients and hessian:
\(\beta\) = [-0.7, 0.1, 0.2, 0.8]
\[ H = \begin{bmatrix} -6000 & 50 & 60 & 70 \\ 50 & -700 & 50 & 100 \\ 60 & 50 & -300 & 20 \\ 70 & 100 & 20 & -6000 \end{bmatrix} \]
mvrnorm()
function from the MASS
library to generate 10,000 draws of the model coefficients.beta <- c(b1 = -0.7, b2 = 0.1, b3 = 0.2, b4 = 0.8)
hessian <- matrix(c(
-6000, 50, 60, 70,
50, -700, 50, 100,
60, 50, -300, 20,
70, 100, 20, -6000),
ncol = 4, byrow = TRUE)
covariance <- -1*(solve(hessian))
draws <- as.data.frame(MASS::mvrnorm(10^5, beta, covariance))
# First, separate draws for each group
draws_A <- draws %>%
select(
price = b1,
cacao = b2)
draws_B <- draws %>%
mutate(
price = b1 + b3,
cacao = b2 + b4) %>%
select(price, cacao)
head(draws_A)
#> price cacao
#> 1 -0.7103629 0.098974739
#> 2 -0.6745787 0.140467828
#> 3 -0.7080966 0.143054548
#> 4 -0.7221457 0.110694519
#> 5 -0.6648105 -0.001851373
#> 6 -0.6820227 0.153294795
head(draws_B)
#> price cacao
#> 1 -0.6121547 0.9114146
#> 2 -0.5288541 0.9348620
#> 3 -0.5550391 0.9629691
#> 4 -0.5400009 0.9266010
#> 5 -0.5059122 0.8043792
#> 6 -0.5446152 0.9636818
# Compute WTP
wtp_A <- draws_A / (-1*draws_A$price)
wtp_B <- draws_B / (-1*draws_B$price)
head(wtp_A)
#> price cacao
#> 1 -1 0.139329835
#> 2 -1 0.208230451
#> 3 -1 0.202026877
#> 4 -1 0.153285579
#> 5 -1 -0.002784813
#> 6 -1 0.224764941
head(wtp_B)
#> price cacao
#> 1 -1 1.488863
#> 2 -1 1.767712
#> 3 -1 1.734957
#> 4 -1 1.715925
#> 5 -1 1.589958
#> 6 -1 1.769473
# Get a confidence interval
ci(wtp_A)
#> mean lower upper
#> price -1.0000000 -1.00000000 -1.000000
#> cacao 0.1429716 0.03711834 0.250247
ci(wtp_B)
#> mean lower upper
#> price -1.000000 -1.000000 -1.000000
#> cacao 1.828346 1.412967 2.403132
logitr
package to estimate the following homogeneous model:\[ \tilde{u}_j = \beta_1 x_j^{\mathrm{price}} + \beta_2 \delta_j^{\mathrm{feat}} + \beta_3 \delta_j^{\mathrm{dannon}} + \beta_4 \delta_j^{\mathrm{hiland}} + \beta_5 \delta_j^{\mathrm{weight}} + \varepsilon_j \]
where the three \(\delta\) coefficients are dummy variables for Dannon, Hiland, and Weight Watchers brands (Yoplait is the reference level).
# First, need to make some dummy coefficients for "brand"
yogurt <- fastDummies::dummy_cols(yogurt, "brand")
model <- logitr(
data = yogurt,
obsID = "obsID",
outcome = "choice",
pars = c("price", "feat", "brand_dannon", "brand_hiland", "brand_weight")
)
summary(model)
#> =================================================
#> Call:
#> logitr(data = yogurt, outcome = "choice", obsID = "obsID", pars = c("price",
#> "feat", "brand_dannon", "brand_hiland", "brand_weight"))
#>
#> Frequencies of alternatives:
#> 1 2 3 4
#> 0.402156 0.029436 0.229270 0.339138
#>
#> Exit Status: 3, Optimization stopped because ftol_rel or ftol_abs was reached.
#>
#> Model Type: Multinomial Logit
#> Model Space: Preference
#> Model Run: 1 of 1
#> Iterations: 19
#> Elapsed Time: 0h:0m:0.1s
#> Algorithm: NLOPT_LD_LBFGS
#> Weights Used?: FALSE
#> Robust? FALSE
#>
#> Model Coefficients:
#> Estimate Std. Error z-value Pr(>|z|)
#> price -0.366508 0.024364 -15.0431 < 2.2e-16 ***
#> feat 0.491744 0.120058 4.0959 4.206e-05 ***
#> brand_dannon -0.734420 0.080639 -9.1075 < 2.2e-16 ***
#> brand_hiland -4.449491 0.187097 -23.7818 < 2.2e-16 ***
#> brand_weight -1.375509 0.088975 -15.4594 < 2.2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Log-Likelihood: -2656.8878885
#> Null Log-Likelihood: -3343.7419990
#> AIC: 5323.7757771
#> BIC: 5352.7168000
#> McFadden R2: 0.2054148
#> Adj McFadden R2: 0.2039195
#> Number of Observations: 2412.0000000
logitr
package to estimate the same model but with the following mixing distributions:model_mxl <- logitr(
data = yogurt,
obsID = "obsID",
outcome = "choice",
pars = c("price", "feat", "brand_dannon", "brand_hiland", "brand_weight"),
randPars = c(price = "n", feat = "n")
)
summary(model_mxl)
#> =================================================
#> Call:
#> logitr(data = yogurt, outcome = "choice", obsID = "obsID", pars = c("price",
#> "feat", "brand_dannon", "brand_hiland", "brand_weight"),
#> randPars = c(price = "n", feat = "n"))
#>
#> Frequencies of alternatives:
#> 1 2 3 4
#> 0.402156 0.029436 0.229270 0.339138
#>
#> Exit Status: 3, Optimization stopped because ftol_rel or ftol_abs was reached.
#>
#> Model Type: Mixed Logit
#> Model Space: Preference
#> Model Run: 1 of 1
#> Iterations: 27
#> Elapsed Time: 0h:0m:2s
#> Algorithm: NLOPT_LD_LBFGS
#> Weights Used?: FALSE
#> Robust? FALSE
#>
#> Model Coefficients:
#> Estimate Std. Error z-value Pr(>|z|)
#> price_mu -0.391475 0.026757 -14.6310 < 2.2e-16 ***
#> feat_mu 0.264899 0.207460 1.2769 0.2017
#> brand_dannon -0.794621 0.086362 -9.2011 < 2.2e-16 ***
#> brand_hiland -4.789186 0.216287 -22.1427 < 2.2e-16 ***
#> brand_weight -1.458237 0.095632 -15.2485 < 2.2e-16 ***
#> price_sigma 0.026399 0.091062 0.2899 0.7719
#> feat_sigma 2.433177 0.504727 4.8208 1.43e-06 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Log-Likelihood: -2645.8312988
#> Null Log-Likelihood: -3343.7419990
#> AIC: 5305.6625975
#> BIC: 5346.1801000
#> McFadden R2: 0.2087215
#> Adj McFadden R2: 0.2066280
#> Number of Observations: 2412.0000000
#>
#> Summary of 10k Draws for Random Coefficients:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> price (normal) -0.4928975 -0.409283 -0.3914787 -0.3915011 -0.3736765 -0.2946348
#> feat (normal) -9.1921630 -1.377908 0.2630395 0.2614922 1.9034268 9.0514876