Fundamental Techniques in Data Science with R Notes

## Load the 'diabetes' data from the tab-delimited file
## '../../../data/diabetes.txt'
diabetes <- read.table(paste0(dataDir, "diabetes.txt"),
    header = TRUE,
    sep = "\t")

write.table(boys,
    paste0(dataDir, "boys.txt"),
    row.names = FALSE,
    sep = "\t",
    na = "-999")

## Load the 2017 UTMB data from the comma-separated file
## '../../../data/utmb_2017.csv'
utmb1 <- read.csv(paste0(dataDir, "utmb_2017.csv"))

write.csv2(boys, paste0(dataDir, "boys.csv"), row.names = FALSE, na = "")

• The read.csv() function assumes the values are seperated by commas
• For EU-formatted CSV files (semicolons instead of commas), we can use read.csv2()

Reading this format is tricky, if we want to read SAV files there are two popular options:

• foreign::read.spss()
• haven::read_spss()

foreign doesn’t have a read equivelant

## Load the foreign package:
library(foreign)

## Use foreign::read.spss() to read '../../../data/mtcars.sav' into a list
mtcars1 <- read.spss(paste0(dataDir, "mtcars.sav"))

## Read '../../../data/mtcars.sav' as a data frame
mtcars2 <- read.spss(paste0(dataDir, "mtcars.sav"), to.data.frame = TRUE)

## Read '../../../data/mtcars.sav' without value labels
mtcars3 <- read.spss(paste0(dataDir, "mtcars.sav"),
    to.data.frame = TRUE,
    use.value.labels = FALSE)

library(labelled)

## Use haven::read_spss() to read '../../../data/mtcars.sav' into a tibble
mtcars4 <- read_spss(paste0(dataDir, "mtcars.sav"))

haven::read_spss() converts any SPSS variables with labels into labelled vectors

• We can use the labbeled::unlabeleld() function to remove the value labels
• [V-Sha~ -> V-Shaped

The best option we have for writing SPSS data is haven::write_sav()

write_sav(mtcars2, paste0(dataDir, "mctars2.sav"))

This will preserve label information provided by factor variables and the haven_labelled class

We have two good options for loading data from Excel spreadsheets:

• readxl::read_excel()
• openxlsx::read.xlsx()

## Load the packages:
library(readxl)
library(openxlsx)

## Use the readxl::read_excel() function to read the data from the 'titanic'
## sheet of the Excel workbook stored at '../../../data/example_data.xlsx'
titanic2 <- read_excel(paste0(dataDir, "example_data.xlsx"),
    sheet = "titanic")

## Use the openxlsx::read.xlsx() function to read the data from the 'titanic'
## sheet of the Excel workbook stored at '../../../data/example_data.xlsx'
titanic3 <- read.xlsx(paste0(dataDir, "example_data.xlsx"),
    sheet = "titanic")

The openxlsx package provides a powerful toolkit for programmatically building Excel workbooks in R and saving the results.

## Use the openxlsx::write.xlsx() function to write the 'diabetes' data to an
## XLSX workbook
write.xlsx(diabetes, paste0(dataDir, "diabetes.xlsx"), overwrite = TRUE)

## Use the openxlsx::write.xlsx() function to write each data frame in a list

## to a separate sheet of an XLSX workbook
write.xlsx(list(titanic = titanic, diabetes = diabetes, mtcars = mtcars),
    paste0(dataDir, "example_data.xlsx"),
    overwrite = TRUE)

We can define our own functions using the “function()” function

square <- function(x) {
    out <- x^2
    out
}

square(5)

25

One-line functions don’t need braces:

square <- function(x) x^2
square(5)

25

There are three types of loops in R: for, while and until. Deze course behandelt alleen for (lol)

for(INDEX in RANGE) { Stuff To Do with the Current INDEX Value }

examples:

val <- 0
for(i in 1:100) {
    val <- val + i
}
val

5050

This loop will compute the mean of every column in the mtcars data

means <- rep(0, ncol(mtcars))
for(j in 1:ncol(mtcars)) {
    means[j] <- mean(mtcars[ , j])
}
means

listThing <- list(6,9,4,2)
sum <- 0
for(i in listThing) {
    sum <- i + sum
}
sum

21

Usually there is a function to achieve what you are trying to do with the loop (sum() for example)

Usually one of two forms:

apply(DATA, MARGIN, FUNCTION, ...)
apply(DATA, FUNCTION, ...) # not sure if this works??? but was in slides

Vectors come in one of six “atomic modes”:

• numeric/double
• logical -> bool
• character
• integer
• complex
  • A complex value in R is defined via the pure imaginary value i.
• raw
  • Intended to hold raw bytes

Not possible, some data types will be converted

a <- c(1,2,3)
b <- c("foo", "bar")
c(a,b)

[1] “1” “2” “3” “4” “5” “foo” “bar”

matrix(1:5, nrow = 5, ncol = 2)

cbind(c(1,2,3),c(4,5,6))

m1 <- matrix(1:5, nrow = 5, ncol = 2)
attributes(m1)

Matrices are populated in column-major order, by default

matrix(1:9, 3, 3)

The byrow = TRUE option allows us to fill by row-major order

matrix(1:9, 3, 3, byrow = TRUE)

Like vectors, matrices can only hold one type of data

l3 <- list(name = "bob",
    alive = TRUE,
    age = 33,
    relationshipStatus = 42+3i)

l3$name

bob

Can also be done post-hoc via the names()

l1 = list(3, "hi", TRUE)
names(l1) <- c("name1", "second", "last")
l1$second

hi

(l4 <- list())
l4$people <- c("Bob", "Alice", "Suzy")
l4$money <- 0
l4$logical <- FALSE
l4

$people [1] “Bob” “Alice” “Suzy” $money [1] 0 $logical [1] FALSE

We can create one using data.frame()

data.frame(1:6, c(-1, 1), seq(0.1, 0.6, 0.1))

Columns can also be given names:

data.frame(x = 1:6, y = c(-1, 1), z = seq(0.1, 0.6, 0.1))

x y z 1 1 -1 0.1 2 2 1 0.2 3 3 -1 0.3 4 4 1 0.4 5 5 -1 0.5 6 6 1 0.6

Creating from matrix

data.frame(matrix(5, 4, 3))

Sample makes a vector of some length randomly sampled from a vector

Runif generates a vector of random numbers

data.frame(a = sample(c(TRUE, FALSE), 5, replace = TRUE),
    b = sample(c("foo", "bar"), 5, replace = TRUE),
    c = runif(5)
    )

x <- factor(sample(1:3, 8, TRUE), labels = c("red", "yellow", "blue"))
attributes(x)

$levels [1] “red” “yellow” “blue” $class [1] “factor”

Even though a factor’s data are represented by an integer vector, R does not consider factors to be interger/numeric data.

Rarely used, does have some differences though, drops any names or dimnames attribute, and only allows selection of a single element

x <- 10:20
x [[2]]

11

The [] method returns objects of class list (or data.frame if foo was a data.frame) while the [[]] method returns objects whose class is determined by the type of their values.

y <- matrix(1:24, 6, 4)
y[4:6,2:3] # 2

We can mix different selecting styles

y <- matrix(1:24, 6, 4)
y[c(2,4),c(FALSE, TRUE, TRUE)] # 2

We can use all three operators to subset lists.

list1 <- list(2, "noo", TRUE)
names(list1) <- c("first", "second", "last")

list1$first

list1[1]

list1[["first"]]

2

l4 <- list(2, "noo", TRUE)
names(l4) <- c("first", "second", "last")
tmp1 <- l4[1]
tmp1$people
# [1] "Bob" "Alice" "Suzy"
class(tmp1)
# [1] "list"
tmp2 <- l4[[1]]
# [1] "Bob" "Alice" "Suzy"
class(tmp2)
# [1] "character"

numeric

We can subset the columns of a data frame using list semantics

d3 <- data.frame(1:6, c(-1, 1), seq(0.1, 0.6, 0.1))
names(d3) <- c("a", "b", "c")
d3$a

d3 <- data.frame(1:6, c(-1, 1), seq(0.1, 0.6, 0.1))
names(d3) <- c("a", "b", "c")
d3[1]

Filter by rows

d3 <- data.frame(1:6, c(-1, 1), seq(0.1, 0.6, 0.1))
names(d3) <- c("a", "b", "c")
d3[(d3$a)%%2==0 & d3$c>0,]

Notic the comma at the end!

This comma is needed because R uses the following notation for subsetting data frames.

data(rows you want, columns you want)

Example use:

library(dplyr)

sqr <- function (x)
{
    x**2
}

plus1 <- function (x)
{
    x + 1
}

test <-
    (1:10) %>%
    sqr() %>%
    plus1()

test

This is equivalent to:

test <-
    sqr(
        plus1(
            (1:10)))

test

Specify where you want the piped value to go

library(dplyr)

power <- function (x, y)
{
    x**y
}

test <-
    (1:10) %>%
    power(3, y=.)

test

[%$%] expose the names in lhs to the rhs expression. This is useful when functions do not have a built-in data argument. or The exposition pipe exposes the contents of an object to the next function in the pipeline.

library(dplyr)

# all equivelant:

cats %$% t.test(Hwt ~ Sex)
cats %>% t.test(Hwt ~ Sex, data = .)
t.test(Hwt ~ Sex, data = cats)

We can create a basic scatterplot using the plot() function

plot(runif(20, 0, 100), runif(20, 0, 10))

diabetes %$% plot(
    y = tc,
    x = bmi,
    ylab = "Total Cholesterol",
    xlab = "Body Mass Index",
    main = "Relation between BMI and Cholesterol",
    ylim = c(0, 350),
    xlim = c(0, 50)
    )

We can create a simple histogram with the hist() function

hist(runif(30,0,10))

We can create simple boxplots via the boxplot function

boxplot(runif(30,40,60),runif(60,10,50))

We start by initializing the object that will define our plot

• This is regardless of what we want to plot
• We must specify a data source and a aesthetic via the aes() function

library(ggplot2)
p1 <- ggplot(data = diabetes, mapping = aes(x = bmi, y = glu))
p1

This doesn’t make any plot yet

We can also add statistical summaries of the data.

p1 + geom_point() + geom_smooth()

p1 + geom_point() + geom_smooth(method = "lm")

Changing the style outside of aes() applies the styling to the entire plot

We can apply canned themes to adjust a plot overall appearance\ Default theme:

p1.1 <- p1 + geom_point()

p1.1 + theme_classic()

You can also modify individual themes

p1.1 + theme_classic() +
    theme(axis.title = element_text(size = 16,
                                    family = "serif",
                                    face = "bold",
                                    color = "blue"),
        aspect.ratio = 1)

Faceting allow us to make arrays of conditional plots.

p7 <- ggplot(titanic, aes(age, survived, color = class)) +
    geom_jitter(height = 0.05) +
    geom_smooth(method = "glm",
                method.args = list(family = "binomial"),
                se = FALSE)
p7 + facet_wrap(vars(sex))

p8 <- ggplot(titanic, aes(age, survived)) +
    geom_jitter(height = 0.05) +
    geom_smooth(method = "glm",
                method.args = list(family = "binomial"),
                se = FALSE)
p8 + facet_grid(vars(sex), vars(class))

If we want to paste several different plots into a single figure (without faceting), we can use the utilities in the gridExtra package.

library(gridExtra)

grid.arrange(p1 + geom_point(),
             p3 + geom_boxplot(),
             p1 + geom_line(),
             p8 + facet_grid(vars(sex), vars(class)),
             ncol = 2)

Missing data patterns represent unique combinations of observerd and missing items

• \(P\) items -> \(2^P\) possible patterns

The ggmice package provides some nice ways to visualize incomplete data and objects created during missing data treatment.

library(ggmice); library(ggplot2)
ggmice(boys, aes(wgt, hgt)) + geom_point()

We can also create a nicer version of the response pattern plot.

plot_pattern(boys, rotate = TRUE)

The naniar package also provides some nice visualization and numerical summary routines for incomplete data.

naniar::gg_miss_upset(boys)

Tukey (1977) described a procedure for flagging potential outliers based on a box-and-whiskers plot.

• Does not require normally distributed X
• Not sensitive to outliers

A fence is an interval defined as the following function of the first quartile, the third quartile, and the inner quartile range \((IQR=Q_3-Q_1)\): \[F=\{Q_1-C\times IQR,\ Q_3+C\times IQR\}\]

• Taking C = 1.5 produces the inner fence
• Taking C = 3.0 produces the outer fence

We can use these fences to indentify potential outliers:

• Any value that falls outside of the inner fence is a possible outlier
• Any value that falls outside of the inner fence is a probable outlier

To compare statistics, we consider their breakdown points.

• The breakdown point is the minimum proportion of cases that must be replaced by ∞ to cause the value of the statistic to go to ∞.

The mean has a breakdown point of 1/N.

• Replacing a single value with ∞ will produce an infinite mean.
• This is why we shouldn’t use basic z-scores to find outliers.

The median has breakdown point of 50%.

• We can replace n < N/2 of the observations with ∞ without producing an infinite median.

Nominal, ordinal, and binary items can have outliers

• Outliers on categorical variables are often more indicative of bad variables than outlying cases.

Ordinal

• Most participant endorse one of the lowest categories on an ordinal item, but a few participants endorse the highest category.
• The participants who endorse the highest category may be outliers.

Nominal

• Groups with very low membership may be outliers on nominal grouping variables.

Binary

• If most endorse the item, the few who do not may be outliers.

If we locate any outliers, they must be treated.

• Outliers cause by errors, mistakes, or malfunctions (i.e., error outliers) should be directly corrected.
• Labeling non-error outliers is a subjective task.
  • A (non-error) outlier must originate from a population separate from the one we care about.
  • Don’t blindly automate the decision process.

The most direct solution is to delete any outlying observation.

• If you delete non-error outliers, the analysis should be reported twice: with outliers and without.

For univariate outliers, we can use less extreme types of deletion.

• Delete outlying values (but not the entire observation).
• These empty cells then become missing data.
• Treat the missing values along with any naturally-occurring nonresponse.

Winsorization:

• Replace the outliers with the nearest non-outlying value

We can also use robust regression procedures to estimate the model directly in the presence of outliers.

Weight the objective function to reduce the impact of outliers

• M-estimation

Trim outlying observations during estimation

• Least trimmed squares, MCD, MVE

Take the median, instead of the mean, of the squared residuals

• Least median of squares

Model some quantile of the DV’s distribution instead of the mean

• Quantile regression

Model the outcome with a heavy-tailed distribution

• Laplacian, Student’s T

The above equation only describes the best fit line. We need to account for the estimation error

Minimize the sum of squared residuals

\(\hat \epsilon_n=Y_n-\hat Y_n\)

We may also want to know how well our model explains the outcome.

• Our model explains some proportion of the outcome’s variability.
• The residual variance ˆ𝜎2 = Var(ˆ𝜀) will be less than Var(Y)

We quantify the proportion of the outcome’s variance that is explained by our model using the R2 statistic

\(\bar y\) is average of y

When assessing predictive performance, we will most often use the mean squared error as our criterion.

MSE is better to use when we want to predict, so we dont care about explenation.

\[RMSE = \sqrt{MSE}\]

We can use information criteria to quickly compare non-nested models while accounting for model complexity.

Information criteria balance two competing forces.

• Blue: The optimized loglikelihood quantifies fit to the data.
• Red: The penalty term corrects for model complexity

K number of parameters, amount of x variables in this case

These statistics are only interesting for comparing models, on their on they don’t say a lot

All R factors have an associated contrasts attribute. • The contrasts define a coding to represent the grouping information. • Modeling functions code the factors using the rules defined by the contrasts.

x <- factor(c("a", "b", "c"))
contrasts(x)

For variables with only two levels, we can test the overall factor’s significance by evaluating the significance of a single dummy code.

out <- lm(Price ~ Man.trans.avail, data = Cars93)
partSummary(out, -c(1, 2))
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 23.841 1.623 14.691 <2e-16
# Man.trans.availYes -6.603 2.004 -3.295 0.0014
# Residual standard error: 9.18 on 91 degrees of freedom
# Multiple R-squared: 0.1066,Adjusted R-squared: 0.09679
# F-statistic: 10.86 on 1 and 91 DF, p-value: 0.001403

For variables with more than two levels, we need to simultaneously evaluate the significance of each of the variable’s dummy codes

partSummary(out4, -c(1, 2))
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 21.7187 2.9222 7.432 6.25e-11
# Man.trans.availYes -5.8410 1.8223 -3.205 0.00187
# DriveTrainFront -0.2598 2.8189 -0.092 0.92677
# DriveTrainRear 10.5169 3.3608 3.129 0.00237
# Residual standard error: 8.314 on 89 degrees of freedom
# Multiple R-squared: 0.2834,Adjusted R-squared: 0.2592
# F-statistic: 11.73 on 3 and 89 DF, p-value: 1.51e-06

summary(out4)$r.squared - summary(out)$r.squared
# [1] 0.1767569
# anova(out, out4)
# Analysis of Variance Table
# Model 1: Price ~ Man.trans.avail
# Model 2: Price ~ Man.trans.avail + DriveTrain
# Res.Df RSS Df Sum of Sq F Pr(>F)
# 1 91 7668.9
# 2 89 6151.6 2 1517.3 10.976 5.488e-05 ***
# ---
# Signif. codes:
# 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

To quantify uncertainty in our predictions, we want to use an appropriate interval estimate: There are two flavors:

In addtive MLR we have the following equation: \(Y=\beta_0+\beta_1X+\beta_2Z+\epsilon\) This equation assumes that X and Z are independent predictors of Y.

When X and Z are independent predictors, the following are true

• X and Z can be correlated
• \(\beta_1\) and \(\beta_2\) are partial regression coefficients
• The effect of X and Y is the same at all levels of Z, and the effect of Z on Y is the same at all levels of X.

When testing the moderation, we hypothesize that the effect of X on Y varies af a function of Z.

• \(Y=\beta_0+f(Z)X+\beta_2Z+\epsilon\)

If we assume that Z linearly affects the relationship between X and Y:

• \(f(Z)=\beta_1+\beta_3Z\)

Combining the two equations:

• \(Y=\beta_0+(\beta_1+\beta_3Z)X+\beta_2Z+\epsilon\)
• \(Y=\beta_0+\beta_1X+\beta_2Z+\beta_3XZ+\epsilon\)

\(Y=\beta_0+\beta_1X+\beta_2Z+\beta_3XZ+\epsilon\)

To test for moderation we simply need to test the significan of the interaction term, XZ

\(\hat\beta_3\) quantifies the effect of Z on the focal effect (the X->Y effect)

• For a unit change in Z, it is the expected change in the effect of X and Y

\(\hat\beta_1\text{ and } \hat\beta_2\) are conditional effects.

• interpreted where the other predictor is zero
• For a unit change in X, ˆ𝛽1 is the expected change in Y, when Z = 0.
• For a unit change in Z, ˆ𝛽2 is the expected change in Y, when X = 0.

This is OK: \(Y=\beta_0+\beta_1X+\beta_2Z+\beta_3XZ+\beta_4X^2+\beta_5X^3\)

This is not: \(Y=\beta_0X^{\beta_1}+\epsilon\)

• power is also a parameter here -> not allowed

If the model is not linear in the parameters, then we’re not even working with linear regression.

• We need to move to entirely different modeling paradigm

Each modeled X must exhibit a linear relation with Y.

• We can define X via nonlinear transformations of the original data.

• \(N>P\)
• No \(X_p\) can be a linear combination of other predictors

If the predictor matrix is not full rank, the model is not estimable.

• The parameter estimates cannot be uniquely determined from the data

The predictors do not correlate with the errors

• \(Cov(\hat Y,\epsilon)=0\)
• \(E[\epsilon_n]=0\)

If the predictors are not exogenous, the estimated regression coefficients will be biased.

• \(Var(\epsilon_n)=\sigma^2<\infty\)

see spherical

\(\Cov(\epsilon_j,\epsilon_j)=0,\ i\ne j\)

see spherical

\(e\sim N(0,\sigma^2)\)

If errors are non-normal, small-sample inferences may be biased.

• The justification for some tests and procedures used in regression analysis may not hold

The assumption of spherical erros combines assumptions 4 and 5 \(Var(\epsilon)=\sigma^2I_N\)

If the errors are not spherical, the standard errors will be biased.

• The estimated regression coefficients will be unbiased, though.
• But the statistics and stuf can be biased

One of the most useful diagnostic graphics is the plot of residuals vs. predicted values.

We can create this plot by plotting the fitted lm() object in R.

out1 <- lm(Price ~ Horsepower,
    data = Cars93)
plot(out1, 1)

Outliers can be identified by scrutinizing the residuals.

• Observations with residuals of large magnitude may be outliers.
• The difficulty arises in quantifying what constitutes a “large” residual.

If the residuals do not have constant variance, then we cannot directly compare them.

• We need to standardize the residuals in some way

We are specifically interested in externally studentized residuals.

We can’t simply standardize the ordinary residuals.

• Internally studentized residuals
• Outliers can pull the regression line towards themselves.
• The internally studentized residuals for outliers will be too small.

Begin by defining the concept of a deleted residual: \[ \hat \epsilon_{(n)}=Y_n-\hat Y_{(n)} \]

• \(\hat\epsilon_{(n)}\) quantifies the distance of \(Y_n\) from the regression line estimated after excluding the nth observation.
• \(\hat Y_{(n)}\) is the predicted value for the nth observatin from the regression line estimated after exluding the nth observation

We identify high-leverage observations through their leverage values.

• An observation’s leverage, \(h_n\), quantifies the extent to which its predictors affect the fitted regression model.
• Observations with \(X\) values very far from the mean, \(\bar X\), affect the fitted model disproportionately.

Index plots of the leverage values can help spotlight high-leverage points.

• Again, look for observations that clearly “stand out from the crowd.”

In simple linear regression, the nth leverage is given by: \[h_n=\frac 1 N + \frac{(X_n-\bar X)^2}{\sum^N_{m=1}(X_m-\bar X)^2}\]

hatvalues(out1) %>% plot()

Observations with high leverage or large (externally) studentized residuals are not necessarily influential.

• High-leverage observations tend to be more influential than outliers.
• The worst problems arise from observations that are both outliers and have high leverage.

Measures of influence simultaneously consider extremity in both X and Y dimensions.

• Observations with high measures of influence are very likely to cause problems.

All measures of influence use the same logic as the deleted residual.

• Compare models estimated from the whole sample to models estimated from samples excluding individual observations.

One of the most common measures of influence is Cook’s Distance \[Cook's D_n=\frac{\sum^N_{n=1}(\hat Y_n-\hat Y_{(n)})^2}{(P+1)\hat \sigma^2}=(P+1)^{-1}t^2_n\frac{h_n}{1-h_n}\]

Index plots of Cook’s distances can help spotlight the influential points.

• Look for observations that clearly “stand out from the crowd.”

cd <- cooks.distance(out1)
plot(cd)

Observation number 28 was the most influential according to Cook’s Distance

(maxD <- which.max(cd))

## Exclude the influential case:
Cars93.2 <- Cars93[-maxD, ]
## Fit model with reduced sample:
out2 <- lm(Price ~ Horsepower, data = Cars93.2)
round(summary(out1)$coefficients, 6)

round(summary(out2)$coefficients, 6)

Removing the outlier did not have a massive effect on the slope but it did lower the t-values by quite a bit

You could also remove the 2 most influential observations

(maxDs <- sort(cd) %>% names() %>% tail(2) %>% as.numeric())

## Exclude influential cases:
Cars93.2 <- Cars93[-maxDs, ]

## Fit model with reduced sample:
out2.2 <- lm(Price ~ Horsepower, data = Cars93.2)

The most common way to address influential observations is simply to delete them and refit the model.

• This approach is often effective—and always simple—but it is not fool-proof.
• Although an observation is influential, we may not be able to justify excluding it from the analysis.

Robust regression procedures can estimate the model directly in the presence of influential observations.

• Observations in the tails of the distribution are weighted less in the estimation process, so outliers and high-leverage points cannot exert substantial influence on the fit.
• We can do robust regression with the rlm() function from the MASS package.

In R you can train a generalized linear model with glm

data(iris)

## Generalized linear model:
glmFit <- glm(Petal.Length ~ Petal.Width + Species,
              family = gaussian(link = "identity"),
              data = iris)

Assumption 1 implies a linear relation between continuous predictors and the logit of the success probability.

• We can basically evaluate the linearity assumption using the same methods we applied with linear regression.

car::crPlots(glmFit, terms = ~ age + fare)

Basically the same as for linear regression Assumption 2 implies two conditions:

1. P > N
2. No severe (multi)collinearity among the predictors

We can quantify multicollinearity with the variance inflation factor (VIF).

car::vif(glmFit)

VIF > 10 indicates severe multicollinearity.

Assumption 3 implies several conditions.

1. The outcome, Y, is binary.
2. The linear predictor, 𝜂, can explain all the systematic trends in \(\pi\).
   • No residual clustering after accounting for X.
   • No important variables omitted from X.

We can easily check the first condition with summary statistics.

levels(titanic$survived)

table(titanic$survived)

If we have a non-binary, categorical outcome, we can use a different type of model.

• Multiclass nominal variables: multinomial logistic regression
  • nnet::multinom()
• Ordinal variables: Proportional odds logistic regression
  • MASS::polr()
• Counts: Poisson regression
  • glm() with family = "poisson"

The binomial distribution is also appropriate for modeling the proportion of successes in N trials.

For the second condition:

We can check for residual clustering by calculating the ICC using deviance residuals.

## Check for residual dependence induced by 'class':
ICC::ICCbare(x = titanic$class, y = resid(glmFit, type = "deviance"))

Logistic regression models are estimated with numerical methods, so we need larger samples than we would for linear regression models.

• The sample size requirements increase with model complexity.

Some suggested rules of thumb:

• 10 cases for each predictor (Agresti, 2018)
• \(N = 10P/\pi_0\) (Peduzzi, Concato, Kemper, Holford, & Feinstein, 1996)
  • \(P\): Number of predictors
  • \(\pi_0\): Proportion of the minority class
• \(N = 100 + 50P\) (Bujang, Omar, & Baharum, 2018)

The logistic regression may not perform well when the outcome classes are severely imbalanced.

• e.g. way more of one class then the other

with(titanic, table(survived) / length(survived))

We have a few possible solutions for problematic imbalance:

• Down-sampling the majority class
  • Randomly remove some samples from the majority class
• Up-sampling the minority class
  • Random select samples from the minority class to use multiple times
• Use weights when estimating the logistic regression model
  • weights argument in glm()

We don’t actually want to perfectly predict class membership.

• The model cannot estimate with perfectly separable classes.

Why?

• Imagine the shape of the logit function to estimate classes which can be seperated perfectly
• The slope would have to be really low, then really high
• To get this the slope would have to be infinite, which obviously can’t be estimated

Model regularization (e.g., ridge or LASSO penalty) may help.

• glmnet::glmnet()

One of the most direct ways to evaluate classification performance is the confusion matrix.

## Add predictions to the dataset:
titanic %<>%
    mutate(piHat = predict(glmFit, type = "response"),
           yHat = as.factor(ifelse(piHat <= 0.5, "no", "yes")))

library(caret)

cMat <- titanic %$% confusionMatrix(data = yHat, reference = survived)

cMat$table

Each cell in the confusion matrix represents a certain classification result:

• TP true positive: Correctly predicted positive
• TN true negative: Correctly predicted negative
• FP false positive: Falsely predicted positive
• FN false negative: Falsely predicted negative

cMat$overall

cMat$byClass

Summaries/statistics of the Confusion matrix:

• \(Accuracy = (TP+TN)/(P+N)\)
• \(Error\ Rate = (FP+FN)/(P+N)=1-Accuracy\)
• \(Sensitivity = TP/(TP+FN)\)
• \(Specificity = TN/(TN+FP)\)
• \(False\ Positive\ Rate (FPR)=FP/(TN+FP)=1-Specificity\)
• \(Positive\ Predictive\ Value\ (PPV) = TP / (TP + FP)\)
• \(Negative\ Predictive\ Value\ (NPV) = TN / (TN + FN)\)

A receiver operating characteristic (ROC) curve illustrates the diagnostic ability of a binary classifier for all possible values of the classification threshold.

• The ROC curve plots sensitivity against specificity at different threshold values.

The area under the ROC curve (AUC) is a one-number summary of the potential performance of the classifier.

• The AUC does not depend on the classification threshold.

pROC::auc(rocData)

Rule of thumb for AUC according to Mandrekar (2010):

• AUC value from 0.7 – 0.8: Acceptable
• AUC value from 0.8 – 0.9: Excellent
• AUC value over 0.9: Outstanding

We can use numerical methods to estimate an optimal threshold value

library(OptimalCutpoints)
ocOut <- optimal.cutpoints(X = "piHat",
                           status = "survived",
                           tag.healthy = "no",
                           data = titanic,
                           method = "ROC01"
                           )

partSummary(ocOut, -1)

Area under the ROC curve (AUC): 0.83 (0.802, 0.858)

CRITERION: ROC01

Number of optimal cutoffs: 1

Not exam material, might be usefull for the practical though.

Measuring classification performance from a confusion matrix can be problematic.

• Sometimes too coarse.

We can also base our error measure on the residual deviance with the Cross-Entropy Error: \[CEE = -N^{-1} \sum^N_{n=1} Y_n\ln(\hat\pi_n) + (1 - Y_n) \ln(1 -\hat\pi_n)\]

• The CEE is sensitive to classification confidence.
• Stronger predictions are more heavily weighted

The misclassification rate is a na¨ıvely appealing option.

• The proportion of cases assigned to the wrong group

Consider two perfect classifiers:

1. P( ˆYn = 1|Yn = 1) = 0.90, P( ˆYn = 1|Yn = 0) = 0.10, n = 1, 2, . . . , N
2. P( ˆYn = 1|Yn = 1) = 0.55, P( ˆYn = 1|Yn = 0) = 0.45, n = 1, 2, . . . , N

Both of these classifiers will have the same misclassification rate.

• Neither model ever makes an incorrect group assignment.

The first model will have a lower CEE.

• The classifications are made with higher confidence.
• CEE1 = 0.105, CEE2 = 0.598