Dense bundles of multicoloured fibres sweeping across a near-black ground

Week: Regression & Response Surfaces

Workshop 10: Regression in R

Modeling for System Improvement

Tim Fraser, PhD

Assistant Teaching Professor, Systems Engineering
Cornell University

SYSEN 5300: Systems Engineering & Six Sigma

for Design and Operation of Reliable Systems

Risk Analysis
- Risk assessment and risk characterization
- Failure Modes and Effects Analysis
- Fault trees and event trees under uncertainty
Modeling Reliability
- Component and system reliability
- Physical acceleration models
Quality Control
- Six sigma and statistical process control
- Catch performance slipping, in time to act
System Improvement
- Design of experiments and response surfaces
- Optimize system design for reliability
- Today: find a best recipe, and predict it with a 95% CI

Today's Class

Today's data
Load the data, write tidier()
Group Activity: Explain It
Response Surface Exercises
Coding Demo

The big question: As X increases, what happens to Y?

Animated scatterplot titled Correlation: -1: a cloud of points starts as a perfect negative line, then shifts to a positive relationship as the line of best fit swings with it

Cities in Tohoku, Japan: some hit by the tsunami, some not
X: how much damage a city took
Y: social capital — do neighbours trust and help each other?
As X increases, what happens to Y?

Aerial view from a U.S. Navy helicopter of flooded neighbourhoods and fields around Sendai the day after the 2011 tsunami, a smoke plume rising in the distance

Photo: U.S. Navy, Sendai, 12 March 2011 (public domain, via Wikimedia Commons)

Reading Questions

Workshop 10

Volunteers in white hard hats and work gloves hammer nails into lumber on a house-building site, a hammer in the foreground

Photo: Brad Weaver on Unsplash

Modeling Exercises

Load Data

Read in the data and set the categories, then glimpse it: what type is each column?

You'll need these packages….

install.packages(c(
  "dplyr", "readr", "broom",
  "texreg", "gtools"))

pip install pandas
  statsmodels plotnine

library(dplyr)
library(readr)
library(broom)
library(texreg)
library(gtools)

cities = read_csv("workshops/jp_matching_experiment.csv") %>%
  # Tell R to treat year and pref as **ordered categories**
  mutate(year = factor(year),
         pref = factor(pref),
         by_tsunami = factor(by_tsunami, levels = c("Not Hit", "Hit")))

cities %>% glimpse()

import pandas as pd
import statsmodels.formula.api as smf
# lm() -> smf.ols(...).fit()

f = "workshops/jp_matching_experiment.csv"
cities = pd.read_csv(f)
# Treat year and pref as ordered categories
cities['year'] = cities['year'].astype('category')
cities['pref'] = cities['pref'].astype('category')
cities['by_tsunami'] = pd.Categorical(
  cities['by_tsunami'],
  categories=["Not Hit", "Hit"], ordered=True)

print(cities.info())

Rows: 1,057
Columns: 13
$ muni_code                 <chr> "02201", "02202", "02203", "02205", "02208",…
$ muni                      <chr> "Aomori", "Hirosaki", "Hachinohe", "Goshogaw…
$ pref                      <fct> Aomori, Aomori, Aomori, Aomori, Aomori, Aomo…
$ year                      <fct> 2011, 2011, 2011, 2011, 2011, 2011, 2011, 20…
$ by_tsunami                <fct> Hit, Not Hit, Hit, Hit, Hit, Not Hit, Not Hi…
$ social_capital            <dbl> 0.393, 0.413, 0.408, 0.397, 0.403, 0.387, 0.…
$ damage_rate               <dbl> 0.000, 0.000, 0.001, 0.000, 0.000, 0.000, 0.…
$ pop_density               <dbl> 1119.2, 602.1, 1137.5, 317.6, 432.4, 337.1, …
$ exp_dis_relief_per_capita <dbl> 0.01, 0.06, 2.44, 0.11, 1.93, 0.78, 4.91, 0.…
$ income_per_capita         <dbl> 1.11, 0.99, 1.13, 0.86, 1.10, 0.75, 0.63, 0.…
$ unemployment              <dbl> 5.9, 4.5, 5.5, 5.9, 5.3, 7.5, 3.7, 4.9, 10.5…
$ pop_women                 <dbl> 53.6, 54.1, 52.2, 54.7, 52.0, 53.7, 50.1, 54…
$ pop_over_age_65           <dbl> 27.9, 29.2, 27.5, 31.6, 29.6, 31.4, 37.9, 30…

<class 'pandas.DataFrame'>
RangeIndex: 1057 entries, 0 to 1056
Data columns (total 13 columns):
 #   Column                     Non-Null Count  Dtype   
---  ------                     --------------  -----   
 0   muni_code                  1057 non-null   int64   
 1   muni                       1057 non-null   str     
 2   pref                       1057 non-null   category
 3   year                       1057 non-null   category
 4   by_tsunami                 1057 non-null   category
 5   social_capital             1057 non-null   float64 
 6   damage_rate                1057 non-null   float64 
 7   pop_density                1057 non-null   float64 
 8   exp_dis_relief_per_capita  1057 non-null   float64 
 9   income_per_capita          1057 non-null   float64 
 10  unemployment               1057 non-null   float64 
 11  pop_women                  1057 non-null   float64 
 12  pop_over_age_65            1057 non-null   float64 
dtypes: category(3), float64(8), int64(1), str(1)

Modeling Exercises

Load Function

Let's write a little tidier function…

code/11_workshop_exercises.R

# Let's write a little tidier function..
tidier = function(model, ci = 0.95, digits = 3){
  model %>% # for a model object
    # get data.frame of coefficients
    # ask for a confidence interval matching the 'ci' above!
    broom::tidy(conf.int = TRUE, conf.level = ci) %>%
    # And round and relabel them
    reframe(
      term = term,
      # Round numbers to a certain number of 'digits'
      estimate = estimate %>% round(digits),
      se = std.error %>% round(digits),
      statistic = statistic %>% round(digits),
      p_value = p.value %>% round(digits),
      # Get stars to show statistical significance
      stars = p.value %>% gtools::stars.pval(),
      # Get better names
      upper = conf.high %>% round(digits),
      lower = conf.low %>% round(digits))
}

# Let's write a little tidier function..
def tidier(model, ci = 0.95, digits = 3):
    # get confidence interval matching the 'ci' above!
    b = model.conf_int(alpha = 1 - ci)
    # And round and relabel the coefficients
    out = pd.DataFrame({
      'term': model.params.index,
      'estimate': model.params.round(digits),
      'se': model.bse.round(digits),
      'statistic': model.tvalues.round(digits),
      'p_value': model.pvalues.round(digits),
      'lower': b[0].round(digits),
      'upper': b[1].round(digits)})
    # Get stars to show statistical significance
    out['stars'] = pd.cut(out['p_value'],
      [-0.001, .001, .01, .05, .1, 1],
      labels = ['***', '**', '*', '.', ''])
    return out.reset_index(drop = True)

Group Activity: Explain It

Regression · submit on Canvas

Team Boards

Question 1

Reporting Effects

# 1. Estimate the effect of being hit by the tsunami on income per capita,
# controlling for damage rates.
m = cities %>% lm(formula = income_per_capita ~ damage_rate + by_tsunami)
tidier(m)

# Report the effect of being hit by the tsunami.
# As [X] increases by 1 [unit], [Y] increases by [BETA] [units].
# This effect has a 95% confidence interval from [A] to [B].
# This effect is statistically [significant? insignificant?]
# with a p-value of [XXX].

# 1. Estimate the effect of being hit by the tsunami on income per capita,
# controlling for damage rates.
m = smf.ols("income_per_capita ~ damage_rate + by_tsunami", cities).fit()
tidier(m)

# Report the effect of being hit by the tsunami.
# As [X] increases by 1 [unit], [Y] increases by [BETA] [units].
# This effect has a 95% confidence interval from [A] to [B].
# This effect is statistically [significant? insignificant?]
# with a p-value of [XXX].

Question 10

Intercept when using Categorical Variables

# 10. Model the effect of each year on income.
# Which year is not represented? The intercept represents that baseline category
cities$year %>% unique()
m = cities %>% lm(formula = income_per_capita ~ year)
m
# Calculate the predicted income per capita from 2011 to 2017.

# 10. Model the effect of each year on income.
# Which year is not represented? The intercept represents that baseline category
cities['year'].unique()
m = smf.ols("income_per_capita ~ year", cities).fit()
print(m.params)
# Calculate the predicted income per capita from 2011 to 2017.

Response Surface Methodology

Two exercises with the gingerbread data

Exercise 1

Find the best recipe

This table shows predicted yum from the interaction model at the 16 tested recipes.

Which recipe gives the highest predicted yum, and what is it?

Table: data · code

Molasses (cups) \ Ginger (tbsp)	0.5	1	1.5	2
0.5	12.0	13.8	15.6	17.5
0.75	14.2	14.8	15.4	16.0
1	16.5	15.8	15.2	14.6
1.25	18.7	16.8	15.0	13.2

Molasses (cups) \ Ginger (tbsp)	0.5	1	1.5	2
0.5	12.0	13.8	15.6	17.5
0.75	14.2	14.8	15.4	16.0
1	16.5	15.8	15.2	14.6
1.25	18.7	16.8	15.0	13.2

cookies = read_csv("workshops/gingerbread_test1.csv")
m2 = cookies %>% lm(formula = yum ~ molasses * ginger)
expand_grid(molasses = c(0.5, 0.75, 1, 1.25),
            ginger = c(0.5, 1, 1.5, 2)) %>%
  mutate(yhat = predict(m2, newdata = tibble(molasses, ginger))) %>%
  filter(yhat == max(yhat))

cookies = pd.read_csv("workshops/gingerbread_test1.csv")
m2 = smf.ols("yum ~ molasses * ginger", cookies).fit()
grid = pd.MultiIndex.from_product(
  [[0.5, 0.75, 1, 1.25], [0.5, 1, 1.5, 2]],
  names = ["molasses", "ginger"]).to_frame(index = False)
grid["yhat"] = m2.predict(grid)
grid[grid.yhat == grid.yhat.max()]

# A tibble: 1 × 3
  molasses ginger  yhat
     <dbl>  <dbl> <dbl>
1     1.25    0.5  18.7

1.25 cups molasses, 0.5 tbsp ginger: yum 18.7, at the edge of what we tested.

    molasses  ginger     yhat
12      1.25     0.5  18.6885

1.25 cups molasses, 0.5 tbsp ginger: yum 18.7, at the edge of what we tested.

Exercise 2

Predict with a 95% confidence interval

We never baked 1.1 cups of molasses with 1.3 tbsp of ginger.

Using the second-order model, predict yum for that recipe with a 95% confidence interval.

Second-order model

y = β₀ + β₁x₁ + β₂x₂ + β₃x₁x₂ + β₄x₁² + β₅x₂²

m3 = cookies %>% lm(formula = yum ~ molasses * ginger +
                      I(molasses^2) + I(ginger^2))

tibble(molasses = 1.1, ginger = 1.3) %>%
  mutate(predict(m3, newdata = ., se.fit = TRUE,
                 interval = "confidence", level = 0.95) %>%
           with(fit) %>%
           as_tibble())

m3 = smf.ols("yum ~ molasses * ginger + I(molasses**2)"
             " + I(ginger**2)", cookies).fit()

new = pd.DataFrame({"molasses": [1.1], "ginger": [1.3]})
new.assign(**m3.get_prediction(new).summary_frame(alpha = 0.05)
           [["mean", "mean_ci_lower", "mean_ci_upper"]])

# A tibble: 1 × 5
  molasses ginger   fit   lwr   upr
     <dbl>  <dbl> <dbl> <dbl> <dbl>
1      1.1    1.3  13.7  12.5  14.9

We predict a yum of 13.7 for this recipe, and are 95% confident its average yum lies between 12.5 and 14.9.

   molasses  ginger      mean  mean_ci_lower  mean_ci_upper
0       1.1     1.3  13.70051      12.489732      14.911288

We predict a yum of 13.7 for this recipe, and are 95% confident its average yum lies between 12.5 and 14.9.

Exercise 3

Draw the response surface

Predict yum across a grid of recipes with the second-order model, then plot it.

Pipe the grid into geom_tile() and geom_contour().

Plot: data · code

expand_grid(molasses = seq(0.5, 1.25, length.out = 50),
            ginger   = seq(0.5, 2,    length.out = 50)) %>%
  mutate(yhat = predict(m3, newdata = tibble(molasses, ginger))) %>%
  ggplot() +
  geom_tile(mapping = aes(x = molasses, y = ginger, fill = yhat)) +
  geom_contour(mapping = aes(x = molasses, y = ginger, z = yhat),
               color = "white", binwidth = 1) +
  scale_fill_viridis_c(option = "plasma")

surface = pd.MultiIndex.from_product(
  [np.linspace(0.5, 1.25, 50), np.linspace(0.5, 2, 50)],
  names = ["molasses", "ginger"]).to_frame(index = False)
surface["yhat"] = m3.predict(surface)
# plotnine has no geom_contour(): tiles only
(ggplot(surface) +
  geom_tile(aes(x = "molasses", y = "ginger", fill = "yhat")) +
  scale_fill_cmap("plasma"))

Heatmap of predicted yum over molasses 0.5 to 1.25 cups and ginger 0.5 to 2 tbsp, with white contour lines one yum apart: a low valley near 0.85 cups of molasses, rising toward both edges, highest at 1.25 cups molasses and 0.5 tbsp ginger

Yum bottoms out near 0.85 cups of molasses and rises toward both edges, peaking at 22.7 with 1.25 cups molasses and 0.5 tbsp ginger.

Coding Demo

Workshop 10 [Code]

A crowded workshop room: a man and a woman at a shared table lean toward one laptop while others type on laptops behind them

Photo: fran innocenti on Unsplash

Scatterplot of social capital against damage rate for Tohoku cities, blue points not hit and orange points hit by the tsunami, with a red upward-sloping line of best fit

data · code

Fit a regression in R — price explained by carat — live, in the browser.

Appendix

Workshop 10

Slides not taught

The same fit in Python — smf.ols() instead of lm().