Structure, Preview and Statistics

Three methods give initial insights about the data:

• str shows the structure of the data;
• head provides a preview the first rows of the data set;
• summary gives a list of basic descriptive statistics.

How many observations and variables are in the “faithful” data set?

str (faithful)

Answer: 272 observations and 2 variables

Display the first six observations?

head   (faithful)

Answer: 3.6 79, …

What is the shortest, expected and third quartile eruption duration?

summary(faithful)

Answer: Min 1.6min, mean: 3.488 and 3rd Quartile is 4.454min.

Visual inspection

Display the data as a scatter plot and add a “trend line”.

# library(ggplot2)
g <- ggplot(data = faithful,aes(waiting,eruptions))+
   geom_point()+stat_smooth(method="lm")
g 

The data seems to have some linear relationship. Although it looks like there are two classes, one with long durations (above 3min) and another one with short durations (below 3min) - more about this in the logistic regression section.

Let us have a bit of fun: plot the means, min and max, etc. and add some colour and size. Let us also add the unit to the x and y axes.

# library(ggplot2)
attach(faithful)
me = mean(eruptions); sde=sd(eruptions);
g +  geom_hline(yintercept=me) +
     geom_text(aes(x=50, y=me, label=sprintf('y=%.2f',me), vjust=-0.4))+
     geom_vline(xintercept=mean(waiting)) +
     geom_hline(yintercept=min(eruptions), color="red", size=2)+ 
     geom_hline(yintercept=max(eruptions), color="blue", size=2)+ 
     geom_hline(yintercept=me+sde, color="red", linetype=4)+
     geom_hline(yintercept=me-sde, color="red", linetype=4)+
     xlab('x = waiting time [minutes]')+
     ylab('y = eruptions [minutes]')

Correlation

We can also see that the data is correlated: cor(faithful) gives us 0.901. That means the Pearson Correlation Coefficient (PCC) \(r\) reveals a strong association between the variables. The PCC is defined as: \[ r = r_{X,Y} = \frac{s_{XY}}{s_X s_Y},\] where \(s_X\) is the sample standard deviation of \(X\) and \(s_{XY} = \text{cov}(X,Y)\) is the sample covariance. The sample standard deviation \(s\) is obtained form the sample variance \(s^2\) by taking the square root and is defined as: \[s^2 = \frac{1}{n-1} \sum_{i=1}^n (x_i-\bar{x})^2.\] The sample covariance is defined as: \[ \text{cov}(x,y) = \frac{\sum_i(x_i-\bar{x})(y_i-\bar{y})}{n-1}. \]

The above can be determined with the following functions in : correlation cor(), standard deviation sd, variance var() and covariance cov().

Verify that the functions are correct. You can use the attach function for easier access of a data.frame.

attach(faithful)
x = eruptions; y = waiting; n = nrow(faithful);
cat('cor(x,y) = ',  cor(x,y),
    '=cov(x,y)/(sd(x)*sd(y)) =', cov(x,y)/(sd(x)*sd(y)))
cat('cov(x,y) = ',  cov(x,y),
    ' = ', sum( (x - mean(x)) * (y - mean(y)) )/(n-1))

The coefficient of determination is defined as: \[ R^2 = \frac{\sum_{i=1}^n (\hat{y}_i - \bar{y})^2}{\sum_{i=1}^n (y_i - \bar{y})^2}\] That means, that even without the prediction model we are already in the position to give an accuracy statement about the goodness of the model. The coefficient of determination is 0.81.

Using the Model

In the previous section we determined the models coefficient’s and saved them as vector c.

What is the eruption duration given that we waited \(w=80\) minutes? Hint: Use the coefficients from before.

w = 80 #min waiting time
d = c[1] + c[2]*w
d

Solution: We expect to observe an eruption lasting for 4.18min (given 80min waiting time).

As you can see in the result the name of c[1] was used to name the variable d. You can change this using names().

names(d) <- 'duration'
d

Multiple predictions: Predict the eruption times for times 40, 60, 80 and 100 minutes. Can you observe any systematic change?

p = data.frame(waiting=c(40, 60,80,100))
predict(model,p)

Answer: every 20 minutes of waiting gives us an additional 1.51 minutes of eruption (see waiting coefficient \(\beta_1\)).

Math Approach(optional)

Here, we will show to approaches to determine the model analytically. Note, this section is optional. It assumes some familiarity with linear algebra (vectors, matrices, calculus).

Our objective is to minimise the error vector \(\epsilon\) between the vectors \(y\) (observed values) and \(\hat{y}\) (modelled values), because \(y = \hat{y} + \epsilon\). If you look at a 2D-vector it is obvious that the error is at a minimum, when \(\hat{y}\) and \(\epsilon\) are orthogonal, i.e. \(\hat{y} \bot \epsilon\). Two vectors are orthogonal when their product is zero: \[ \hat{y} \bot \epsilon \Rightarrow\hat{y}^\intercal \epsilon=0.\]

The modelled (predicted) values are determined by \(\hat{y}=Xb\). The error is \(\epsilon = y - \hat{y} = y - Xb\). Using the above two formulas allows us to substitute \(\hat{y}\) and \(\epsilon\): \[ \hat{y}^\intercal \epsilon=0 \Rightarrow(Xb)^\intercal (y-Xb)=0.\] Applying the transposition, associative law, distributive law, and realising that we can drop \(b^\intercal\) gives us: \[ (Xb)^\intercal (y-Xb) = b^\intercal X^\intercal (y-Xb) = 0 \Rightarrow X^\intercal y = (X^\intercal X)b \] Since, \(X^\intercal X\) is a squared matrix we can compute the inverse. \[ X^\intercal y = (X^\intercal X)b \Rightarrow(X^\intercal X)^{-1}X^\intercal y = (X^\intercal X)^{-1}(X^\intercal X)b\] Hence, the prediction model coefficients are: \[ b = (X^\intercal X)^{-1}X^\intercal y \]

Now let us use the above formula to derive the coefficients for the geyser prediction model. The interesting aspect is that the matrix inverse is obtained using solve rather than ^(-1).

y = faithful$eruptions
x = faithful$waiting
e = rep(1,nrow(faithful)) # one vector
X = cbind(e,x)

b = solve((t(X)%*%X)) %*% t(X) %*% y
b

Note that sometimes there are duplicated records, which may exist on purpose (e.g. same scenario observed multiple times). In this case the pseudo-inverse (generalised-inverse) needs to be used, e.g. the Moore-Penrose inverse.

library(MASS) # obtain ginv function

bp = ginv((t(X)%*%X)) %*% t(X) %*% y
bp

Training and Test Data

It is common to divide the data into training and test data. This gives reassurance that the fitted model will work for unseen data (simulated by the test data). The training data is used to train and learn the model. The test data ‘’validates’’ the learned model.

Let us do a 70:30 data split, i.e. 70% of the data will be used for learning the model randomly; and the remaining 30% of the data will be used to test the accuracy of the model. We will make the split reproducible by fixing the seed of the random number with set.seed(1).

n = nrow(faithful)  # number of observations
ns = floor(.7*n)    # 70% times n = number of training data records
set.seed(1)         # seed of random generator (for reproducibility)
I = sample(1:n, ns) # sample of indexes
train = faithful[I,]    # training data (=70% of data)
test  = faithful[-I,]   # test data (=30% of data)

Use the training data to learn the linear regression model, and provide the RMSE for the training data.

lr <- lm (eruptions ~ waiting, data=train)

yh = predict(lr)
y  = train$eruptions
rmse  <- function(y,yh) {sqrt(mean((y-yh)^2))}
rmse(y,yh)

Now compute the test RMSE for the model to see how well it performs for unseen data.

yh = predict(lr, newdata = test)
y  = test$eruptions
rmse(y,yh)

You should observe that the RMSE for the test data is higher than the RMSE for the training data.

Note a different seed will lead to a different accuracy. Hence, it would be better to repeat the splitting multiple times and use the average accuracy instead.

Practice makes perfect

Use the code written previously to specify the training and test data by embedding it in a function and show its usage.

dataSplit <- function(data=faithful, p=.7, seed=1){
  n = nrow(data) # number of observations
  ns = floor(p*n) # 70% times n = number of training data records
  set.seed(seed) # seed of random generator (for reproducibility)
  I = sample(1:n, ns) # sample of indexes
  train = data[I,] # training data (=70% of data)
  test = data[-I,] # test data (=30% of data)
  D = list(train=train, test=test)
  return(D)
}

D = dataSplit()
str(D)
# attach(D) # to have access to train and test directly

Analyse the intercountry life-cycle savings data.

• Provide basic descriptive statistics (what is the expected personal saving?);
• Analyse the correlations (which two features have the highest correlation?);
• What is the root mean squared error of the LR model?
• What are the aggregated personal savings for pop15=35, pop75=2.3, dpi=800, ddpi=3.5? (use linear regression)

summary(LifeCycleSavings)
r = cor(LifeCycleSavings)
r[r==1]=0; max(r)

lr <- lm(sr ~ ., data = LifeCycleSavings)

yh = predict(lr)
y  = LifeCycleSavings$sr
sqrt(mean(sum((y-yh)^2))) # RMSE

df <- data.frame(pop15=35, pop75=2.3, dpi=800, ddpi=3.5)
predict(lr, newdata=df)

Answer: (a) mean sr = 10.51, (b) pop75 and dpi, (c) RMSE = 25.1, (d) sr = 9.70

Data source

The data is obtained using the Titanic library. This is based on the Kaggle data set. Check the content of the titanic library. Does it contain a training and test data set?

# library(titanic)
ls("package:titanic")

Please inspect the training and test data set. Does the test data set contain passenger information?

#?titanic_train
?titanic_test

Answer: no passenger information is available, i.e. we cannot determine the test accuracy of the model. Note: the information was not provided, since it is a competition on Kaggle.

By default there is another data set about the Titanic in the R environment. However, it contains only aggregated information.

There is a Titanic website, which contains even more data about the disaster.

However, we will use the previous train data set, and we will assume that this is the entire data available.

Explore data

Similar, to the linear regression explore the data using names, head and summary.

1. Does the data contain a Survived and Sex field?
2. What is the name in the sixth record?
3. How many records does the data set contain?
4. How many elements are not available in the Age field?

T = titanic_train # assume from now on this is training and test data.
names(T)
head(T)
summary(T)

Answers: 1) yes; 2) James Moran; 3) 891; 4) 177.

Data cleaning

The field Age has several “not available(na)” entries. Your task is to set those to the average age. What is the average age?

# sum(is.na(T$Age))
mv = mean(T$Age, na.rm = TRUE)
T$Age[is.na(T$Age)] = mv
mv

Answer: 29.7 years.

Create a training set D that contains 80% of the data (round down). How many records has D?

n = length(T$PassengerId) # number of passengers
ns = floor(0.8*n) # #records for training data set 

set.seed(1) # fix random number generator
I = sample(1:n, ns) # create sample indizis (no repetitions)

I[1:10] # show the first 10 records
length(I); # or a bit more precise: length(unique(I)) 
D = T[I,]

Answer: 712.

Two models

Let’s derive two simple prediction models.

# library(dplyr)
aD <- D %>% group_by(Survived) %>% summarise(nb = length(PassengerId)) 
aD %>% mutate(pnb = nb/ns)

test = T$Survived[-I]
mean(0 == test) 

aD2 <- D %>% group_by(Survived, Sex) %>% summarise(nb = length(PassengerId)) 
# aD2
# library (tidyr) # reshape aD2
aD2 %>% pivot_wider(Sex, names_from = Survived, values_from=nb)

pm <- rep(0,length(test)) # no one survives (as model 1)
pm[T$Sex[-I] =="female"]=1 # all females survive
mean(pm == test) 

Logistic Regression Model

Let’s derive a logistic regression model. Previously, we found the existing field names.

names(D)

From the previous models we have learned that Sex is an important feature. Let us use this field to build the logistic regression model f.

f <- glm(Survived ~ Sex, data = D, family = binomial ) # would result in the same result as above
#f <- glm(Survived ~ Sex + Age + Pclass+ SibSp + Fare + Embarked, data = D, family = binomial )
#summary(f)
summary(f)$coef

f_prob = predict(f, type = "response")
f_prob[1:10]

yh = rep(0, length(f_prob))
yh[f_prob>=0.5] = 1
yh[1:5]

C = table(yh, D$Survived)
C

# correctly predicted in training set
C[1,1]+C[2,2] # true negatives + true positives
# cat( 'ns = ', ns, ' = ', length(D$Survived), ' = ', sum(C), '\n')
(C[1,1]+C[2,2]) / ns
# mean(yh == D$Survived) # same as before 

# correctly predicted in test set
ft_prob = predict(f, newdata= T[-I,], type = "response") # note newdata
# length(ft_prob)
# ft_prob[1:10]

yht = rep(0, length(ft_prob)) # init
yht[ft_prob>=0.5] = 1 # set modeled response

Ct = table(yht, T[-I, "Survived"])
# Ct
mean(yht == T$Survived[-I])

fe <- glm(Survived ~ Sex + Age + Pclass+ SibSp + Fare + Embarked, data = D, family = binomial )
fe

fte_prob = predict(fe, newdata= T[-I,], type = "response") # test data probs.

yhte = rep(0, length(fte_prob)) # init
yhte[fte_prob>=0.5] = 1 # set response

mean(yhte == T$Survived[-I])