3.1 Loading the Raw Data

I will ingest my data.

starbucks_raw <- read_csv(here("data", "starbucks-nutrition.csv")) %>%
  mutate(id = row_number()) %>%
  clean_names()

starbucks_raw

# A tibble: 1,147 × 16
   product_name             size   milk  whip serv_size_m_l calories total_fat_g
   <chr>                    <chr> <dbl> <dbl>         <dbl>    <dbl>       <dbl>
 1 brewed coffee - dark ro… short     0     0           236        3         0.1
 2 brewed coffee - dark ro… tall      0     0           354        4         0.1
 3 brewed coffee - dark ro… gran…     0     0           473        5         0.1
 4 brewed coffee - dark ro… venti     0     0           591        5         0.1
 5 brewed coffee - decaf p… short     0     0           236        3         0.1
 6 brewed coffee - decaf p… tall      0     0           354        4         0.1
 7 brewed coffee - decaf p… gran…     0     0           473        5         0.1
 8 brewed coffee - decaf p… venti     0     0           591        5         0.1
 9 brewed coffee - medium … short     0     0           236        3         0.1
10 brewed coffee - medium … tall      0     0           354        4         0.1
# … with 1,137 more rows, and 9 more variables: saturated_fat_g <dbl>,
#   trans_fat_g <chr>, cholesterol_mg <dbl>, sodium_mg <dbl>,
#   total_carbs_g <dbl>, fiber_g <chr>, sugar_g <dbl>, caffeine_mg <dbl>,
#   id <int>

3.2 Cleaning the Data

I will check the variables types and check for implausible values in my variables of interest.

starbucks_raw %>%
  select(product_name, size, milk, whip, calories, sugar_g, caffeine_mg, id) %>%
  summary()

 product_name           size                milk            whip       
 Length:1147        Length:1147        Min.   :0.000   Min.   :0.0000  
 Class :character   Class :character   1st Qu.:1.000   1st Qu.:0.0000  
 Mode  :character   Mode  :character   Median :2.000   Median :0.0000  
                                       Mean   :2.513   Mean   :0.2467  
                                       3rd Qu.:4.000   3rd Qu.:0.0000  
                                       Max.   :5.000   Max.   :1.0000  
    calories        sugar_g       caffeine_mg           id        
 Min.   :  0.0   Min.   : 0.00   Min.   :  0.00   Min.   :   1.0  
 1st Qu.:130.0   1st Qu.:18.00   1st Qu.: 30.00   1st Qu.: 287.5  
 Median :220.0   Median :34.00   Median : 75.00   Median : 574.0  
 Mean   :228.4   Mean   :34.99   Mean   : 91.86   Mean   : 574.0  
 3rd Qu.:320.0   3rd Qu.:49.00   3rd Qu.:150.00   3rd Qu.: 860.5  
 Max.   :640.0   Max.   :89.00   Max.   :475.00   Max.   :1147.0  

I need to change the variable type for some as well as recode and relevel. There are no implausible values.

starbucks <- starbucks_raw %>%
  select(id, product_name, size, milk, whip, calories, sugar_g, caffeine_mg) %>%
  mutate(id = as.character(id),
         size = factor(size),
         size = fct_relevel(size, "tall", 
                            "grande", 
                            "venti"),
         milk = factor(milk),
         milk = fct_recode(milk, "no_milk" = "0",
                           "nonfat" = "1",
                           "2_pct" = "2",
                           "soy" = "3",
                           "coconut" = "4",
                           "whole" = "5"),
         milk = fct_relevel(milk, "no_milk", 
                            "nonfat", 
                            "2_pct", 
                            "soy", 
                            "coconut", 
                            "whole"),
         whip = factor(whip),
         whip = fct_recode(whip, "whip" = "1",
                           "no_whip" = "0"),
         whip = fct_relevel(whip, "whip", 
                            "no_whip"),
         sug = sugar_g,
         caff = caffeine_mg,
         name = product_name,
         cal = calories) %>%
  select(-product_name, -calories, -sugar_g, -caffeine_mg) %>%
  droplevels()

summary(starbucks)

      id                 size          milk          whip          sug       
 Length:1147        grande :334   no_milk:165   whip   :283   Min.   : 0.00  
 Class :character   venti  :320   nonfat :222   no_whip:864   1st Qu.:18.00  
 Mode  :character   tall   :318   2_pct  :190                 Median :34.00  
                    short  :123   soy    :190                 Mean   :34.99  
                    trenta : 21   coconut:190                 3rd Qu.:49.00  
                    doppio :  7   whole  :190                 Max.   :89.00  
                    (Other): 24                                              
      caff            name                cal       
 Min.   :  0.00   Length:1147        Min.   :  0.0  
 1st Qu.: 30.00   Class :character   1st Qu.:130.0  
 Median : 75.00   Mode  :character   Median :220.0  
 Mean   : 91.86                      Mean   :228.4  
 3rd Qu.:150.00                      3rd Qu.:320.0  
 Max.   :475.00                      Max.   :640.0  
                                                    

This looks better. size, milk, and whip are factor variables and id and name are characters.

I need to change the categories for size.

summary(starbucks$size) %>%
  kable()

I will combine all the categories except tall, grande, and venti into a category called ‘other’ because the current categories are very small relative to the primary three.

starbucks <- starbucks %>%
  mutate(size = fct_recode(size, "other" = "short",
                           "other" = "1 scoop",
                           "other" = "1 shot",
                           "other" = "doppio",
                           "other" = "quad",
                           "other" = "solo",
                           "other" = "trenta",
                           "other" = "triple")) %>%
  droplevels()

summary(starbucks$size) %>%
  kable()

This looks better.

3.3 Missingness Check

I will check for missing values in my data set.

miss_var_summary(starbucks)

# A tibble: 8 × 3
  variable n_miss pct_miss
  <chr>     <int>    <dbl>
1 id            0        0
2 size          0        0
3 milk          0        0
4 whip          0        0
5 sug           0        0
6 caff          0        0
7 name          0        0
8 cal           0        0

I have no missing values.

4.1 Listing the Tibble

I will list my tibble.

starbucks

# A tibble: 1,147 × 8
   id    size   milk    whip      sug  caff name                             cal
   <chr> <fct>  <fct>   <fct>   <dbl> <dbl> <chr>                          <dbl>
 1 1     other  no_milk no_whip     0   130 brewed coffee - dark roast         3
 2 2     tall   no_milk no_whip     0   193 brewed coffee - dark roast         4
 3 3     grande no_milk no_whip     0   260 brewed coffee - dark roast         5
 4 4     venti  no_milk no_whip     0   340 brewed coffee - dark roast         5
 5 5     other  no_milk no_whip     0    15 brewed coffee - decaf pike pl…     3
 6 6     tall   no_milk no_whip     0    20 brewed coffee - decaf pike pl…     4
 7 7     grande no_milk no_whip     0    25 brewed coffee - decaf pike pl…     5
 8 8     venti  no_milk no_whip     0    30 brewed coffee - decaf pike pl…     5
 9 9     other  no_milk no_whip     0   155 brewed coffee - medium roast       3
10 10    tall   no_milk no_whip     0   235 brewed coffee - medium roast       4
# … with 1,137 more rows

4.2 Size and Identifiers

My tibble has 1147 rows and 8 columns. The name of the variable that provides the identifier for each row is id, which is a character variable. I will demonstrate that each row has a unique value for that variable here.

n_distinct(starbucks$id)

[1] 1147

4.3 Saving the R data set

I will save the data set as an .Rds file.

# saveRDS(starbucks, file = here("data", "starbucks"))

5.1 Defining the Variables

5.2 Numerical Description

Here is a numerical description of my tibble.

Hmisc::describe(starbucks)

starbucks 

 8  Variables      1147  Observations
--------------------------------------------------------------------------------
id 
       n  missing distinct 
    1147        0     1147 

lowest : 1    10   100  1000 1001, highest: 995  996  997  998  999 
--------------------------------------------------------------------------------
size 
       n  missing distinct 
    1147        0        4 
                                      
Value        tall grande  venti  other
Frequency     318    334    320    175
Proportion  0.277  0.291  0.279  0.153
--------------------------------------------------------------------------------
milk 
       n  missing distinct 
    1147        0        6 

lowest : no_milk nonfat  2_pct   soy     coconut
highest: nonfat  2_pct   soy     coconut whole  
                                                          
Value      no_milk  nonfat   2_pct     soy coconut   whole
Frequency      165     222     190     190     190     190
Proportion   0.144   0.194   0.166   0.166   0.166   0.166
--------------------------------------------------------------------------------
whip 
       n  missing distinct 
    1147        0        2 
                          
Value         whip no_whip
Frequency      283     864
Proportion   0.247   0.753
--------------------------------------------------------------------------------
sug 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
    1147        0       89    0.999    34.99     25.6        0        4 
     .25      .50      .75      .90      .95 
      18       34       49       67       76 

lowest :  0  1  2  3  4, highest: 85 86 87 88 89
--------------------------------------------------------------------------------
caff 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
    1147        0       60    0.994    91.86       84        0        0 
     .25      .50      .75      .90      .95 
      30       75      150      190      225 

lowest :   0  10  15  20  25, highest: 410 425 445 470 475
--------------------------------------------------------------------------------
name 
       n  missing distinct 
    1147        0       93 

lowest : Blended Strawberry Lemonade                   brewed coffee - dark roast                    brewed coffee - decaf pike place roast        brewed coffee - medium roast                  brewed coffee - True North Blend Blonde roast
highest: Vanilla Sweet Cream Cold Brew                 Very Berry Hibiscus Starbucks Refreshers      White Chocolate Mocha                         White Hot Chocolate                           Youthberry Brewed Tea                        
--------------------------------------------------------------------------------
cal 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
    1147        0       71    0.999    228.4    156.3        5       50 
     .25      .50      .75      .90      .95 
     130      220      320      420      470 

lowest :   0   2   3   4   5, highest: 590 600 620 630 640
--------------------------------------------------------------------------------

6.1 Question

The question I hope to answer with the proposed linear model is:

How effectively does milk type, amount of sugar, amount of caffeine, and beverage size predict the number of calories in a Starbucks beverage?

6.2 My Quantitative Outcome

My outcome is cal. I am interested in this because I enjoy Starbucks and tend to watch my caloric intake.

I will display a count of the number of rows in starbucks with complete information on cal.

starbucks %>%
  filter(complete.cases(.)) %>%
  select(cal) %>%
  nrow()

[1] 1147

Here is a histogram of cal.

ggplot(starbucks, aes(x = cal)) +
  geom_histogram(aes(y = stat(density)),
                 bins = 15,
                 col = "white",
                 fill = "darkgreen") +
  stat_function(fun = dnorm, 
                args = list(mean = mean(starbucks$cal), 
                            sd = sd(starbucks$cal)),
                col = "tan3", lwd = 1.5) +
  labs(title = "Distribution of `cal` with minor right skew",
       x = "# of calories",
       y = "subjects (density)",
       caption = "1147 beverages in `starbucks`")

While the distribution roughly approximates a normal distribution, there is a minor right skew. It is fairly continuous and I have displayed it in 15 bins. I considered using the natural log, but it caused a much more significant left skew than the right skew in this untransformed outcome.

Here I will show that there are at least 10, ordered, distinct values.

n_distinct(starbucks$cal)

[1] 71

mosaic::favstats(starbucks$cal) %>% 
  kable()

Registered S3 method overwritten by 'mosaic':
  method                           from   
  fortify.SpatialPolygonsDataFrame ggplot2

There are 71 distinct values, ranging from 0 to 640 calories, with a mean of 228 calories.

6.3 My Planned Predictors (Linear Model)

The predictors I will use are size, sug, caff, and milk.

I will demonstrate that I have an input which is quantitative, sug, with at least 10 different, ordered, observed values.

n_distinct(starbucks$sug)

[1] 89

mosaic::favstats(starbucks$sug) %>% 
  kable()

sug has 89 distinct values ranging from 0 to 89 grams, with a mean of 35 grams.

I will demonstrate that I have a categorical input which has between 3 and 6 categories, that will be used as a factor, and that has at least 30 observations in each level of the factor.

n_distinct(starbucks$milk)

[1] 6

tabyl(starbucks$milk) %>% 
  kable(digits = 3)

glimpse(starbucks$milk)

 Factor w/ 6 levels "no_milk","nonfat",..: 1 1 1 1 1 1 1 1 1 1 ...

milk has 6 categories with at least 30 observations in each. 2_pct, soy, coconut, and whole each have 190 beverages, nonfat contains 222 beverages, and no_milk contains 165 beverages. The variable is a factor.

The total number of candidate predictors is no more than 4 + (N-100)/100 candidate regression inputs, rounding down, where N is the sample size.

4 + (1147-100)/100 = 4 + 1047/100 = 4 + 10.47 = 14.47

7.1 Question

The question I hope to answer with the proposed logistic model is:

How effectively does beverage size, type of milk, and amount of sugar and caffeine predict whether a beverage has whipped cream?

7.2 My Binary Outcome

I will use my logistic regression model to explain whip. I am interested in examining this because I like whipped cream in my coffee sometimes and I am curious about how well these factors predict whether a beverage has whipped cream.

I will display a count of the number of rows with each of the two possible values of this outcome and a missingness check.

tabyl(starbucks$whip) %>%
  kable(digits = 3)

starbucks %>% select(whip) %>% 
  miss_var_summary() %>% 
  kable()

283 of the beverages, or 25%, have whipped cream and 864 beverages, or 75%, do not have whipped cream. I do not have any missing values.

7.3 My Planned Predictors (Logistic Model)

The predictors I will use in this model are size, sug, caff, and milk. In section 6.3 I demonstrated that I have a quantitative input with at least 10 different, ordered, observed values and a categorical input with between 3 and 6 categories that is a factor and has at least 30 observations in each level. The sample size and number of predictors are the same.

8.1 Missingness

My sample has no missing data.

8.2 Outcome Transformation

My quantitative outcome is fairly normally distributed. There is a minor right skew as shown in this histogram.

ggplot(starbucks, aes(x = cal)) +
  geom_histogram(aes(y = stat(density)),
                 bins = 15,
                 col = "white",
                 fill = "darkgreen") +
  stat_function(fun = dnorm, 
                args = list(mean = mean(starbucks$cal), 
                            sd = sd(starbucks$cal)),
                col = "tan3", lwd = 1.5) +
  labs(title = "Distribution of `cal` has minor right skew",
       x = "# of calories",
       y = "subjects (density)",
       caption = "1147 beverages from `starbucks`")

My outcome can be 0, so I will add 1 to every subject’s cal value in order to do my transformation assessment.

starbucks %$%
  boxCox((cal + 1) ~ sug + caff + milk + size)

My result is very close to 1, so it well approximates a normal distribution. I considered a logarithm transformation but it caused a left skew that was more substantial that my current minor right skew, so I will not transform my outcome.

8.3 Scatterplot Matrix and Collinearity

I will create a scatterplot matrix to evaluate the associations of my predictors and outcome.

ggpairs(starbucks %>%
          select(size, sug, caff, milk, cal))

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

It looks like sug is very closely correlated with cal which makes sense because sugar is known to increase calories of a beverage. It looks like caff does not have a strong relationship with cal which also makes sense, as caffeine does not contain calories. The relationship between the milk categories and cal show that no milk has much fewer calories than the rest of the categories. The size category ‘other’ has the fewest calories compared to the remaining categories, which makes sense because this category consists of smaller beverages such as an espresso shot.

I will look at each of my predictors’ ability to predict cal.

p_sug <- ggplot(starbucks, aes(x = sug,
                                      y = cal)) +
                  geom_point(alpha = 0.3,
                             col = "forestgreen") +
  labs(title = "Sugar",
       x = "sugar (g)",
       y = "# calories")
  
p_size <- ggplot(starbucks, aes(x = size, y = cal, 
                               fill = size)) +
    geom_bar(stat = "identity") +
  labs(title = "Size",
       x = "size",
       y = "# calories")
  
p_caff <- ggplot(starbucks, aes(x = caff,
                                      y = cal)) +
                  geom_point(alpha = 0.6,
                             col = "goldenrod") +
  labs(title = "Caffeine",
       x = "caffeine (mg)",
       y = "# calories")
  
p_milk <- ggplot(starbucks, aes(x = milk, y = cal,
                                fill = milk)) +
  geom_bar(stat = "identity") +
  labs(title = "Milk",
       x = "milk",
       y = "# calories")

(p_sug + p_size) / (p_caff + p_milk) +
  plot_annotation(title = "Sugar is Most Predictive for Number of Calories",
                  subtitle = "",
                  caption = "1147 beverages from `starbucks`")

This plot confirms what I saw in my scatterplot matrix, however in a more attractive format. Sugar in grams is highly correlated with calories. Caffeine does not have a linear association. The amount of calories increases with increasing size of the beverage. Milk type does not have an association, but the presence of any milk or no milk does look different.

I will look at variance inflation factors for collinearity.

vif(lm(cal ~ size + sug + caff + milk, data = starbucks)) %>% 
  kable()

All of the variance inflation factors are below 2, so I am not worried about collinearity between predictors.

8.4 Model A

modA_lm <- lm(cal ~ size + sug + caff + milk, 
              data = starbucks)

tidy(modA_lm, conf.int = TRUE, conf.level = 0.90) %>%
    select(term, estimate, conf.low, conf.high, std.error, p.value) %>%
    kable(digits = c(0,2,2,2,2,3))

glance(modA_lm) %>%
  select(r.squared, adj.r.squared, AIC, BIC) %>%
  kable(digits = c(4,4,0,0))

d <- datadist(starbucks)
options(datadist = "d")

modA_ols <- ols(cal ~ sug + size + caff + milk, 
                data = starbucks, x = TRUE, y = TRUE)

set.seed(2022)
validate(modA_ols)

          index.orig  training      test optimism index.corrected  n
R-square      0.8578    0.8585    0.8566   0.0019          0.8559 40
MSE        2693.0025 2655.2618 2715.8641 -60.6023       2753.6048 40
g           145.1467  144.4762  145.0659  -0.5896        145.7363 40
Intercept     0.0000    0.0000   -0.0310   0.0310         -0.0310 40
Slope         1.0000    1.0000    1.0015  -0.0015          1.0015 40

par(mfrow=c(2,2))
plot(modA_lm)

8.5 Non-Linearity

I will assess whether a nonlinear term might improve my model by starting with a Spearman rho-squared plot.

sp <- spearman2(cal ~ size + sug + caff + milk, 
                data = starbucks)

plot(sp)

The Spearman rho squared plot suggests the strongest predictor is sug. Since this is a quantitative variable, I will spend my degrees of freedom on trying a restricted cubic spline with 5 knots, which will add 3 degrees of freedom. I want to use about 6 degrees of freedom beyond the main effects to account for non-linearity as per the project requirements, so I will also create an interaction term between size and caff which will add an additional 3 degrees of freedom. I will use size in the interaction term because it has the second largest rho-squared. I do not want to create an interaction with milk as this will add far too many degrees of freedom (15 df). caff has a low rho-squared, but I need a variable with 1 degree of freedom and I am already using sug in a nonlinear term, so I will use this variable in my interaction term.

8.6 Model B

modB_ols <- ols(cal ~ rcs(sug, 5) + size * caff + milk, 
                data = starbucks, x = TRUE, y = TRUE)

summary(modB_ols) %>% 
  kable(digits = c(2))

plot(summary(modB_ols))

plot(nomogram(modB_ols))

anova(modB_ols)

                Analysis of Variance          Response: cal 

 Factor                                     d.f. Partial SS   MS         
 sug                                           4  7409886.678 1852471.669
  Nonlinear                                    3     8838.133    2946.044
 size  (Factor+Higher Order Factors)           6   105837.593   17639.599
  All Interactions                             3    11710.417    3903.472
 caff  (Factor+Higher Order Factors)           4    29977.029    7494.257
  All Interactions                             3    11710.417    3903.472
 milk                                          5   765638.136  153127.627
 size * caff  (Factor+Higher Order Factors)    3    11710.417    3903.472
 TOTAL NONLINEAR + INTERACTION                 6    22581.446    3763.574
 REGRESSION                                   16 18654155.290 1165884.706
 ERROR                                      1130  3066292.377    2713.533
 F      P     
 682.68 <.0001
   1.09 0.3541
   6.50 <.0001
   1.44 0.2299
   2.76 0.0265
   1.44 0.2299
  56.43 <.0001
   1.44 0.2299
   1.39 0.2165
 429.66 <.0001
              

modB_lm <- lm(cal ~ rcs(sug, 5) + size * caff + milk, data = starbucks)

par(mfrow=c(2,2))
plot(modB_lm)

starbucks %>%
  filter(id == "80") %>%
  select(id, name, size, milk, sug, cal) %>%
  kable()

8.7 Validating the Models

I will start by partitioning my data in training and test sets.

set.seed(2022)

split_samples <- starbucks$cal %>%
  createDataPartition(p = 0.7, list = FALSE)

starb_train <- starbucks[split_samples,]
starb_test <- starbucks[-split_samples,]

dim(starb_train)

[1] 804   8

dim(starb_test)

[1] 343   8

804 beverages, or 70%, are in the training subset and 343 beverages, or 30%, are in the test subset. My 8 columns remained in each subset.

modA_lm_train <- 
  lm(cal ~ size + sug + caff + milk, 
     data = starb_train)

tidy(modA_lm_train) %>%
  kable(digits = 3)

modA_pred <- modA_lm_train %>%
  predict(starb_test)

modA_pred[1:3]

        1         2         3 
 13.61833 -13.41668 102.90898 

postResample(modA_pred, starb_test$cal)

      RMSE   Rsquared        MAE 
55.8912927  0.8312869 44.2853299 

set.seed(2022)
validate(modB_ols, method = "boot", B = 40)

          index.orig  training      test optimism index.corrected  n
R-square      0.8588    0.8601    0.8569   0.0032          0.8556 40
MSE        2673.3151 2624.7217 2709.9485 -85.2268       2758.5419 40
g           145.4735  144.7860  145.2518  -0.4659        145.9394 40
Intercept     0.0000    0.0000    0.2788  -0.2788          0.2788 40
Slope         1.0000    1.0000    1.0004  -0.0004          1.0004 40

starb_linear_summaries <- tibble(
  model = c("main_effects", "augmented"),
  RMSE = c(RMSE(modA_pred, starb_test$cal), sqrt(2755.8552)),
  R2 = c(R2(modA_pred, starb_test$cal), 0.8559))

starb_linear_summaries %>%
  kable(digits = 3)

8.8 Final Model

I prefer the main effects model A. This is a list of the point estimates and confidence intervals for the model fit in the entire data set.

summary(modA_ols, conf.type=c('simultaneous')) %>% 
  kable(digits = 2)

This is the estimates shown graphically in a plot of effect size of the main effects model in the entire data set.

plot(summary(modA_ols))

The effect of sug on my outcome cal in my model is meaningful. If I have 2 beverages, Beverage A and Beverage B, that are the same size with the same amount of caffeine and same type of milk, but Beverage A has 18 grams of sugar and Beverage B has 49 grams of sugar, then my main effects model predicts that Beverage B’s number of calories will be 151.07 calories higher than Beverage A’s.

This is the validated R-squared for the main effects model.

set.seed(2022)
validate(modA_ols, method = "boot", B = 40)

          index.orig  training      test optimism index.corrected  n
R-square      0.8578    0.8585    0.8566   0.0019          0.8559 40
MSE        2693.0025 2655.2618 2715.8641 -60.6023       2753.6048 40
g           145.1467  144.4762  145.0659  -0.5896        145.7363 40
Intercept     0.0000    0.0000   -0.0310   0.0310         -0.0310 40
Slope         1.0000    1.0000    1.0015  -0.0015          1.0015 40

The estimated R-squared across n = 40 training samples was 0.8585, and was similar in the resulting test samples at 0.8566, suggesting optimism of 0.0019. The augmented model’s R-squared, corrected for overfitting, is 0.8559.

This is a nomogram for the main effects model.

plot(nomogram(modA_ols))

For a tall beverage with 50 grams of sugar, 50 mg of caffeine, and whole milk added, the total points added up on the nomogram would be 50 (sug) + 1 (size) + 1 (caff) + 22 (milk) = 79 total points which correlates with a prediction of about 325 calories.

I will calculate this more exactly and obtain a prediction interval.

new_sub <- tibble(size = "tall", sug = 50, 
                  caff = 50, milk = "whole")

predict(modA_lm, newdata = new_sub, 
        interval = "prediction", level = 0.95)

     fit      lwr      upr
1 326.86 224.1175 429.6025

The predicted calories for this beverage is 326.86 calories with a 95% uncertainty interval of 224.12 - 429.60 calories.

9.1 Missingness

I have no missing data.

starbucks %>%
  tabyl(whip, milk) %>%
  kable()

starbucks %>%
  tabyl(whip, size) %>%
  kable()

9.2 Model Y

modY_lrm <- lrm((whip == "whip") ~ size + sug + caff + milk, 
               data = starbucks, x = TRUE, y = TRUE)

modY_lrm

Logistic Regression Model
 
 lrm(formula = (whip == "whip") ~ size + sug + caff + milk, data = starbucks, 
     x = TRUE, y = TRUE)
 
                        Model Likelihood    Discrimination    Rank Discrim.    
                              Ratio Test           Indexes          Indexes    
 Obs          1147    LR chi2     258.19    R2       0.300    C       0.801    
  FALSE        864    d.f.            10    g        1.782    Dxy     0.602    
  TRUE         283    Pr(> chi2) <0.0001    gr       5.944    gamma   0.603    
 max |deriv| 2e-08                          gp       0.222    tau-a   0.224    
                                            Brier    0.147                     
 
              Coef    S.E.   Wald Z Pr(>|Z|)
 Intercept    -4.5105 0.4986 -9.05  <0.0001 
 size=grande  -0.8679 0.2172 -4.00  <0.0001 
 size=venti   -1.8521 0.2872 -6.45  <0.0001 
 size=other    0.3195 0.2760  1.16  0.2470  
 sug           0.0677 0.0058 11.67  <0.0001 
 caff         -0.0012 0.0013 -0.88  0.3769  
 milk=nonfat   1.3971 0.4520  3.09  0.0020  
 milk=2_pct    1.5123 0.4530  3.34  0.0008  
 milk=soy      1.6532 0.4566  3.62  0.0003  
 milk=coconut  1.6869 0.4542  3.71  0.0002  
 milk=whole    1.5123 0.4530  3.34  0.0008  
 

modY_glm <- glm((whip == "whip") ~ size + sug + caff + milk, 
               data = starbucks, family = "binomial"(link = "logit"))

tidy(modY_glm, conf.int = TRUE, conf.level = 0.90) %>%
    select(term, estimate, conf.low, conf.high, std.error, p.value) %>%
    kable(digits = c(0,2,2,2,2,3))

plot(summary(modY_lrm))

plot(nomogram(modY_lrm))

modY_aug <- augment(modY_glm, starbucks, type.predict = "response")

modY_aug %$%
  confusionMatrix(
    data = factor(.fitted >= 0.5),
    reference = factor(whip == "whip"),
    positive = "TRUE")

Confusion Matrix and Statistics

          Reference
Prediction FALSE TRUE
     FALSE   791  188
     TRUE     73   95
                                          
               Accuracy : 0.7724          
                 95% CI : (0.7471, 0.7964)
    No Information Rate : 0.7533          
    P-Value [Acc > NIR] : 0.06948         
                                          
                  Kappa : 0.291           
                                          
 Mcnemar's Test P-Value : 1.708e-12       
                                          
            Sensitivity : 0.33569         
            Specificity : 0.91551         
         Pos Pred Value : 0.56548         
         Neg Pred Value : 0.80797         
             Prevalence : 0.24673         
         Detection Rate : 0.08282         
   Detection Prevalence : 0.14647         
      Balanced Accuracy : 0.62560         
                                          
       'Positive' Class : TRUE            
                                          

9.3 Non-Linearity

I will assess whether a nonlinear term might improve my model by starting with a Spearman rho-squared plot.

sp <- spearman2((whip == "whip") ~ size + sug + caff + milk,
               data = starbucks)

plot(sp)

The Spearman rho squared plot suggests the strongest predictor is sug. Since this is a quantitative variable, I will spend my degrees of freedom on trying a restricted cubic spline with 5 knots, which will add 3 degrees of freedom.

9.4 Model Z

modZ_lrm <- lrm((whip == "whip") ~ rcs(sug, 5) + size + caff + milk, 
               data = starbucks, x = TRUE, y = TRUE)

modZ_lrm

Logistic Regression Model
 
 lrm(formula = (whip == "whip") ~ rcs(sug, 5) + size + caff + 
     milk, data = starbucks, x = TRUE, y = TRUE)
 
                        Model Likelihood    Discrimination    Rank Discrim.    
                              Ratio Test           Indexes          Indexes    
 Obs          1147    LR chi2     271.76    R2       0.314    C       0.806    
  FALSE        864    d.f.            13    g        1.925    Dxy     0.611    
  TRUE         283    Pr(> chi2) <0.0001    gr       6.856    gamma   0.611    
 max |deriv| 4e-05                          gp       0.226    tau-a   0.227    
                                            Brier    0.147                     
 
              Coef    S.E.   Wald Z Pr(>|Z|)
 Intercept    -4.6971 0.7879 -5.96  <0.0001 
 sug          -0.0732 0.0783 -0.93  0.3499  
 sug'          2.9307 1.2449  2.35  0.0186  
 sug''        -4.3537 1.8239 -2.39  0.0170  
 sug'''        2.0294 0.8774  2.31  0.0207  
 size=grande  -0.8489 0.2246 -3.78  0.0002  
 size=venti   -1.6380 0.2860 -5.73  <0.0001 
 size=other    0.6803 0.3170  2.15  0.0319  
 caff         -0.0008 0.0013 -0.60  0.5516  
 milk=nonfat   1.5510 0.4682  3.31  0.0009  
 milk=2_pct    1.6585 0.4687  3.54  0.0004  
 milk=soy      1.8123 0.4736  3.83  0.0001  
 milk=coconut  1.8403 0.4709  3.91  <0.0001 
 milk=whole    1.6585 0.4687  3.54  0.0004  
 

modZ_glm <- glm((whip == "whip") ~ rcs(sug, 5) + size + caff + milk, 
               data = starbucks, family = "binomial"(link = "logit"))

tidy(modZ_glm, conf.int = TRUE, conf.level = 0.90) %>%
    select(term, estimate, conf.low, conf.high, std.error, p.value) %>%
    kable(digits = c(0,2,2,2,2,3))

plot(summary(modZ_lrm))

plot(nomogram(modZ_lrm))

anova(modZ_lrm)

                Wald Statistics          Response: (whip == "whip") 

 Factor     Chi-Square d.f. P     
 sug        134.10      4   <.0001
  Nonlinear  12.17      3   0.0068
 size        44.63      3   <.0001
 caff         0.35      1   0.5516
 milk        16.56      5   0.0054
 TOTAL      152.63     13   <.0001

modZ_aug <- augment(modZ_glm, starbucks, type.predict = "response")

modZ_aug %$%
  confusionMatrix(
    data = factor(.fitted >= 0.5),
    reference = factor(whip == "whip"),
    positive = "TRUE"  )

Confusion Matrix and Statistics

          Reference
Prediction FALSE TRUE
     FALSE   793  191
     TRUE     71   92
                                          
               Accuracy : 0.7716          
                 95% CI : (0.7462, 0.7956)
    No Information Rate : 0.7533          
    P-Value [Acc > NIR] : 0.07927         
                                          
                  Kappa : 0.2833          
                                          
 Mcnemar's Test P-Value : 1.955e-13       
                                          
            Sensitivity : 0.32509         
            Specificity : 0.91782         
         Pos Pred Value : 0.56442         
         Neg Pred Value : 0.80589         
             Prevalence : 0.24673         
         Detection Rate : 0.08021         
   Detection Prevalence : 0.14211         
      Balanced Accuracy : 0.62146         
                                          
       'Positive' Class : TRUE            
                                          

9.5 Validating the Models

I will validate the Nagelkerke R-square and the C-statistic for the “augmented model” Z and compare them to that of the “main effects” model Y through the validate function.

set.seed(2022)
validate(modY_lrm, method = "boot", B = 40)

          index.orig training    test optimism index.corrected  n
Dxy           0.6024   0.6095  0.5931   0.0164          0.5859 40
R2            0.2996   0.3066  0.2905   0.0161          0.2835 40
Intercept     0.0000   0.0000 -0.0305   0.0305         -0.0305 40
Slope         1.0000   1.0000  0.9580   0.0420          0.9580 40
Emax          0.0000   0.0000  0.0147   0.0147          0.0147 40
D             0.2242   0.2300  0.2166   0.0134          0.2108 40
U            -0.0017  -0.0017  0.0003  -0.0020          0.0003 40
Q             0.2260   0.2317  0.2163   0.0154          0.2106 40
B             0.1472   0.1453  0.1489  -0.0036          0.1508 40
g             1.7824   1.8144  1.7278   0.0866          1.6958 40
gp            0.2225   0.2243  0.2189   0.0054          0.2171 40

validate(modZ_lrm, method = "boot", B = 40)

          index.orig training    test optimism index.corrected  n
Dxy           0.6112   0.6309  0.5994   0.0315          0.5798 40
R2            0.3135   0.3355  0.3021   0.0334          0.2801 40
Intercept     0.0000   0.0000 -0.0748   0.0748         -0.0748 40
Slope         1.0000   1.0000  0.9097   0.0903          0.9097 40
Emax          0.0000   0.0000  0.0342   0.0342          0.0342 40
D             0.2361   0.2556  0.2264   0.0292          0.2068 40
U            -0.0017  -0.0017  0.0013  -0.0031          0.0013 40
Q             0.2378   0.2573  0.2250   0.0323          0.2055 40
B             0.1467   0.1435  0.1488  -0.0053          0.1519 40
g             1.9251   2.0877  1.8704   0.2173          1.7078 40
gp            0.2263   0.2347  0.2220   0.0127          0.2136 40

starb_logistic_summaries <- tibble(
  model = c("main_effects", "augmented"),
  R2 = c(0.2835, 0.2801),
  C = c((0.5 + (0.5859/2)), (0.5 + (0.5798/2))),
  )

starb_logistic_summaries %>%
  kable(digits = 4)

9.6 Final Model

I prefer the main effects model Y. Here is a list of the odds ratios and confidence intervals for the model fit in the entire data set. Since a binary logistic regression model predicts the log odds of whipped cream (whip == “whip”), I will exponentiate to interpret the odds ratios and confidence intervals.

tidy(modY_glm, exponentiate = TRUE, conf.int = TRUE) %>% 
  kable(digits = 3)

I can also run a summary of the lrm model to see a more detailed summary.

summary(modY_lrm) %>%
  kable(digits = 3)

This is a graphical display of the the effect sizes of the main effects model in the entire data set.

plot(summary(modY_lrm))

The effect of sug on my outcome whip in my model is meaningful, as shown by the odds ratio interval being entirely greater than 1. If I have 2 beverages, Beverage A and Beverage B, that are the same size with the same amount of caffeine and same type of milk, but Beverage A has 18 grams of sugar and Beverage B has 49 grams of sugar, then my augmented model predicts that Beverage B’s odds of having whipped cream will be 8.164 times Beverage A’s odds of having whipped cream. Beverage B’s odds are 816.4% as large as Beverage A’s.

This is a plot of the ROC curve for the main effects model in the entire data set.

prob <- predict(modY_glm, starbucks, type="response")
pred <- prediction(prob, starbucks$whip)
perf <- performance(pred, measure = "tpr", x.measure = "fpr")
auc <- performance(pred, measure="auc")

auc <- round(auc@y.values[[1]],3)
roc.data <- data.frame(fpr=unlist(perf@x.values),
                       tpr=unlist(perf@y.values),
                       model="LRM")

ggplot(roc.data, aes(x=fpr, ymin=0, ymax=tpr)) +
    geom_ribbon(alpha=0.2, fill = "blue") +
    geom_line(aes(y=tpr), col = "blue") +
    geom_abline(intercept = 0, slope = 1, lty = "dashed") +
    labs(title = paste0("Model Y: ROC Curve w/ AUC=", auc))

This plot displays that my model predicts the odds of a beverage having whipped cream better than guessing.

This is the validated Nagelkerke R-square and C-statistic for the main effects model.

set.seed(2022)
validate(modY_lrm, method = "boot", B = 40)

          index.orig training    test optimism index.corrected  n
Dxy           0.6024   0.6095  0.5931   0.0164          0.5859 40
R2            0.2996   0.3066  0.2905   0.0161          0.2835 40
Intercept     0.0000   0.0000 -0.0305   0.0305         -0.0305 40
Slope         1.0000   1.0000  0.9580   0.0420          0.9580 40
Emax          0.0000   0.0000  0.0147   0.0147          0.0147 40
D             0.2242   0.2300  0.2166   0.0134          0.2108 40
U            -0.0017  -0.0017  0.0003  -0.0020          0.0003 40
Q             0.2260   0.2317  0.2163   0.0154          0.2106 40
B             0.1472   0.1453  0.1489  -0.0036          0.1508 40
g             1.7824   1.8144  1.7278   0.0866          1.6958 40
gp            0.2225   0.2243  0.2189   0.0054          0.2171 40

The Nagelkerke R-squared for the main effect logistic regression model in the entire dataset is 0.2835. The C-statistic corrected for optimism is 0.5+(0.5859/2) = 0.79295

This is a nomogram for the main effects model.

plot(nomogram(modY_lrm))

There are two new subjects of interest, Beverage A and Beverage B, which are the same size and have the same type of milk and milligrams of caffeine, but Beverage A has 18 grams of sugar and Beverage B has 49 grams of sugar. Beverage B’s odds of having whipped cream will be 8.164 times Beverage A’s odds of having whipped cream. The odds of Beverage B having whipped cream is 816.4% as large as the odds of Beverage A having whipped cream.

There are two new subjects of interest, Beverage C and Beverage D, which are the same size and have the same type of milk and grams of sugar, but Beverage C has 30 mg of caffeine and Beverage D has 150 mg of caffeine. Beverage D’s odds of having whipped cream will be 0.868 times Beverage C’s odds of having whipped cream. The odds of Beverafe D having whipped cream is 86.8% as large as the odds of Beverage C having whipped cream.