2  A logistic regression model for predicting survival on the titanic

Can we predict which passengers would survive the sinking of the Titanic?

See:

https://www.kaggle.com/c/titanic/overview/evaluation

https://michaelallen1966.github.io/titanic/front_page.html

The original data comes from:

https://www.kaggle.com/c/titanic/data

Though we will download and use a data set that has been pre-processed ready for machine learning. IMPORTANT - real data that you will work with will need preprocessing first. There’s a handy guide to preprocessing data here (based on the Titanic data we’re using) : https://michaelallen1966.github.io/titanic/01_preprocessing.html

The data includes:

Variable	Definition
survival	Survival (0 = No, 1 = Yes)
pclass	Ticket class
sex	Sex
Age	Age in years
sibsp	# of siblings / spouses aboard the Titanic
parch	# of parents / children aboard the Titanic
ticket	Ticket number
fare	Passenger fare
cabin	Cabin number
embarked	Port of Embarkation(C=Cherbourg, Q=Queenstown, S=Southampton)

2.1 Logistic regression

In this example we will use logistic regression (see https://en.wikipedia.org/wiki/Logistic_regression).

For an introductory video on logistic regression see: https://www.youtube.com/watch?v=yIYKR4sgzI8

Logistic regression takes a range of features (which we will normalise/standardise to put on the same scale) and returns a probability that a certain classification (survival in this case) is true.

We will go through the following steps:

• Download and save pre-processed data
• Split data into features (X) and label (y)
• Split data into training and test sets (we will test on data that has not been used to fit the model)
• Standardise data
• Fit a logistic regression model (from sklearn)
• Predict survival of the test set, and assess accuracy
• Review model coefficients (weights) to see importance of features
• Show probability of survival for passengers

2.2 Load modules

A standard Anaconda install of Python (https://www.anaconda.com/distribution/) contains all the necessary modules. Use your base environment if in doubt.

import numpy as np
import pandas as pd
# Import machine learning methods
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

2.3 Load data

The section below downloads pre-processed data, and saves it to a subfolder (from where this code is run). If data has already been downloaded that cell may be skipped.

Code that was used to pre-process the data ready for machine learning may be found at: https://github.com/MichaelAllen1966/1804_python_healthcare/blob/master/titanic/01_preprocessing.ipynb

download_required = True

if download_required:

    # Download processed data:
    address = 'https://raw.githubusercontent.com/MichaelAllen1966/' + \
                '1804_python_healthcare/master/titanic/data/processed_data.csv'

    data = pd.read_csv(address)

    # Create a data subfolder if one does not already exist
    import os
    data_directory ='./data/'
    if not os.path.exists(data_directory):
        os.makedirs(data_directory)

    # Save data
    data.to_csv(data_directory + 'processed_data.csv', index=False)

data = pd.read_csv('data/processed_data.csv')
# Make all data 'float' type
data = data.astype(float)

2.4 Examine loaded data

The data is in the form of a Pandas DataFrame, so we have column headers providing information of what is contained in each column.

We will use the DataFrame .head() method to show the first few rows of the imported DataFrame. By default this shows the first 5 rows. Here we will look at the first 10.

data.head(10)

	PassengerId	Survived	Pclass	Age	SibSp	Parch	Fare	AgeImputed	EmbarkedImputed	CabinLetterImputed	...	Embarked_missing	CabinLetter_A	CabinLetter_B	CabinLetter_C	CabinLetter_D	CabinLetter_E	CabinLetter_F	CabinLetter_G	CabinLetter_T	CabinLetter_missing
0	1.0	0.0	3.0	22.0	1.0	0.0	7.2500	0.0	0.0	1.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	1.0
1	2.0	1.0	1.0	38.0	1.0	0.0	71.2833	0.0	0.0	0.0	...	0.0	0.0	0.0	1.0	0.0	0.0	0.0	0.0	0.0	0.0
2	3.0	1.0	3.0	26.0	0.0	0.0	7.9250	0.0	0.0	1.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	1.0
3	4.0	1.0	1.0	35.0	1.0	0.0	53.1000	0.0	0.0	0.0	...	0.0	0.0	0.0	1.0	0.0	0.0	0.0	0.0	0.0	0.0
4	5.0	0.0	3.0	35.0	0.0	0.0	8.0500	0.0	0.0	1.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	1.0
5	6.0	0.0	3.0	28.0	0.0	0.0	8.4583	1.0	0.0	1.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	1.0
6	7.0	0.0	1.0	54.0	0.0	0.0	51.8625	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	1.0	0.0	0.0	0.0	0.0
7	8.0	0.0	3.0	2.0	3.0	1.0	21.0750	0.0	0.0	1.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	1.0
8	9.0	1.0	3.0	27.0	0.0	2.0	11.1333	0.0	0.0	1.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	1.0
9	10.0	1.0	2.0	14.0	1.0	0.0	30.0708	0.0	0.0	1.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	1.0

10 rows × 26 columns

Note that in the above, “Imputed” column names indicate whether that feature was blank in the original data and was filled in (“imputed”) in the preprocessing (so “AgeImputed” value of 1.0 means the age was missing, 0.0 means it wasn’t). You can’t have missing cells in Machine Learning, so you have to decide what to do when you have missing data. This is part of the pre-processing step you’ll need to do with real data, and there’s more information about this in the pre-processing notebook here : https://michaelallen1966.github.io/titanic/01_preprocessing.html

We can also show a summary of the data with the .describe() method.

data.describe()

	PassengerId	Survived	Pclass	Age	SibSp	Parch	Fare	AgeImputed	EmbarkedImputed	CabinLetterImputed	...	Embarked_missing	CabinLetter_A	CabinLetter_B	CabinLetter_C	CabinLetter_D	CabinLetter_E	CabinLetter_F	CabinLetter_G	CabinLetter_T	CabinLetter_missing
count	891.000000	891.000000	891.000000	891.000000	891.000000	891.000000	891.000000	891.000000	891.000000	891.000000	...	891.000000	891.000000	891.000000	891.000000	891.000000	891.000000	891.000000	891.000000	891.000000	891.000000
mean	446.000000	0.383838	2.308642	29.361582	0.523008	0.381594	32.204208	0.198653	0.002245	0.771044	...	0.002245	0.016835	0.052750	0.066218	0.037037	0.035915	0.014590	0.004489	0.001122	0.771044
std	257.353842	0.486592	0.836071	13.019697	1.102743	0.806057	49.693429	0.399210	0.047351	0.420397	...	0.047351	0.128725	0.223659	0.248802	0.188959	0.186182	0.119973	0.066890	0.033501	0.420397
min	1.000000	0.000000	1.000000	0.420000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	...	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
25%	223.500000	0.000000	2.000000	22.000000	0.000000	0.000000	7.910400	0.000000	0.000000	1.000000	...	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	1.000000
50%	446.000000	0.000000	3.000000	28.000000	0.000000	0.000000	14.454200	0.000000	0.000000	1.000000	...	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	1.000000
75%	668.500000	1.000000	3.000000	35.000000	1.000000	0.000000	31.000000	0.000000	0.000000	1.000000	...	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	1.000000
max	891.000000	1.000000	3.000000	80.000000	8.000000	6.000000	512.329200	1.000000	1.000000	1.000000	...	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000

8 rows × 26 columns

The first column is a passenger index number. We will remove this, as this is not part of the original Titanic passenger data, and will not help us train our model.

# Drop Passengerid (axis=1 indicates we are removing a column rather than a row)
# We drop passenger ID as it is not original data
# inplace=True means change the dataframe itself - don't create a copy with this column dropped

data.drop('PassengerId', inplace=True, axis=1)

2.5 Looking at a summary of passengers who survived or did not survive

Before running machine learning models, it is always good to have a look at your data. Here we will separate passengers who survived from those who died, and we will have a look at differences in features.

We will use a mask to select and filter passengers. The mask applies Boolean values (True and False) to entries depending on a given condition. Below, we use this to create a mask that has True values for any rows where the Survived value is 1.0 (ie where the patient survived), and then store only those rows in a separate dataframe called “survived”. Then we do the same thing but for those who died.

mask = data['Survived'] == 1 # Mask for passengers who survive
survived = data[mask] # filter using mask

mask = data['Survived'] == 0 # Mask for passengers who died
died = data[mask] # filter using mask

Now let’s look at average (mean) values for each feature for those who survived and those who died. We can make comparing them easier by putting these values in a new DataFrame so we can look at them side by side. What do you notice? What features do you think might have influenced survival?

summary = pd.DataFrame() # New empty DataFrame
summary['survived'] = survived.mean()
summary['died'] = died.mean()

summary

	survived	died
Survived	1.000000	0.000000
Pclass	1.950292	2.531876
Age	28.291433	30.028233
SibSp	0.473684	0.553734
Parch	0.464912	0.329690
Fare	48.395408	22.117887
AgeImputed	0.152047	0.227687
EmbarkedImputed	0.005848	0.000000
CabinLetterImputed	0.602339	0.876138
CabinNumber	18.961988	6.074681
CabinNumberImputed	0.611111	0.885246
male	0.318713	0.852459
Embarked_C	0.271930	0.136612
Embarked_Q	0.087719	0.085610
Embarked_S	0.634503	0.777778
Embarked_missing	0.005848	0.000000
CabinLetter_A	0.020468	0.014572
CabinLetter_B	0.102339	0.021858
CabinLetter_C	0.102339	0.043716
CabinLetter_D	0.073099	0.014572
CabinLetter_E	0.070175	0.014572
CabinLetter_F	0.023392	0.009107
CabinLetter_G	0.005848	0.003643
CabinLetter_T	0.000000	0.001821
CabinLetter_missing	0.602339	0.876138

2.6 Divide into X (features) and y (labels)

We will separate out our features (the data we use to make a prediction) from our label (what we are trying to predict - survival here). By convention our features are called X (usually upper case to denote multiple features), and the label (survived or not) y.

X = data.drop('Survived',axis=1) # X = all 'data' except the 'survived' column
y = data['Survived'] # y = 'survived' column from 'data'

2.7 Divide into training and tets sets

When we test a machine learning model we should always test it on data that has not been used to train the model. We will use sklearn’s train_test_split method to randomly split the data: 75% for training, and 25% for testing.

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.25)

Let’s have a look at the standard deviation and the mean of the feature values in the training data (Standard Deviation for each feature shown first, then mean for each feature)

X_train.std(), X_train.mean()

(Pclass                  0.831431
 Age                    12.963059
 SibSp                   1.147940
 Parch                   0.803459
 Fare                   46.558286
 AgeImputed              0.411525
 EmbarkedImputed         0.054677
 CabinLetterImputed      0.421493
 CabinNumber            25.933916
 CabinNumberImputed      0.413581
 male                    0.477421
 Embarked_C              0.377799
 Embarked_Q              0.288276
 Embarked_S              0.442443
 Embarked_missing        0.054677
 CabinLetter_A           0.143351
 CabinLetter_B           0.219955
 CabinLetter_C           0.240201
 CabinLetter_D           0.186250
 CabinLetter_E           0.197088
 CabinLetter_F           0.121524
 CabinLetter_G           0.066915
 CabinLetter_T           0.038691
 CabinLetter_missing     0.421493
 dtype: float64,
 Pclass                  2.318862
 Age                    29.352171
 SibSp                   0.555389
 Parch                   0.387725
 Fare                   31.497535
 AgeImputed              0.215569
 EmbarkedImputed         0.002994
 CabinLetterImputed      0.769461
 CabinNumber            10.673653
 CabinNumberImputed      0.781437
 male                    0.649701
 Embarked_C              0.172156
 Embarked_Q              0.091317
 Embarked_S              0.733533
 Embarked_missing        0.002994
 CabinLetter_A           0.020958
 CabinLetter_B           0.050898
 CabinLetter_C           0.061377
 CabinLetter_D           0.035928
 CabinLetter_E           0.040419
 CabinLetter_F           0.014970
 CabinLetter_G           0.004491
 CabinLetter_T           0.001497
 CabinLetter_missing     0.769461
 dtype: float64)

2.8 Standardise data

We can see above that there are quite different scales across different features in this data. For example, passenger class is on a very different scale numerically than the fare paid.

We want all of our features to be on roughly the same scale. This generally leads to a better model, and also allows us to more easily compare the importance of different features.

We will use standardisation to scale our feature values, where we use the mean and standard deviation of the training set of data to normalise the data. We subtract the mean of the training set values, and divide by the standard deviation of the training data. Note that the mean and standard deviation of the training data are used to standardise the test set data as well.

Here we will use sklearn’s StandardScaler method. This method also copes with problems we might otherwise have (such as if one feature has zero standard deviation in the training set). We write a little function so whenever we need to standardise some data, we can just call this function, pass in the training and test feature (X) data, and it’ll return the same data but standardised.

def standardise_data(X_train, X_test):

    # Initialise a new scaling object for normalising input data
    sc = StandardScaler()

    # Apply the scaler to the training and test sets
    train_std=sc.fit_transform(X_train)
    test_std=sc.fit_transform(X_test)

    return train_std, test_std

Now let’s call this function and use it to standardise our training and test data.

X_train_std, X_test_std = standardise_data(X_train, X_test)

2.9 Fit logistic regression model

Now we will fit a logistic regression model, using sklearn’s LogisticRegression method. Our machine learning model fitting (training) is only two lines of code! By using the name model for our logistic regression model we will make our model more interchangeable later on if we wanted to try a different kind of model, for example.

model = LogisticRegression()
model.fit(X_train_std,y_train)

LogisticRegression()

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

2.10 Predict values

Now we can use the trained model to predict survival. We will test the accuracy of both the training and test data sets.

# Predict training and test set labels
y_pred_train = model.predict(X_train_std)
y_pred_test = model.predict(X_test_std)

2.11 Calculate accuracy

In this example we will measure accuracy simply as the proportion of passengers where we make the correct prediction. In later examples we will look at other measures of accuracy which explore false positives and false negatives in more detail.

# The shorthand below says to check each predicted y value against the actual
# y value in the training data.  This gives a list of True and False values
# for each prediction, where True indicates the predicted value matches the
# actual value.  Then we take the mean of these Boolean values, which gives
# us a proportion (where if all values were True, the proportion would be 1.0)
# If you want to see why that works, just uncomment the following line of code
# to see what y_pred_train == y_train is doing.
# print (y_pred_train == y_train)
accuracy_train = np.mean(y_pred_train == y_train)
accuracy_test = np.mean(y_pred_test == y_test)

print (f'Accuracy of predicting training data = {accuracy_train}')
print (f'Accuracy of predicting test data = {accuracy_test}')

Accuracy of predicting training data = 0.8203592814371258
Accuracy of predicting test data = 0.7982062780269058

Not bad - about 80% accuracy. You will probably see that accuracy of predicting the training set is usually higher than the test set. Because we are only testing one random sample, you may occasionally see otherwise. In later examples we will look at the best way to repeat multiple tests, and look at what to do if the accuracy of the training set is significantly higher than the test set (a problem called ‘overfitting’).

2.12 Examining the model coefficients (weights)

Not all features are equally important. And some may be of little or no use at all, unnecessarily increasing the complexity of the model. In later examples we will look at selecting features which add value to the model (or removing features that don’t).

Here we will look at the importance of features – how they affect our estimation of survival. These are known as the model coefficients (if you come from a traditional statistics background), or model weights (if you come from a machine learning background).

Because we have standardised our input data the magnitude of the weights may be compared as an indicator of their influence in the model. Weights with higher negative numbers mean that that feature correlates with reduced chance of survival. Weights with higher positive numbers mean that that feature correlates with increased chance of survival. Those weights with values closer to zero (either positive or negative) have less influence in the model.

We access the model weights my examining the model coef_ attribute. The model may predict more than one outcome label, in which case we have weights for each label. Because we are predicting a single label here (survive or not), the weights are found in the first element ([0]) of the coef_ attribute.

co_eff = model.coef_[0]
co_eff

array([-0.84475058, -0.55123105, -0.40925875, -0.01028495,  0.01404105,
       -0.21700132,  0.09049331,  0.01086364, -0.04374518, -0.38007729,
       -1.27635266,  0.09730836,  0.14421872, -0.18824044,  0.09049331,
       -0.04918584,  0.04665105, -0.13362438,  0.0071312 ,  0.13648839,
        0.10288958, -0.14729349, -0.16976434,  0.01086364])

So we have an array of model weights.

Not very readable for us mere humans is it?!

We will transfer the weights array to a Pandas DataFrame. The array order is in the same order of the list of features of X, so we will put that those into the DataFrame as well. And we will sort by influence in the model. Because both large negative and positive values are more influential in the model we will take the absolute value of the weight (ie remove any negative sign), and then sort by that absolute value. That will give us a more readable table of most influential features in the model.

co_eff_df = pd.DataFrame() # create empty DataFrame
co_eff_df['feature'] = list(X) # Get feature names from X
co_eff_df['co_eff'] = co_eff
co_eff_df['abs_co_eff'] = np.abs(co_eff)
co_eff_df.sort_values(by='abs_co_eff', ascending=False, inplace=True)

Let’s look at the DataFrame. What do you conclude?

co_eff_df

	feature	co_eff	abs_co_eff
10	male	-1.276353	1.276353
0	Pclass	-0.844751	0.844751
1	Age	-0.551231	0.551231
2	SibSp	-0.409259	0.409259
9	CabinNumberImputed	-0.380077	0.380077
5	AgeImputed	-0.217001	0.217001
13	Embarked_S	-0.188240	0.188240
22	CabinLetter_T	-0.169764	0.169764
21	CabinLetter_G	-0.147293	0.147293
12	Embarked_Q	0.144219	0.144219
19	CabinLetter_E	0.136488	0.136488
17	CabinLetter_C	-0.133624	0.133624
20	CabinLetter_F	0.102890	0.102890
11	Embarked_C	0.097308	0.097308
6	EmbarkedImputed	0.090493	0.090493
14	Embarked_missing	0.090493	0.090493
15	CabinLetter_A	-0.049186	0.049186
16	CabinLetter_B	0.046651	0.046651
8	CabinNumber	-0.043745	0.043745
4	Fare	0.014041	0.014041
7	CabinLetterImputed	0.010864	0.010864
23	CabinLetter_missing	0.010864	0.010864
3	Parch	-0.010285	0.010285
18	CabinLetter_D	0.007131	0.007131

Now - the actual order of features, and the weights for each, will be different each time you fit the model, as it’s random and we haven’t put in place any controls. So don’t worry if yours differs a bit. But, you’ll likely find that amongst the top weighted features are :

• male (being male reduces probability of survival)
• Pclass (lower class passengers, who have a higher class number, reduces probability of survival)
• age (being older reduces probability of survival)
• CabinNumberImputed (cabin number is missing, which may mean they didn’t have a cabin - likely lower class)

2.13 Show predicted probabilities

The predicted probabilities are for the two alternative classes 0 (does not survive) or 1 (survive).

Ordinarily we do not see these probabilities - the predict method used above applies a cut-off of 0.5 to classify passengers into survived or not, but we can see the individual probabilities for each passenger if desired.

In a later example we will see how we can use such probabilities to adjust the sensitivity of our model to detecting survivors or non-survivors.

Each passenger has two values. These are the probability of not surviving (first value) or surviving (second value). Because we only have two possible classes we only need to look at one. Multiple values are important when there are more than one class being predicted.

# Show first ten predicted classes
classes = model.predict(X_test_std)
classes[0:10]

array([0., 0., 0., 0., 1., 0., 0., 0., 0., 1.])

# Show first ten predicted probabilities
# (note how the values relate to the classes predicted above)
probabilities = model.predict_proba(X_test_std)
probabilities[0:10]

array([[0.9585624 , 0.0414376 ],
       [0.98719408, 0.01280592],
       [0.83298478, 0.16701522],
       [0.94523796, 0.05476204],
       [0.04687281, 0.95312719],
       [0.5156704 , 0.4843296 ],
       [0.91270338, 0.08729662],
       [0.89453433, 0.10546567],
       [0.5000722 , 0.4999278 ],
       [0.08915059, 0.91084941]])