Q1a

Create the data.frame that keeps all the observations whose value of container is either ‘CAN’ or ‘NON REFILLABLE BOTTLE’ in the given data.frame, beer_markets.

table(beer_markets$container)

## 
##                            CAN                            KEG                       KEG BALL 
##                          53015                            118                              2 
##          NON REFILLABLE BOTTLE NON REFILLABLE BOTTLE ALUMINUM  NON REFILLABLE BOTTLE PLASTIC 
##                          19095                            186                            543 
##              REFILLABLE BOTTLE 
##                            156

beer_markets <- filter(beer_markets, 
                       container == 'CAN' | 
                       container =='NON REFILLABLE BOTTLE')

Q1b

Create the data.frame that has factor-type variables of container and market using the resulting data.frame from Q1a.

beer_markets <- beer_markets %>% 
  mutate(container = factor(container),
         market = factor(market))

levels(beer_markets$container)

## [1] "CAN"                   "NON REFILLABLE BOTTLE"

levels(beer_markets$market)

##  [1] "ALBANY"               "ATLANTA"              "BALTIMORE"           
##  [4] "BIRMINGHAM"           "BOSTON"               "BUFFALO-ROCHESTER"   
##  [7] "CHARLOTTE"            "CHICAGO"              "CINCINNATI"          
## [10] "CLEVELAND"            "COLUMBUS"             "DALLAS"              
## [13] "DENVER"               "DES MOINES"           "DETROIT"             
## [16] "EXURBAN NY"           "GRAND RAPIDS"         "HARTFORD-NEW HAVEN"  
## [19] "HOUSTON"              "INDIANAPOLIS"         "JACKSONVILLE"        
## [22] "KANSAS CITY"          "LITTLE ROCK"          "LOS ANGELES"         
## [25] "LOUISVILLE"           "MEMPHIS"              "MIAMI"               
## [28] "MILWAUKEE"            "MINNEAPOLIS"          "NASHVILLE"           
## [31] "NEW ORLEANS-MOBILE"   "OKLAHOMA CITY-TULSA"  "OMAHA"               
## [34] "ORLANDO"              "PHILADELPHIA"         "PHOENIX"             
## [37] "PITTSBURGH"           "PORTLAND, OR"         "RALEIGH-DURHAM"      
## [40] "RICHMOND"             "RURAL ALABAMA"        "RURAL ARKANSAS"      
## [43] "RURAL CALIFORNIA"     "RURAL COLORADO"       "RURAL FLORIDA"       
## [46] "RURAL GEORGIA"        "RURAL IDAHO"          "RURAL ILLINOIS"      
## [49] "RURAL INDIANA"        "RURAL IOWA"           "RURAL KANSAS"        
## [52] "RURAL KENTUCKY"       "RURAL LOUISIANA"      "RURAL MAINE"         
## [55] "RURAL MICHIGAN"       "RURAL MINNESOTA"      "RURAL MISSISSIPPI"   
## [58] "RURAL MISSOURI"       "RURAL MONTANA"        "RURAL NEBRASKA"      
## [61] "RURAL NEVADA"         "RURAL NEW HAMPSHIRE"  "RURAL NEW MEXICO"    
## [64] "RURAL NEW YORK"       "RURAL NORTH CAROLINA" "RURAL NORTH DAKOTA"  
## [67] "RURAL OHIO"           "RURAL OKLAHOMA"       "RURAL OREGON"        
## [70] "RURAL PENNSYLVANIA"   "RURAL SOUTH CAROLINA" "RURAL SOUTH DAKOTA"  
## [73] "RURAL TENNESSEE"      "RURAL TEXAS"          "RURAL VERMONT"       
## [76] "RURAL VIRGINIA"       "RURAL WASHINGTON"     "RURAL WEST VIRGINIA" 
## [79] "RURAL WISCONSIN"      "RURAL WYOMING"        "SACRAMENTO"          
## [82] "SALT LAKE CITY"       "SAN ANTONIO"          "SAN DIEGO"           
## [85] "SAN FRANCISCO"        "SEATTLE"              "ST. LOUIS"           
## [88] "SURBURBAN NY"         "SYRACUSE"             "TAMPA"               
## [91] "URBAN NY"             "WASHINGTON DC"

table(beer_markets$container)

## 
##                   CAN NON REFILLABLE BOTTLE 
##                 53015                 19095

table(beer_markets$market)

## 
##               ALBANY              ATLANTA            BALTIMORE           BIRMINGHAM 
##                  487                 1279                  374                 1137 
##               BOSTON    BUFFALO-ROCHESTER            CHARLOTTE              CHICAGO 
##                  872                  607                 1246                 1879 
##           CINCINNATI            CLEVELAND             COLUMBUS               DALLAS 
##                 1270                 1226                 1862                 2098 
##               DENVER           DES MOINES              DETROIT           EXURBAN NY 
##                  796                  716                 1731                  321 
##         GRAND RAPIDS   HARTFORD-NEW HAVEN              HOUSTON         INDIANAPOLIS 
##                  739                  370                 1673                 1213 
##         JACKSONVILLE          KANSAS CITY          LITTLE ROCK          LOS ANGELES 
##                  501                  663                  452                 1564 
##           LOUISVILLE              MEMPHIS                MIAMI            MILWAUKEE 
##                  833                  530                 2616                  728 
##          MINNEAPOLIS            NASHVILLE   NEW ORLEANS-MOBILE  OKLAHOMA CITY-TULSA 
##                  801                  989                  852                  800 
##                OMAHA              ORLANDO         PHILADELPHIA              PHOENIX 
##                 1017                 1135                  433                 2263 
##           PITTSBURGH         PORTLAND, OR       RALEIGH-DURHAM             RICHMOND 
##                  352                  552                 1126                 1063 
##        RURAL ALABAMA       RURAL ARKANSAS     RURAL CALIFORNIA       RURAL COLORADO 
##                  305                  160                  848                   21 
##        RURAL FLORIDA        RURAL GEORGIA          RURAL IDAHO       RURAL ILLINOIS 
##                  522                  460                  154                 1195 
##        RURAL INDIANA           RURAL IOWA         RURAL KANSAS       RURAL KENTUCKY 
##                  481                 1060                  179                  225 
##      RURAL LOUISIANA          RURAL MAINE       RURAL MICHIGAN      RURAL MINNESOTA 
##                  381                  353                  754                  138 
##    RURAL MISSISSIPPI       RURAL MISSOURI        RURAL MONTANA       RURAL NEBRASKA 
##                  354                  640                  354                  110 
##         RURAL NEVADA  RURAL NEW HAMPSHIRE     RURAL NEW MEXICO       RURAL NEW YORK 
##                  557                   25                  427                   13 
## RURAL NORTH CAROLINA   RURAL NORTH DAKOTA           RURAL OHIO       RURAL OKLAHOMA 
##                  909                  129                  257                   54 
##         RURAL OREGON   RURAL PENNSYLVANIA RURAL SOUTH CAROLINA   RURAL SOUTH DAKOTA 
##                   38                  298                 1295                  153 
##      RURAL TENNESSEE          RURAL TEXAS        RURAL VERMONT       RURAL VIRGINIA 
##                  423                 1771                  139                  185 
##     RURAL WASHINGTON  RURAL WEST VIRGINIA      RURAL WISCONSIN        RURAL WYOMING 
##                  330                  265                 1306                   39 
##           SACRAMENTO       SALT LAKE CITY          SAN ANTONIO            SAN DIEGO 
##                  981                  320                 2615                  656 
##        SAN FRANCISCO              SEATTLE            ST. LOUIS         SURBURBAN NY 
##                  871                  903                 1347                  872 
##             SYRACUSE                TAMPA             URBAN NY        WASHINGTON DC 
##                  294                 3180                  735                  863

Q1c

Create the data.frame that has a factor-type variable of market whose reference level is “BUFFALO-ROCHESTER” using the resulting data.frame from Q1b.

beer_markets$market <- relevel(beer_markets$market, 
                               "BUFFALO-ROCHESTER")

Q1d

Split the resulting data.frame of Q1c into training and testing data.frames such that approximately 67% of observations in the resulting data.frame of Q1c belong to the training data.frame and the rest observations belong to the testing data.frame.

set.seed(1234)
gp <- runif(nrow(beer_markets))
train <- filter(beer_markets, gp > .33)  # training data
test <- filter(beer_markets, gp <= .33)  # test data

Q1e

Consider the three linear regression models with the following three different formulas:

formula_1 <- log(price_per_floz) ~ market + container + brand + log(beer_floz)

formula_2 <- log(price_per_floz) ~ market + container + brand * log(beer_floz)

formula_3 <- log(price_per_floz) ~ market + container + brand * promo * log(beer_floz)

Conduct exploratory data analysis regarding the model with formula_3.

sum_train <- skimr::skim(
  select(train,
         price_per_floz, beer_floz, 
         container, brand, promo, market)
)

library(GGally)
ggpairs( select(train,
                price_per_floz, beer_floz, 
                container, brand, promo)  )

- More visualization should be considered.

Q1f

• Train the three linear regression models using the training data.frame from Q1d with the three different formulas, formula_1, formula_2, and formula_3, provided in Q1e.

• What are the main goals of these linear regression models? Provide the summary of the result for each linear regression model.

• Create the data.frame of all the beta estimates for all variables and their t-statistics and p-values from each model.

model_1 <- lm(formula_1, data = train)
model_2 <- lm(formula_2, data = train)
model_3 <- lm(formula_3, data = train)

# summary(model_1)
# summary(model_2)
# summary(model_3)

library(stargazer)
stargazer(model_1, model_2, model_3, type = 'html')

stargazer(model_1, model_2, model_3, type = 'html', 
          omit = c("market", "container", "Constant"))

b1 <- coef(model_1)
b2 <- coef(model_2)
b3 <- coef(model_3)

library(broom)
df_b1 <- tidy(model_1)
df_b2 <- tidy(model_2)
df_b3 <- tidy(model_3)

• The goal of these linear regression is to answer the following questions:
• 1. How are the explanatory variables—sales volume, promotion, and brand characteristics (e.g, consumer loyalty, beer taste)—related with the beer price?
• 2. Given the linear relationship between the outcome and the explanatory variables, what are the predicted beer prices for unseen data?
• Intuition behind the model:
  • I choose the pricing model. The beer industry is close to oligopoly, which means that each beer company has enough market power to set the price level above its competitive level of price.
  • Each beer company sets the price of a brand’s beer based on the degree of consumers’ loyalty, sales volume, promotion, container type, and market characteristics (as well as various consumer characteristics).

Q1g

Interpret the beta estimates of the following variables from the model with formula_3. - (1) marketALBANY - (2) marketEXURBAN NY - (3) marketRURAL NEW YORK - (4) marketSURBURBAN NY - (5) marketSYRACUSE - (6) marketURBAN NY

q1g_df_b3 <- df_b3 %>% 
  filter(str_detect(term, "marketALBANY") |
           str_detect(term, "marketEXURBAN NY") |
           str_detect(term, "marketRURAL NEW YORK") |
           str_detect(term, "marketSURBURBAN NY") |
           str_detect(term, "marketSYRACUSE") |
           str_detect(term, "marketURBAN NY") ) %>% 
  mutate(sinificance = ifelse(p.value <= .1, T, F)) %>% 
  select(sinificance, everything())
  # beta for marketRURAL NEW YORK is not statistically significant

b3['marketALBANY']  # ***

## marketALBANY 
##    0.0358602

exp(b3['marketALBANY'])  # 1.036511 

## marketALBANY 
##     1.036511

• 1. All else being equal, being in ALBANY relative to being in BUFFALO-ROCHESTER is associated with an increase in price_per_floz by 3.65%.

b3['marketEXURBAN NY']  # ***

## marketEXURBAN NY 
##        0.1808736

exp(b3['marketEXURBAN NY'])

## marketEXURBAN NY 
##         1.198264

• 2. All else being equal, being in EXURBAN NY relative to being in BUFFALO-ROCHESTER is associated with an increase in price_per_floz by 19.83%.

b3['marketRURAL NEW YORK']  # there is no star

## marketRURAL NEW YORK 
##          -0.02045906

• 3. We cannot reject the null hypothesis that the beta coefficient for the variable marketRURAL NEW YORK is zero.

b3['marketSURBURBAN NY'] # ***

## marketSURBURBAN NY 
##         0.05639214

exp(b3['marketSURBURBAN NY'])

## marketSURBURBAN NY 
##           1.058012

• 4. All else being equal, being in SURBURBAN NY relative to being in BUFFALO-ROCHESTER is associated with an increase in price_per_floz by 5.80%

b3['marketSYRACUSE'] # ***

## marketSYRACUSE 
##    -0.04721861

exp(b3['marketSYRACUSE'])

## marketSYRACUSE 
##      0.9538788

1 - exp(b3['marketSYRACUSE'])

## marketSYRACUSE 
##     0.04612116

• 5. All else being equal, being in SYRACUSE relative to being in BUFFALO-ROCHESTER is associated with a decrease in price_per_floz by 4.61%.

b3['marketURBAN NY'] # ***

## marketURBAN NY 
##      0.1652982

exp(b3['marketURBAN NY'])

## marketURBAN NY 
##       1.179745

• 6. All else being equal, being in URBAN NY relative to being in BUFFALO-ROCHESTER is associated with an increase in price_per_floz by 17.97%.

Q1h

• Across the three models in Q1f, how is the percentage change in the price of beer sensitive to the percentage change in the volume of beer purchases for each brand?

• In the model with formula_3, how does such sensitivity vary by promo?

q1h_df_b1 <- tidy(model_1) %>% 
  filter( !str_detect(term, "market"),
          !str_detect(term, "container") ) %>% 
  mutate(sinificance = ifelse(p.value <= .1, T, F)) %>% 
  select(sinificance, everything())

q1h_df_b2 <- tidy(model_2) %>% 
  filter( !str_detect(term, "market"),
          !str_detect(term, "container") ) %>% 
  mutate(sinificance = ifelse(p.value <= .1, T, F)) %>% 
  select(sinificance, everything()) 

q1h_df_b3 <- tidy(model_3) %>% 
  filter( !str_detect(term, "market"),
          !str_detect(term, "container")  ) %>% 
  mutate(sinificance = ifelse(p.value <= .1, T, F)) %>% 
  select(sinificance, everything())

# Model 1
q1h_df_b1

(sensitivity1_ALL <- b1['log(beer_floz)'])

## log(beer_floz) 
##     -0.1408958

• All else being equal, 1% increase in beer_floz is associated with a decrease in price_per_floz by 0.141%.

# Model 2
q1h_df_b2

## BUD LIGHT
(sensitivity2_BUD <- b2['log(beer_floz)'])

## log(beer_floz) 
##     -0.1457584

• All else being equal, an increase in BUD LIGHT beer_floz by 1% is associated with a decrease in its price_per_floz by 0.146%.

q1h_df_b2

## BUSCH LIGHT
(sensitivity2_BUSCH <- 
    b2['log(beer_floz)'] + b2['brandBUSCH LIGHT:log(beer_floz)'])

## log(beer_floz) 
##     -0.1533481

• All else being equal, an increase in BUSCH LIGHT beer_floz by 1% is associated with a decrease in its price_per_floz by 0.153%.

q1h_df_b2

## COORS LIGHT 
(sensitivity2_COORS <- b2['log(beer_floz)'])

## log(beer_floz) 
##     -0.1457584

• All else being equal, an increase in COORS LIGHT beer_floz by 1% is associated with a decrease in its price_per_floz by 0.146%.

• cf) There are no stars (or dots) next to b2['brandCOORS LIGHT:log(beer_floz)'], meaning that it is not statistically significant.

• So we cannot reject the null hypotheses that the beta coefficient for the variable, brandCOORS LIGHT:log(beer_floz) is zero.

q1h_df_b2

## MILLER LITE
(sensitivity2_MILLER <- 
    b2['log(beer_floz)'] + b2['brandMILLER LITE:log(beer_floz)'])

## log(beer_floz) 
##     -0.1673114

• All else being equal, an increase in MILLER LITE beer_floz by 1% is associated with a decrease in its price_per_floz by 0.167%.

q1h_df_b2

## NATURAL LIGHT
(sensitivity2_NATURAL <- 
    b2['log(beer_floz)'] + b2['brandNATURAL LIGHT:log(beer_floz)'])

## log(beer_floz) 
##    -0.09202679

• All else being equal, an increase in NATURAL LIGHT beer_floz by 1% is associated with a decrease in its price_per_floz by 0.092%.

# Model 3 
q1h_df_b3

## BUD LIGHT
(sensitivity3_BUD_NOpromo <- b3['log(beer_floz)'])

## log(beer_floz) 
##     -0.1410682

(sensitivity3_BUD_promo <- b3['log(beer_floz)'])

## log(beer_floz) 
##     -0.1410682

• All else being equal, an increase in BUD LIGHT sales by 1% with or without promo is associated with a decrease in its price_per_floz by 0.141%

• cf) b3['promoTRUE:log(beer_floz)'] is not statistically significant.

q1h_df_b3

## BUSCH LIGHT
(sensitivity3_BUSCH_NOpromo <- 
    b3['log(beer_floz)'] + b3['brandBUSCH LIGHT:log(beer_floz)'])

## log(beer_floz) 
##     -0.1561043

(sensitivity3_BUSCH_promo <- 
    b3['log(beer_floz)'] + b3['brandBUSCH LIGHT:log(beer_floz)'] + 
  b3['brandBUSCH LIGHT:promoTRUE:log(beer_floz)'])

## log(beer_floz) 
##     -0.1152271

• All else being equal, an increase in BUSCH LIGHT sales by 1% without promo is associated with a decrease in its price_per_floz by 0.156%.

• All else being equal, an increase in BUSCH LIGHT sales by 1% with promo is associated with a decrease in its price_per_floz by 0.115%.

q1h_df_b3

## COORS LIGHT 
(sensitivity3_COORS_NOpromo <- b3['log(beer_floz)'])

## log(beer_floz) 
##     -0.1410682

(sensitivity3_COORS_promo <- 
    b3['log(beer_floz)']+ b3['brandCOORS LIGHT:promoTRUE:log(beer_floz)'])

## log(beer_floz) 
##     -0.1173336

• All else being equal, an increase in COORS LIGHT sales by 1% without promo is associated with a decrease in its price_per_floz by 0.141%.

• All else being equal, an increase in COORS LIGHT sales by 1% with promo is associated with a decrease in its price_per_floz by 0.117%.

q1h_df_b3

## MILLER LITE
(sensitivity3_MILLER_NOpromo <- 
    b3['log(beer_floz)'] + b3['brandMILLER LITE:log(beer_floz)'])

## log(beer_floz) 
##     -0.1654951

(sensitivity3_MILLER_promo <- 
    b3['log(beer_floz)'] + b3['brandMILLER LITE:log(beer_floz)'] + 
  b3['brandMILLER LITE:promoTRUE:log(beer_floz)'])

## log(beer_floz) 
##     -0.1352884

• All else being equal, an increase in MILLER LITE sales by 1% without promo is associated with a decrease in its price_per_floz by 0.165%.

• All else being equal, an increase in MILLER LITE sales by 1% with promo is associated with a decrease in its price_per_floz by 0.135%.

q1h_df_b3

## NATURAL LIGHT
(sensitivity3_NATURAL_NOpromo <- 
    b3['log(beer_floz)'] + b3['brandNATURAL LIGHT:log(beer_floz)'])

## log(beer_floz) 
##    -0.09924425

(sensitivity3_NATURAL_promo <- 
    b3['log(beer_floz)'] + b3['brandNATURAL LIGHT:log(beer_floz)'] + 
  b3['brandNATURAL LIGHT:promoTRUE:log(beer_floz)'])

## log(beer_floz) 
##    -0.04637121

• All else being equal, an increase in NATURAL LIGHT sales by 1% without promo is associated with a decrease in its price_per_floz by 0.099%.

• All else being equal, an increase in NATURAL LIGHT sales by 1% with promo is associated with a decrease in its price_per_floz by 0.046%.

# all elasticities in one data.frame:
sensitivity <- data.frame(
  brand = c('BUD LIGHT', 'BUSCH LIGHT', 'COORS LIGHT', 'MILLER LITE', 'NATURAL LIGHT'),
  model1 = round(sensitivity1_ALL, digits = 4),
  model2 = round(c(sensitivity2_BUD, sensitivity2_BUSCH, sensitivity2_COORS, sensitivity2_MILLER, sensitivity2_NATURAL), digits = 4),
  model3_NOpromo = round(c(sensitivity3_BUD_NOpromo, sensitivity3_BUSCH_NOpromo, sensitivity3_COORS_NOpromo, sensitivity3_MILLER_NOpromo, sensitivity3_NATURAL_NOpromo), digits = 4),
  model3_promo = round(c(sensitivity3_BUD_promo, sensitivity3_BUSCH_promo, sensitivity3_COORS_promo, sensitivity3_MILLER_promo, sensitivity3_NATURAL_promo), digits = 4)
)

sensitivity

• In model_1, the volume sensitivity to changes in price does not vary by brand.

• In model_2, the volume sensitivity to changes in price varies by brand.

  • MILLER LITE is the most sensitive one.
  • NATURAL LIGHT is the least sensitive one.
  • BUD LIGHT and COORS LIGHT have the same sensitivity.

• In model_3, the brand-specific volume sensitivity to changes in price can vary by promo.

  • For all beer brands except for BUD LIGHT, promo == TRUE made their demands less sensitive!

Q1i

• For each model in Q1f, draw a residual plot.
• On average, are the prediction correct? Are there systematic errors?
• Which model do you prefer most and why?

test <- test %>% 
  mutate( pred_1 = predict(model_1, newdata = test),
          pred_2 = predict(model_2, newdata = test),
          pred_3 = predict(model_3, newdata = test) )


# residual plot
p1 <- ggplot(data = test, aes(x = pred_1, 
                        y = log(price_per_floz) - pred_1 )) +
  geom_point(color = 'grey60', alpha = .1) +
  geom_smooth(color = 'blue') +
  geom_hline(aes(yintercept = 0), color = 'red') +  # y = 0
  theme_minimal()
p1

p2 <- ggplot(data = test, aes(x = pred_2, 
                        y = log(price_per_floz) - pred_2 )) +
  geom_point(color = 'grey60', alpha = .1) +
  geom_smooth(color = 'blue') +
  geom_hline(aes(yintercept = 0), color = 'red') +  # y = 0
  theme_minimal()
p2

p3 <- ggplot(data = test, aes(x = pred_3, 
                        y = log(price_per_floz) - pred_3 )) +
  geom_point(color = 'grey60', alpha = .1) +
  geom_smooth(color = 'blue') +
  geom_hline(aes(yintercept = 0), color = 'red') +  # y = 0
  theme_minimal()
p3

- When answering the question—Which model do you prefer most and why?—, we can consider intuition behind formulas across models and the significance of beta estimates.

• Which model do you prefer most and why?
  • I prefer model 3.
  • The beta estimates for promo-related interaction terms except for BUD LIGHT are statistically significant.
  • Intuitively, promotion affects purchasing decision of consumers who are sufficiently sensitive to changes in beer volumes.
  • So, promo is likely to be a relevant variable to explain the variation of the beer price.