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

Week: Regression

Workshop 11: 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
- Required function, stated conditions, specified time
- One element of the broader risk analysis

Quality Control
- Six sigma and statistical process control
- Detect when performance is deteriorating
- Take corrective action in time
System Improvement
- Optimize system design for reliability
- Design of experiments
- Response surfaces — innovation, problem solving, and discovery

Today's Class

Reading Questions
Coding Demo

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

Animated scatterplot: a cloud of points shifts from a negative to a positive relationship as the line of best fit swings with it

Reading Questions

Workshop 11

Activity: Explain It

Please complete the following steps:

Pick 2 reading questions from slides below.
Copy the slide matching that function, and add it to our Slide Deck [Link]!
Update it with your Group Number and your names.
Fill in the slide with your answers to these questions.
Take a screenshot of your slide and submit it in the CANVAS text submit box!

Question 1 card, Reporting Effects: estimate the effect of being hit by the tsunami on income per capita, controlling for damage rates, with m equals cities piped into lm of formula income per capita against damage rate plus by tsunami, then tidier of m; report the effect using the sentence frame as X increases by 1 unit, Y increases by BETA units, with a 95 percent confidence interval from A to B, statistically significant or insignificant with a p-value of XXX

• • •

Question 10 card, Intercept when using Categorical Variables: model the effect of each year on income; which year is not represented? The intercept represents that baseline category; cities dollar year piped into unique, then m equals cities piped into lm of formula income per capita against year; calculate the predicted income per capita from 2011 to 2017

https://github.com/timothyfraser/sysen/blob/main/code/11_workshop_exercises.R

Load Data

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())

You'll need these packages….

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

pip install pandas statsmodels plotnine

Load Function

# 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 = statistic %>% 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))
}

# The same little tidier function..
# broom::tidy -> .params and .conf_int()
def tidier(model, ci = 0.95, digits = 3):
    # confidence interval matching the 'ci' above
    b = model.conf_int(alpha = 1 - ci)
    out = pd.DataFrame({
      'term': model.params.index,
      # Round numbers to 'digits' places
      'estimate': model.params.round(digits),
      'se': model.bse.round(digits),
      'statistic': model.tvalues.round(digits),
      'p_value': model.pvalues.round(digits),
      # Get better names
      'lower': b[0].round(digits),
      'upper': b[1].round(digits)})
    # Stars to show statistical significance
    out['stars'] = pd.cut(out['p_value'],
      [0, .001, .01, .05, .1, 1],
      labels = ['***', '**', '*', '.', ''])
    return out.reset_index(drop = True)

https://github.com/timothyfraser/sysen/blob/main/code/11_workshop.R

Coding Demo

Workshop 11

https://github.com/timothyfraser/sysen/blob/main/code/11_workshop.R

Fit the multivariate model in R — live, in the browser.

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