Solutions: Exercises 7–8

library(tidyverse)
library(broom)

hotel <- read_csv("data/hotel_upgrades.csv", show_col_types = FALSE) |>
  mutate(
    loyalty_status = factor(loyalty_status, levels = c("None", "Silver", "Gold", "Platinum")),
    special_event = factor(special_event),
    direct_booking = factor(direct_booking)
  )

mean(hotel$upgrade)

[1] 0.2411111

upgrade_model <- glm(
  upgrade ~ loyalty_status + prior_stays + occupancy_rate +
    total_spend_eur + special_event + direct_booking,
  family = binomial,
  data = hotel
)
tidy(upgrade_model, exponentiate = TRUE, conf.int = TRUE)

# A tibble: 9 × 7
  term                   estimate std.error statistic p.value conf.low conf.high
  <chr>                     <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 (Intercept)               0.955  0.572      -0.0806 9.36e-1    0.310     2.93 
2 loyalty_statusSilver      1.66   0.237       2.14   3.26e-2    1.04      2.63 
3 loyalty_statusGold        1.90   0.320       2.01   4.44e-2    1.01      3.54 
4 loyalty_statusPlatinum    4.44   0.490       3.04   2.35e-3    1.70     11.7  
5 prior_stays               1.22   0.0557      3.51   4.43e-4    1.09      1.36 
6 occupancy_rate            0.966  0.00699    -4.99   6.08e-7    0.952     0.979
7 total_spend_eur           1.00   0.000425    1.50   1.34e-1    1.000     1.00 
8 special_eventYes          0.558  0.257      -2.27   2.32e-2    0.330     0.909
9 direct_bookingYes         1.67   0.194       2.64   8.18e-3    1.15      2.46 

profiles <- tibble(
  loyalty_status = factor(c("None", "Gold"), levels = levels(hotel$loyalty_status)),
  prior_stays = 3,
  occupancy_rate = 75,
  total_spend_eur = 500,
  special_event = factor("No", levels = levels(hotel$special_event)),
  direct_booking = factor("Yes", levels = levels(hotel$direct_booking))
)

profiles$probability <- predict(
  upgrade_model,
  newdata = profiles,
  type = "response"
)

profiles

# A tibble: 2 × 7
  loyalty_status prior_stays occupancy_rate total_spend_eur special_event
  <fct>                <dbl>          <dbl>           <dbl> <fct>        
1 None                     3             75             500 No           
2 Gold                     3             75             500 No           
# ℹ 2 more variables: direct_booking <fct>, probability <dbl>

prepare <- function(path) {
  read_csv(path, show_col_types = FALSE) |>
    mutate(
      engagement_4 = 6 - engagement_4_reverse,
      engagement = rowMeans(across(c(engagement_1, engagement_2, engagement_3, engagement_4)), na.rm = TRUE)
    )
}

original <- prepare("data/employee_survey.csv")
replication <- prepare("data/employee_replication.csv")

original_model <- lm(performance ~ engagement + workload + training_group, data = original)
replication_model <- lm(performance ~ engagement + workload + training_group, data = replication)

comparison <- bind_rows(
  tidy(original_model, conf.int = TRUE) |> mutate(sample = "Original"),
  tidy(replication_model, conf.int = TRUE) |> mutate(sample = "Replication")
)
comparison

# A tibble: 8 × 8
  term           estimate std.error statistic  p.value conf.low conf.high sample
  <chr>             <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl> <chr> 
1 (Intercept)       43.9      1.74      25.3  1.91e-89   40.5      47.3   Origi…
2 engagement         4.86     0.282     17.3  4.43e-52    4.31      5.42  Origi…
3 workload          -1.54     0.312     -4.94 1.07e- 6   -2.16     -0.930 Origi…
4 training_grou…     2.45     0.622      3.93 9.60e- 5    1.22      3.67  Origi…
5 (Intercept)       45.4      2.23      20.3  3.45e-61   41.0      49.8   Repli…
6 engagement         4.11     0.357     11.5  2.70e-26    3.41      4.82  Repli…
7 workload          -1.43     0.405     -3.54 4.57e- 4   -2.23     -0.636 Repli…
8 training_grou…     2.46     0.786      3.13 1.88e- 3    0.917     4.01  Repli…

write_csv(comparison, "replication_model_comparison.csv")