Let’s say our utility function is:
\[v_j = \beta_1 x_j^{\mathrm{price}} + \beta_2 x_j^{\mathrm{cacao}} + \beta_3 \delta_j^{\mathrm{hershey}} + \beta_4 \delta_j^{\mathrm{lindt}}\]
And we estimate the following coefficients:
Parameter | Coefficient |
---|---|
\(\beta_1\) | -0.1 |
\(\beta_2\) | 0.1 |
\(\beta_3\) | -2.0 |
\(\beta_4\) | -0.1 |
What are the expected probabilities of choosing each of these bars using a logit model?
Attribute | Bar 1 | Bar 2 | Bar 3 |
---|---|---|---|
Price | $1.20 | $1.50 | $3.00 |
% Cacao | 10% | 60% | 80% |
Brand | Hershey | Lindt | Ghirardelli |
# Create the model coefficient variables
b1 <- -0.1
b2 <- 0.1
b3 <- -2.0
b4 <- -0.1
# Compute the observed utility for each chocolate bar
v1 <- b1*1.2 + b2*0.1 + b3
v2 <- b1*1.5 + b2*0.6 + b4
v3 <- b1*3 + b2*0.8
# Compute the probabilities using the logit fraction
denom <- exp(v1) + exp(v2) + exp(v3)
p1 <- exp(v1) / denom
p2 <- exp(v2) / denom
p3 <- exp(v3) / denom
p1
#> [1] 0.06925051
p2
#> [1] 0.4723548
p3
#> [1] 0.4583946
Observations - Height of students (inches):
#> [1] 65 69 66 67 68 72 68 69 63 70
# Load the data
x <- c(65, 69, 66, 67, 68, 72, 68, 69, 63, 70)
# Compute the value of f(x) for each x
f_x <- dnorm(x, 68, 4)
# The likelihood is just the product of the probabilities in f_x
prod(f_x)
#> [1] 1.447528e-11
# But this is a really tiny number, so computing the log-likelihood is helpful
log(prod(f_x))
#> [1] -24.95858
# Of course, the way we typically compute the log-likelihood is by summing up
# the log of the values in f_x
sum(log(f_x))
#> [1] -24.95858
library(tidyverse)
# Create a vectors of values for the mean
means <- c(66, 67, 68, 69, 70)
# Compute the likelihood using different values for the mean:
L1 <- sum(log(dnorm(x, means[1], 4)))
L2 <- sum(log(dnorm(x, means[2], 4)))
L3 <- sum(log(dnorm(x, means[3], 4)))
L4 <- sum(log(dnorm(x, means[4], 4)))
L5 <- sum(log(dnorm(x, means[5], 4)))
logLiks <- c(L1, L2, L3, L4, L5)
# Create a data frame of the results
df <- data.frame(means, logLiks)
df
#> means logLiks
#> 1 66 -25.83358
#> 2 67 -25.08358
#> 3 68 -24.95858
#> 4 69 -25.45858
#> 5 70 -26.58358
# Filter out the row with the maximum likelihood value:
df %>%
filter(logLiks == max(logLiks))
#> means logLiks
#> 1 68 -24.95858
# Plot the result:
df %>%
ggplot(aes(x = means, y = logLiks)) +
geom_line() +
geom_point() +
theme_bw() +
labs(x = "Mean Value", y = "Log-likelihood Values")
Suppose we estimate the following utility model describing preferences for cars:
\[ u_j = \alpha p_j + \beta_1 x_j^{mpg} + \beta_2 x_j^{elec} + \varepsilon_j \]
The estimated model produces the following results:
Parameter | Coefficient |
---|---|
\(\alpha\) | -0.7 |
\(\beta_1\) | 0.1 |
\(\beta_2\) | -0.4 |
Hessian:
\[ \begin{bmatrix} -6000 & 50 & 60 \\ 50 & -700 & 50 \\ 60 & 50 & -300 \end{bmatrix} \]
Compute a 95% confidence interval around the coefficients using:
# 1. Get coefficients
beta <- c(price = -0.7, mpg = 0.1, elec = -4.0)
# 2. Get covariance matrix
hessian <- matrix(c(
-6000, 50, 60,
50, -700, 50,
60, 50, -300),
ncol = 3, byrow = TRUE)
covariance <- -1*solve(hessian)
se <- sqrt(diag(covariance))
lower <- beta - 2*se
upper <- beta + 2*se
coef_ci <- data.frame(
mean = beta,
lower = lower,
upper = upper
)
coef_ci
#> mean lower upper
#> price -0.7 -0.72585699 -0.674143
#> mpg 0.1 0.02392002 0.176080
#> elec -4.0 -4.11629585 -3.883704
draws <- as.data.frame(MASS::mvrnorm(10^4, beta, covariance))
coef_ci <- maddTools::ci(draws, ci = 0.95)
coef_ci
#> mean lower upper
#> price -0.7001425 -0.72497169 -0.675090
#> mpg 0.1004593 0.02665909 0.174476
#> elec -3.9998885 -4.11361564 -3.883845