-
Notifications
You must be signed in to change notification settings - Fork 0
/
Lab4.R
523 lines (366 loc) · 15.9 KB
/
Lab4.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
####################################################
## ##
## 204B Lab 4 ##
## LOGISTIC REGRESSION AND GLM ##
## ##
## Kristine O'Laughlin ##
## Winter 2018 ##
## ##
####################################################
# Install and load the 'titanic' package into your R
# library using the following code:
#install.packages('titanic')
library(titanic)
# You will also need to combine the test and training
# datasets in order to get the complete dataset for
# this example:
data('titanic_train')
data('titanic_test')
titanic_test$Survived = NA
complete_data = rbind(titanic_train, titanic_test)
# And then subsetting the variables we will be using
# for this homework assignment:
titanic_dat = complete_data[, c('Survived', 'Sex', 'Age',
'Fare')]
# You will use this dataset to answer the following
# questions.
head(titanic_dat)
# Question 1
############
# Report the summary statistics for the titanic_dat
# dataset (i.e., report means, standard deviations,
# number of missing values). (Hint: this may be
# important for a later question on this homework!)
# Get descriptives for continuous variables:
means = sapply(titanic_dat[, 3:4], mean, na.rm = T)
sds = sapply(titanic_dat[, 3:4], sd, na.rm = T)
nas = colSums(sapply(titanic_dat[3:4], is.na))
descrip = matrix(c(means, sds, nas), nrow = 2,
ncol = 3)
colnames(descrip) = c('Mean', 'SD', 'Missing')
rownames(descrip) = c('Age', 'Fare')
# Get desriptives for categorical variables:
survivalcount = table(titanic_dat$Survived)
survivalnas = sum(sapply(titanic_dat[1], is.na))
survival = matrix(c(survivalcount, survivalnas),
nrow = 3, ncol = 1)
rownames(survival) = c('Perished', 'Survived',
'Missing')
colnames(survival) = 'Count'
gendercount = table(titanic_dat$Sex)
gendernas = sum(sapply(titanic_dat[2], is.na))
gender = matrix(c(gendercount, gendernas),
nrow = 3, ncol = 1)
rownames(gender) = c('Female', 'Male', 'Missing')
colnames(gender) = 'Count'
# Print all descriptives in a list
list(descrip, survival, gender)
# Question 2
############
# Run a logistic regression to predict whether
# passengers' survival (Survived [levels: 1 =
# survived; 0 = perished]) is predicted by cost of
# fare paid by passengers (Fare). Write out the
# logistic regression model that you have just fit.
fit0 = glm(Survived ~ Fare, family = binomial,
data = titanic_dat)
summary(fit0)
# Regression Model:
# log(p_hat / (1 - p_hat)) = b_0 + b_1 * Fare
# Question 3
############
# Interpret each coefficient from the analysis
# above in terms of both odds-ratios and
# probabilities.
# Setting plotting parameters
par(mfrow = c(1, 3))
# Creating a range of x values
x_range = seq(min(titanic_dat$Fare, na.rm = T),
max(titanic_dat$Fare, na.rm = T),
.01)
logit = coef(fit0)[1] + coef(fit0)[2] * x_range # Logits prediction
odds = exp(logit) # Odds prediction
prob = odds / (1 + odds) # Predicted Probability
# Plot predicted probabilities (sigmoidal)
plot(titanic_dat$Fare, titanic_dat$Survived,
ylab = "P(Survival)", xlab = 'Fare',
main = "Predicted Probability")
points(x_range, prob, pch = 16, col = "green")
# Plot odds (multiplicative)
plot(x_range, odds, pch = 16,col = "red",
xlab = 'Fare', ylab = "Odds of Survival",
main = 'Odds')
# Plot log odds (linear)
plot(x_range, logit, col = "black",
xlab = 'Fare', ylab = 'Log Odds of Survival',
main = 'Log Odds')
# To get to odds-ratios, we can take the exponential
# of our coefficients:
exp(coef(fit0))
# B0: When cost of fare is 0 (probably not meaningful),
# the odds of survival are 0.39.
# B1: For every one-unit increase in the price of fare
# to be a passenger on the Titanic, the odds of
# survival increase by a factor of 1.02. So for
# each unit increase in fare, the odds
# of survival are multiplied by 1.02.
# To interpret using probabilities:
exp(coef(fit0)) / (1 + exp(coef(fit0)))
# B0: Predicted probability of survival is .28 if fare
# is equal to 0 (again, probably not meaningful)
# B1: We could intepret this different ways, for example
# we might interpret the predicted probability for the
# mean value of fare:
mean(titanic_dat$Fare, na.rm = TRUE)
predict(fit0, list(Fare = mean(titanic_dat$Fare, na.rm = TRUE)),
type = 'resp')
# type = 'resp' will return predicted probabilities
# And we see that the predicted probability for
# a passenger who paid the average price of fare
# had a .39 probability of survival
# Question 4
############
# Plot the observed data with the logistic curve
# overlaid on the scatterplot.
# Setting plotting parameters
par(mfrow = c(1, 1))
with(titanic_dat, plot(Fare, Survived))
curve(predict(fit0, data.frame(Fare = x), type =
'resp'), add = T, col = 'red', lwd = 2)
# Question 5
############
# Evaluate the fit of your logistic regression
# model (hint: see slides 50-54).
# Error rate is defined as the proportion of cases
# for which the deterministic prediction is wrong
# So, if the model predicts a probability greater
# than .5 (greater than chance survival rate), but
# the passenger actually perished, then this is
# counted as error.
# Likewise, if the model predicts a probability
# less than .5 (less than chance survival rate),
# but the passenger actually survived, then this
# is also counted as error.
# Here, we run into an issue due to missing values.
# If you try to run this using the analyses conducted
# above:
mean((fitted(fit0) > 0.5 & titanic_dat$Survived == 0) |
(fitted(fit0) < 0.5 & titanic_dat$Survived == 1))
# ...you'll get NA because missing values were dropped
# when you fit the model, but they are still included
# in the original data set
# To get around this, we are going to create a complete
# dataset with no missing values and then re-run our
# analysis:
titanic_sub = titanic_dat[complete.cases(titanic_dat), ]
# Now rerun the model with complete data:
mod0 = glm(Survived ~ Fare, family = binomial,
data = titanic_sub)
# Now, we can evaluate the fit of our model:
mean((fitted(mod0) > 0.5 & titanic_sub$Survived == 0) |
(fitted(mod0) < 0.5 & titanic_sub$Survived == 1))
# So our error rate is 0.33. This is not great,
# but is better than chance (0.5)
# We can also compare this to a null model:
modnull = glm(Survived ~ 1, family = binomial, data = titanic_sub)
mean((fitted(modnull) > 0.5 & titanic_sub$Survived == 0) |
(fitted(modnull) < 0.5 & titanic_sub$Survived == 1))
# ...and we're doing slightly better than the
# null model which has an error rate of .41.
# We can also evaluate model deviance--a statistical
# summary of model fit, defined for generalized
# linear models to be an analogy to residual
# standard deviation
# Deviance is a measure of error: lower deviance
# means better fit to the data
# If a predictor that is simply random noise
# is added to a model, we expected deviance to
# decrease by 1, on average
# When informative predictor is added to a model,
# we expect deviance to decrease by more than 1.
# When k predictors are added to a model, we
# expect deviance to decrease by more than k.
# For models that are not multilevel, deviance
# is equal to -2 * log(likelihood function)
summary(mod0)
# Null deviance (model including intercept only)
# is 964.52 on 713 degrees of freedom
# Residual deviance (model adding in fare) is
# 901.25 on 712 degrees of freedom
# Deviance has decreased by
964.52 - 901.25 # which is more than we would
# expect if the predictor
# had been just noise.
# We can also evaluate AIC, which allows us to
# compare models while applying a penalty
# for model complexity (penalty is equal to
# +2p with p = k + 1)
summary(mod0)$deviance + 2 * (1 + 1)
# And this matches what we saw in our output:
summary(mod0)
summary(modnull)
# When evaluating AIC, lower is better and we
# see that the model including the predictor
# has a lower AIC value than does the null model
# Question 6
############
# Imagine that you have just remembered that women
# and children were given priority in boarding
# lifeboats. Repeat your logistic regression, this
# time adding in additional predictors of Age and
# Sex [levels: male; female]. Include an appropriate
# interaction term.
fit1 = glm(Survived ~ Fare + Age + Sex + Age*Sex,
family = binomial, data = titanic_dat)
summary(fit1)
# Question 7
############
# Interpret all coefficients from the analysis above
# in terms of both odds-ratios and probabilities.
# Odds Ratios:
exp(coef(fit1))
# B0: The predicted odds of survival for a female at age 0
# and Fare equal to 0 is 1.37 (probably not meaningful)
# B1: For every one-unit increase in the price of fare,
# the odds of survival increase by a factor of 1.01,
# when holding the other predictors constant.
# B2: For every one-year increase in age, the odds of of
# survival increase by a factor of 1.01, when holding
# the other predictors constant
# B3: This is the ratio of survival between males and
# females. Odds of Male Survival/Odds of Female Survival
# is about .27.
# So, males are about 1/4 as likely to survival as
# females, when holding the other predictors constant.
# We can reverse this ratio to make it more interpretable
# (i.e., compute the odds of female survival/odds of male
# survival) by dividing 1 by the coefficient like so:
1 / exp(coef(fit1)[4]) # And we see that females were about
# 3.69 times more likely to survive
# than males
# B4: Interaction between Sex and Age. This coefficient indicates
# that for each additional year of Age, the odds of survival
# change by a factor of 0.96 for males. Stated differently,
# for every additional year of age, males odds of survival
# decrease by...
1 / exp(coef(fit1)[5]) # ...a factor of 1.04, when holding the
# other predictors constant
# Probabilities:
exp(coef(fit1)) / (1 + exp(coef(fit1)))
# B0: The predicted probability of survival for a female at age 0
# and Fare equal to 0 is .58 (probably not meaningful)
# B1 & B2: We could just get the predicted probabilities when
# Fare and Age are set to their mean:
predict(fit1, list(Fare = mean(titanic_dat$Fare, na.rm = TRUE),
Age = mean(titanic_dat$Age, na.rm = TRUE),
Sex = 'male'), type = 'resp')
# ...and we get a predicted probability of Survival
# of .21 for a male of average age and cost of fare
# We could also look at differences in predicted
# probabilities between different points along the
# range of predictors
# Let's compute the predictive difference with respect
# to Age at the average value of Age and 0 Fare for males
age.mean.pred = predict(fit1, list(Fare = 0,
Age = mean(titanic_dat$Age, na.rm = TRUE),
Sex = 'male'), type = 'resp')
age.plus1.pred = predict(fit1, list(Fare = 0,
Age = (mean(titanic_dat$Age, na.rm = TRUE) + 1),
Sex = 'male'), type = 'resp')
age.mean.pred - age.plus1.pred
# In the Gelman and Hill text, they point out that the divide by 4
# rule gives you a good estimate of the predicted difference for a
# one unit change in the predictor
coef(fit1)[3] / 4 # and we see that the estimate is fairly similar to
# what we got above
# So, for males who paid a fare of 0, each
# additional year of age made a 0.3% positive difference
# in probability of survival.
# B3: We can get the difference in predicted probabilites betwee
# males and females who paid 0 fare, at age 0:
Male = predict(fit1, list(Fare = 0,
Age = 0,
Sex = 'male'), type = 'resp')
Female = predict(fit1, list(Fare = 0,
Age = 0,
Sex = 'female'), type = 'resp')
Male - Female
# And we see that males who paid 0 fare, and are age 0,
# decreased by about 30.7% in probability of survival
# relative to females
# B4: Interaction between Sex and Age. This makes my head
# explode, so we're probably best sticking to explaining
# this relationship in terms of odds here :)
# (Also, see plot at end)
# Question 8
############
# Use the predict() function to determine the odds
# of your own survival had you been a passenger on
# the Titanic (you can use the mean value for Fare).
# Note that predict() returns the predicted values
# at the linear scale, i.e., on the logit scale.
predict(fit1, list(Fare = mean(titanic_dat$Fare, na.rm = T),
Age = 31,
Sex = 'female'), type = 'resp')
# Your TA's predicted probability of survival on
# the Titanic would be .76. (Woo?)
# Question 9
############
# We mentioned that missing values may be important
# to note early in this assignment...this is because
# the likelihood ratio test will not work if your
# models
# Let's look at how missing data is handled in lm()
?glm # So glm() by default will use na.omit
# So different numbers are being excluded from each of
# our two models
# Now rerun the models with complete data:
mod0 = glm(Survived ~ Fare, family = binomial,
data = titanic_sub)
mod1 = glm(Survived ~ Fare + Age + Sex + Age*Sex,
family = binomial, data = titanic_sub)
# install.packages('lmtest')
library(lmtest)
lrtest(mod0, mod1)
anova(mod0, mod1, test = 'Chisq')
# BONUS
#######
# Create a plot to illustrate the interaction effect
# found in question 6.
summary(fit1)
# Save coefficients
b0 = fit1$coef[1] # intercept
Fare = fit1$coef[2] # slope for Fare
Age = fit1$coef[3] # slope for Age
male = fit1$coef[4] # slope for Sex
AgeSex = fit1$coef[5]# slope for interaction
# Create a range of ages for plotting
Age_range = seq(min(titanic_dat$Age, na.rm = T),
max(titanic_dat$Age, na.rm = T), .01)
# Set Fare to be equal to the mean
Fare_val = mean(titanic_dat$Fare, na.rm = T)
# Equation for females in logits:
female_logits = b0 + Age * Age_range + Fare * Fare_val +
male * 0 + AgeSex * Age_range * 0
# Note that these last two terms drop
# out of the equation
# Equation for males in logits:
male_logits = b0 + Age * Age_range + Fare * Fare_val +
male * 1 + AgeSex * Age_range * 1
# Now the last two terms are included
# Convert from logits to probability scale
Female = exp(female_logits) / (1 + exp(female_logits))
Male = exp(male_logits) / (1 + exp(male_logits))
# Plot the age range against predicted probabilities
# for females
plot(Age_range, Female, ylim = c(0, 1), type = 'l',
lwd = 3, lty = 2, col = 'orangered', xlab = 'Age',
ylab = 'P(Survival)',
main = 'Probability of Survival on the Titanic')
# Add in the predicted probabilities for males
lines(Age_range, Male, type = 'l',
lwd = 3, lty = 3, col = 'turquoise2')
# Add a legend to the plot
legend('topleft', c('Female', 'Male'), col =
c('orangered', 'turquoise2'), lwd = 3,
lty = c(2, 3), horiz = T)