Intro to Data Science (Non-Technical Background, R) - Session14

DataKolektiv, Chief Scientist & Owner; Data Scientist for Wikidata, WMDE

---

Session 14: Introduction to Estimation Theory. Multiple Linear Regression. Model diganostics. The role of part correlation in this model. Dummy coding of categorical variables in R. Nested models.

Feedback should be send to goran.milovanovic@datakolektiv.com. These notebooks accompany the Intro to Data Science: Non-Technical Background course 2020/21.

---

What do we want to do today?

In today’s session we need to go beyond the discussion of a simple relationship between two variables: one predictor and one criterion. While Simple Linear Regression is extremely useful in a didactic perspective - to introduce statistical learning methods as such - is is only seldom used in practice. The real world is complex way beyond exploring only simple relationships, and mathematical models in practice almost necessarily deal with many predictors and the existing mutual information among them in an attempt to predict the state of the outcome variable. We will introduce the Multiple Linear Regression model in this session and discuss a more complex scenario involving several predictors. We will also introduce the method of dummy coding when categorical predictors are present in the Multiple Linear Regression scenario. We are also laying ground for even more complex Generalized Linear Models that can handle categorization problems beyond regression. Finally: comparing nested regression models.

0. Prerequisits

Install:

install.packages('QuantPsyc')
install.packages('lattice')

Setup:

dataDir <- paste0(getwd(), "/_data/")
library(tidyverse)
library(Hmisc)
library(ppcor)
library(car)
library(datasets)
library(broom)
library(QuantPsyc)
library(lattice)

1. Multiple Linear Regression: the exposition of the problem

1.1 Iris, again

The `iris’ dataset again offers everything that we need to introduce a new statistical model. One might wonder, but that small and (manifestly!) simple dataset can be used as well to introduce even more complex statistical models than Multiple Linear Regression!

data(iris)
str(iris)

'data.frame':   150 obs. of  5 variables:
 $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
 $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
 $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
 $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
 $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...

And back to the problem with Petal.Length ~ Sepal.Length regression that we have already discussed:

ggplot(data = iris,
       aes(x = Sepal.Length, 
           y = Petal.Length)) +
  geom_point(size = 2, colour = "blue") +
  geom_smooth(method='lm', size = .25) +
  ggtitle("Sepal Length vs Petal Length") +
  xlab("Sepal Length") + ylab("Petal Length") + 
  theme_bw() + 
  theme(panel.border = element_blank()) + 
  theme(plot.title = element_text(hjust = .5))

… where everything looks nice until…

ggplot(data = iris,
       aes(x = Sepal.Length, 
           y = Petal.Length, 
           color = Species, 
           group = Species)) +
  geom_point(size = 2) +
  geom_smooth(method='lm', size = .25) +
  ggtitle("Sepal Length vs Petal Length") +
  xlab("Sepal Length") + ylab("Petal Length") + 
  theme_bw() + 
  theme(panel.border = element_blank()) + 
  theme(plot.title = element_text(hjust = .5))

Not to mention the model assumptions that we have discussed in the previous Session.

Let’s study the problem a bit.

plot(iris[ , c(1,3,5)],
     main = 
       "Inspect: Sepal vs. Petal Length \nfollowing the discovery of the Species...",
     cex.main = .75,
     cex = .6)

Did we ever mention {lattice}, the hero of data visualization in R before {ggplot2}?

# - {latice} xyplot
xyplot(Petal.Length ~ Sepal.Length | Species,
       data = iris,
       xlab = "Sepal Length", 
       ylab = "Petal Length"
)

Consider the conditional densities of Petal.Length given Species (using {lattice} again):

# - {latice} densityplot
densityplot(~ Petal.Length | Species,
            data = iris,
            plot.points = FALSE,
            xlab = "Petal Length", 
            ylab = "Density",
            main = "P(Petal Length|Species)",
            col.line = 'red'
)

And now consider the conditional densities of Sepal.Length given Species:

# - {latice} densityplot
densityplot(~ Sepal.Length | Species,
            data = iris,
            plot.points=FALSE,
            xlab = "Sepal Length", ylab = "Density",
            main = "P(Sepal Length|Species)",
            col.line = 'blue'
)

I have and idea: why not run a series of separate Simple Linear Regressions in the subgroups defined by Species and inspect the results? Let’s do it:

# - setosa
species <- unique(iris$Species)
w1 <- which(iris$Species == species[1])
reg <- lm(Petal.Length ~ Sepal.Length, 
          data = iris[w1,]) 
tidy(reg)

# - versicolor
w2 <- which(iris$Species == species[2])
reg <- lm(Petal.Length ~ Sepal.Length, data = iris[w2,]) 
tidy(reg)

# - virginica
w3 <- which(iris$Species == species[3])
reg <- lm(Petal.Length ~ Sepal.Length, data = iris[w3,]) 
tidy(reg)

I have used broom::tidy to tidy up the model summaries. The {broom} package offers many useful functions to deal with potentially messy outputs of R’s modeling functions such as lm().

So, Species obviously has some effect on Petal.Length, and that effect possibly goes even beyond the effect of Sepal.Length. How do we incorporate another predictor into the regression model?

1.2 The predictor is categorical: Dummy Coding

We will now try to predict Petal.Length from Species alone in a Simple Linear Regression model. First:

is.factor(iris$Species)

[1] TRUE

Ok, and the levels?

levels(iris$Species)

[1] "setosa"     "versicolor" "virginica" 

Regression with one categorical predictor:

reg <- lm(Petal.Length ~ Species, 
          data = iris) 
tidy(reg)

What effects are present? Let’s see: Speciesversicolor, Speciesvirginica, ok, but what happened to Speciessetosa..? It is our baseline, see:

levels(iris$Species)

[1] "setosa"     "versicolor" "virginica" 

The broom::glance() function is similar to summary() but gives us the model overview all tidy:

broom::glance(reg)

Never forget what the regression coefficient of a dummy variable means: it tells us about the effect of moving from the baseline towards the respective reference level. Here: baseline = setosa (cmp. levels(iris$Species) vs. the output of tidy(reg)). Hence: always look after the order of levels in linear models!

For example, we can change the baseline in Species to versicolor:

# - Levels: setosa versicolor virginica
levels(iris$Species)

[1] "setosa"     "versicolor" "virginica" 

iris$Species <- factor(iris$Species, 
                       levels = c("versicolor", 
                                  "virginica",
                                  "setosa")
                       )
levels(iris$Species)

[1] "versicolor" "virginica"  "setosa"    

Regression again:

# - baseline is now: versicolor
reg <- lm(Petal.Length ~ Species, 
          data = iris) 
tidy(reg)

1.3 Understanding dummy coding

Here is another way to perform dummy coding of categorical variables in R:

# - ...just to fix the order of Species back to default
rm(iris); data(iris)
levels(iris$Species)

[1] "setosa"     "versicolor" "virginica" 

In order to understand what dummy coding really is:

contr.treatment(3, base = 1)

  2 3
1 0 0
2 1 0
3 0 1

And then specifically applied to iris$Species:

contrasts(iris$Species) <- contr.treatment(3, base = 1)
contrasts(iris$Species)

           2 3
setosa     0 0
versicolor 1 0
virginica  0 1

Do not forget that:

class(iris$Species)

[1] "factor"

Now let’s play with level ordering:

iris$Species <- factor(iris$Species, 
                       levels = c ("virginica", 
                                   "versicolor", 
                                   "setosa"))
levels(iris$Species)

[1] "virginica"  "versicolor" "setosa"    

contrasts(iris$Species) = contr.treatment(3, base = 1)
# - baseline is now: virginica
contrasts(iris$Species)

           2 3
virginica  0 0
versicolor 1 0
setosa     0 1

# - consider carefully what you need to do

levels(iris$Species)

[1] "virginica"  "versicolor" "setosa"    

2. Multiple Linear Regression: the problem solved

2.1 One categorical + one continuous predictor

Now we run a multiple linear regression model with Sepal.Length and Species (dummy coded) as predictors of Petal.Length:

# - Petal.Length ~ Species (Dummy Coding) + Sepal.Length 
rm(iris); data(iris) # ...just to fix the order of Species back to default
reg <- lm(Petal.Length ~ Species + Sepal.Length, 
          data = iris)
tidy(reg)

glance(reg)

N.B. Since is.factor (iris$Species) == T - R does the dummy coding in lm() internally for us!

Let’s now compare these results with the simple linear regression model:

reg <- lm(Petal.Length ~ Sepal.Length, data=iris) 
tidy(reg)

glance(reg)

2.2 Nested models

We will now specify two regression models: reg1 defined as Petal.Length ~ Sepal.Length and reg2 defined as Petal.Length ~ Species + Sepal.Length. Obviously, reg2 encompasses reg1 in some way, right? Of course: the predictors in one model are a subset of predictors in another. Such models are called nested models. In this terminological game, reg2 would also be called a full model: a terminology will be used quite often in Binary Logistic Regression, the first Generalized Linear Model that we will meet in our next session.

Note on nested models: There is always a set of coefficients for the nested model (e.g. reg1) such that it can be expressed in terms of the full model (reg2). Can you figure it out?

# - reg1 is nested under reg2
reg1 <- lm(Petal.Length ~ Sepal.Length, 
           data = iris)
reg2 <- lm(Petal.Length ~ Species + Sepal.Length, 
           data = iris)

We can use the partial F-test to compare nested models:

anova(reg1, reg2) # partial F-test; Species certainly has an effect beyond Sepal.Length

Analysis of Variance Table

Model 1: Petal.Length ~ Sepal.Length
Model 2: Petal.Length ~ Species + Sepal.Length
  Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
1    148 111.459                                  
2    146  11.657  2    99.802 624.99 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

2.3 Model diagnostics

We can use the same kind of Influence Plot to search for influential cases in Multiple Linear Regression as we did in the case of Simple Linear Regression (except for this time the computation of the relevant indicators is way more complicated):

infReg <- as.data.frame(influence.measures(reg)$infmat)
head(infReg)

This time we will use broom:augment() to obtain the influence measures:

regFrame <- broom::augment(reg2)
head(regFrame)

Produce the Influence Chart:

plotFrame <- data.frame(residual = regFrame$.std.resid,
                        leverage = regFrame$.hat,
                        cookD = regFrame$.cooksd)
ggplot(plotFrame,
       aes(y = residual,
           x = leverage)) +
  geom_point(size = plotFrame$cookD * 100, 
             shape = 1, color = "blue") +
  ylab("Standardized Residual") + 
  xlab("Leverage") +
  ggtitle("Influence Plot\nSize of the circle corresponds to Cook's distance") +
  theme_bw() + 
  theme(panel.border = element_blank()) + 
  theme(plot.title = element_text(size = 8, 
                                  face = "bold", 
                                  hjust = .5))

3. Several continuous predictors

3.1 The stackloss problem

The following example is a modification of the multiple-linear-regression section from R Tutorial.

data(stackloss)
str(stackloss)

'data.frame':   21 obs. of  4 variables:
 $ Air.Flow  : num  80 80 75 62 62 62 62 62 58 58 ...
 $ Water.Temp: num  27 27 25 24 22 23 24 24 23 18 ...
 $ Acid.Conc.: num  89 88 90 87 87 87 93 93 87 80 ...
 $ stack.loss: num  42 37 37 28 18 18 19 20 15 14 ...

The description of the stackloss dataset is found in the documentation:

• Water Temp is the temperature of cooling water circulated through coils in the absorption tower;
• Air Flow is the flow of cooling air;
• Acid Conc. is the concentration of the acid circulating;
• stack.loss (the outcome variable) is 10 times the percentage of the ingoing ammonia to the plant that escapes from the absorption column unabsorbed; that is, an (inverse) measure of the overall efficiency of the plant.

stacklossModel = lm(stack.loss ~ Air.Flow + Water.Temp + Acid.Conc., 
                    data = stackloss)
summary(stacklossModel)

Call:
lm(formula = stack.loss ~ Air.Flow + Water.Temp + Acid.Conc., 
    data = stackloss)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.2377 -1.7117 -0.4551  2.3614  5.6978 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -39.9197    11.8960  -3.356  0.00375 ** 
Air.Flow      0.7156     0.1349   5.307  5.8e-05 ***
Water.Temp    1.2953     0.3680   3.520  0.00263 ** 
Acid.Conc.   -0.1521     0.1563  -0.973  0.34405    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.243 on 17 degrees of freedom
Multiple R-squared:  0.9136,    Adjusted R-squared:  0.8983 
F-statistic:  59.9 on 3 and 17 DF,  p-value: 3.016e-09

glance(stacklossModel)

tidy(stacklossModel)

Prediction for one single new data point:

# predict new data
obs = data.frame(Air.Flow = 72, 
                 Water.Temp = 20, 
                 Acid.Conc. = 85)
predict(stacklossModel, obs)

       1 
24.58173 

The confint() functions works as usual, for 95% CI…

confint(stacklossModel, level = .95) # 95% CI

                  2.5 %      97.5 %
(Intercept) -65.0180339 -14.8213150
Air.Flow      0.4311143   1.0001661
Water.Temp    0.5188228   2.0717495
Acid.Conc.   -0.4818741   0.1776291

… as well as for %99 CI:

confint(stacklossModel, level = .99) # 99% CI

                  0.5 %     99.5 %
(Intercept) -74.3970156 -5.4423333
Air.Flow      0.3247901  1.1064903
Water.Temp    0.2286670  2.3619053
Acid.Conc.   -0.6050987  0.3008536

3.2 Multicolinearity in Multiple Regression

That crazy thing with multiple regression: if the predictors are not correlated at all, why not run a series of simple linear regressions? On the other hand, if the predictors are highly correlated, problems with the estimates arise… John Fox’s {car} package allows us to compute the Variance Inflation Factor quite easily:

VIF <- vif(stacklossModel)
VIF

  Air.Flow Water.Temp Acid.Conc. 
  2.906484   2.572632   1.333587 

The Variance Inflation Factor (VIF) measures the increase in the variance of a regression coefficient due to colinearity. It’s square root (sqrt(VIF)) tells us how much larger a standard error of a regression coefficient is compared to a hypothetical situation where there were no correlations with any other predictors in the model. NOTE: The lower bound of VIF is 1; there is no upper bound, but VIF > 2 typically indicates that one should be concerned.

sqrt(VIF)

  Air.Flow Water.Temp Acid.Conc. 
  1.704841   1.603943   1.154811 

3.3 Part correlation in multiple regression

In multiple regression, it is the semi-partial (or part) correlation that you need to inspect:

• assume a model with X1, X2, X3 as predictors, and Y as a criterion;

• you need a semi-partial of X1 and Y following the removal of X2 and X3 from Y;

• it goes like this: in Step 1, you perform a multiple regression Y ~ X2 + X3;

• In Step 2, you take the residuals of Y, call them RY;

• in Step 3, you regress (correlate) RY ~ X1: the correlation coefficient that you get from Step 3 is the part correlation that you’re looking for!

Recall our model…

stacklossModel = lm(stack.loss ~ Air.Flow + Water.Temp + Acid.Conc.,
                    data = stackloss)
summary(stacklossModel)

Call:
lm(formula = stack.loss ~ Air.Flow + Water.Temp + Acid.Conc., 
    data = stackloss)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.2377 -1.7117 -0.4551  2.3614  5.6978 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -39.9197    11.8960  -3.356  0.00375 ** 
Air.Flow      0.7156     0.1349   5.307  5.8e-05 ***
Water.Temp    1.2953     0.3680   3.520  0.00263 ** 
Acid.Conc.   -0.1521     0.1563  -0.973  0.34405    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.243 on 17 degrees of freedom
Multiple R-squared:  0.9136,    Adjusted R-squared:  0.8983 
F-statistic:  59.9 on 3 and 17 DF,  p-value: 3.016e-09

What is the semi-partial (part correlation) of stack.loss and Air.Flow? Remember the {pcorr} package?

spartCor1 <- spcor.test(x = stackloss$Air.Flow, 
                        y = stackloss$stack.loss,
                        z = stackloss[ , c("Water.Temp", "Acid.Conc.")],
                        method = "pearson")
print(spartCor1)

The unique contribution of Air.Flow is then:

spartCor1$estimate

[1] 0.4631864

spartCor1$statistic

[1] 2.154858

spartCor1$p.value

[1] 0.04580372

---

Further Readings

• Regression, by David M. Lane
• {broom} package: Vignette

R Markdown

• R Markdown is what I have used to produce this beautiful Notebook. We will learn more about it near the end of the course, but if you already feel ready to dive deep, here’s a book: R Markdown: The Definitive Guide, Yihui Xie, J. J. Allaire, Garrett Grolemunds.

---

Goran S. Milovanović

DataKolektiv, 2020/21

contact: goran.milovanovic@datakolektiv.com

---

License: GPLv3 This Notebook is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This Notebook is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this Notebook. If not, see http://www.gnu.org/licenses/.

---