-
-
Notifications
You must be signed in to change notification settings - Fork 2
/
Module-6-Example-3.R
110 lines (67 loc) · 2.75 KB
/
Module-6-Example-3.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
# setwd("SET THE Working Director to THE PATH TO THIS DIRECTORY")
# install.packages("aod")
# install.packages("pROC")
library(aod)
library(stats)
library(pROC)
data<-read.csv("Datasets/cevent.csv")
attach(data)
# print a small part of the data
head(data)
m <- glm(data$event ~ data$chol, family = "binomial")
summary(m)
# ORs per 1 unit increase
# same as calculation by hand (OR): exp(0.02119)
exp(m$coefficients[2])
exp(cbind(OR = coef(m), confint.default(m)))
# OR per 10 unit increase
exp(m$coefficients[2]*10)
exp((m$coefficients[2]-qnorm(0.975)*summary(m)$coefficients[2,2])*10)
exp((m$coefficients[2]+qnorm(0.975)*summary(m)$coefficients[2,2])*10)
# The same cab be produced by using the following
exp(10 * cbind(OR = coef(m), confint.default(m)))
# predict risk for each patient
risk <-predict(m, type=c("response"))
risk
# predict risk for patient with chol of 190: exp(-3.12716+0.02119*190)/(1+exp(-3.12716+0.02119*190))
risk[41]
exp(m$coefficients[1] + m$coefficients[2]*190) / (1+exp(m$coefficients[1] + m$coefficients[2]*190))
# multiple logistic regression
data$male <- ifelse(data$sex =="M", 1, 0)
m2<-glm(data$event ~ data$chol + data$male + data$age, family="binomial")
summary(m2)
# ORs per 1 unit increase
exp(cbind(OR = coef(m2), confint.default(m2)))
# ROC curve
# install.package("pROC")
library(pROC)
# using model with chol and sex and age
data$prob <-predict(m2, type=c("response"))
# ROC Curve
g <- roc(data$event ~ data$prob)
# Get the Area under the curve
# c-statistics
g
print(g)
# Plot the ROC Curve.
plot(g)
# Just print the results
roc(data$event ~ data$prob)
# or plot the graphs as well
roc(data$event ~ data$prob, plot=TRUE)
# To get ride of the padding of graph, you can use the par function to set some enviromental variables.
par(pty="s")
roc(data$event ~ data$prob, plot=TRUE)
# Another way would be to see the x axis as 1-Specificity
plot(1- g$specificities, g$sensitivities, type="l", xlab="1-Specifity", ylab="Sensivity", main="ROC Curve")
abline(a=0, b=1)
grid()
# If you want to have a normal x-axis from zero to one
roc(data$event ~ data$prob, plot=TRUE, legacy.axes=TRUE)
# if you want to get values in precentages
roc(data$event ~ data$prob, plot=TRUE, legacy.axes=TRUE, percent=TRUE)
# If you want to understand better specificities and sensitivities.
# These are just false positive and true negatives.
roc(data$event ~ data$prob, plot=TRUE, legacy.axes=TRUE, percent=TRUE, xlab="False Positive (%)", ylab="True Positive (%)")
# chaning the color
roc(data$event ~ data$prob, plot=TRUE, legacy.axes=TRUE, percent=TRUE, xlab="False Positive (%)", ylab="True Positive (%)", col="blue", lwd=4)