forked from KimmoVehkalahti/IODS-project
-
Notifications
You must be signed in to change notification settings - Fork 0
/
chapter4.Rmd
235 lines (167 loc) · 6.93 KB
/
chapter4.Rmd
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
# Clustering and Classification
This week we are going to see the clustering and classification of Boston data from the MASS package.
## Source
Harrison, D. and Rubinfeld, D.L. (1978) Hedonic prices and the demand for clean air. J. Environ. Economics and Management 5, 81–102.
Belsley D.A., Kuh, E. and Welsch, R.E. (1980) Regression Diagnostics. Identifying Influential Data and Sources of Collinearity. New York: Wiley.
## load the Boston data and explore
```{r}
library(MASS)
library(tidyverse)
library(corrplot)
library(ggplot2)
library(GGally)
# access the MASS package
library(MASS)
# load the data
data("Boston")
# explore the dataset
str(Boston)
summary(Boston)
# plot matrix of the variables
pairs(Boston)
```
The Boston housing dataset consists of 506 rows with 14 columns from suburbs of Boston.These compare the age of the buildings, prices of the buildings, area nitrogen oxide concentrations, average number of rooms for dwelling, proportion of non-retail business acres per town etc.
##Graphical Overview
```{r}
library(tidyverse)
library(corrplot)
# calculate the correlation matrix and round it
cor_matrix<-cor(Boston) %>% round(digits = 2)
# print the correlation matrix
cor_matrix
# visualize the correlation matrix
corrplot(cor_matrix, method="circle", type="upper", cl.pos="b", tl.pos="d", tl.cex = 0.6)
```
Above is the correlation chart of the values. Size of the circle varies according to correlation coefficents. The color of the circle indicates whether it is negatively or positively correlating.
In here it's visible that rad (index of accessibility to radial highways) correlates positively to tax (full-value property-tax rate per \$10,000.) and lstat(lower status of the population (percent)) correlates negatively with medv (median value of owner-occupied homes in \$1000s). High positive correlation can be seen between nox and indus.
##Standardizing
```{r}
# center and standardize variables
boston_scaled <- scale(Boston)
# summaries of the scaled variables
summary(boston_scaled)
# class of the boston_scaled object
class(boston_scaled)
# change the object to data frame
boston_scaled <- as.data.frame(boston_scaled)
```
All the data is standardized now using scale function. Compared to earlier summary, the values have changed from positive to negative and minimized in size.
```{r}
# summary of the scaled crime rate
summary(boston_scaled$crim)
# create a quantile vector of crim and print it
bins <- quantile(boston_scaled$crim)
bins
# create a categorical variable 'crime'
crime <- cut(boston_scaled$crim, breaks = bins, include.lowest = TRUE, labels = c("low", "med_low", "med_high", "high"))
# look at the table of the new factor crime
table(crime)
# remove original crim from the dataset
boston_scaled <- dplyr::select(boston_scaled, -crim)
# add the new categorical value to scaled data
boston_scaled <- data.frame(boston_scaled, crime)
```
##Creating training and test datasets
```{r}
# number of rows in the Boston dataset
n <- nrow(boston_scaled)
# choose randomly 80% of the rows
ind <- sample(n, size = n * 0.8)
# create train set
train <- boston_scaled[ind,]
# create test set
test <- boston_scaled[-ind,]
# save the correct classes from test data
correct_classes <- test$crime
# remove the crime variable from test data
test <- dplyr::select(test, -crime)
```
##Fitting linear discriminat analysis on training set
```{r}
# linear discriminant analysis
lda.fit <- lda(crime ~ ., data = train)
# print the lda.fit object
lda.fit
# the function for lda biplot arrows
lda.arrows <- function(x, myscale = 1, arrow_heads = 0.1, color = "orange", tex = 0.75, choices = c(1,2)){
heads <- coef(x)
arrows(x0 = 0, y0 = 0,
x1 = myscale * heads[,choices[1]],
y1 = myscale * heads[,choices[2]], col=color, length = arrow_heads)
text(myscale * heads[,choices], labels = row.names(heads),
cex = tex, col=color, pos=3)
}
# target classes as numeric
classes <- as.numeric(train$crime)
# plot the lda results
plot(lda.fit, dimen = 2, col = classes, pch = classes)
lda.arrows(lda.fit, myscale = 2)
```
The plot shows that crime rate stands out from the rest of the data.
```{r}
# predict classes with test data
lda.pred <- predict(lda.fit, newdata = test)
# cross tabulate the results
table(correct = correct_classes, predicted = lda.pred$class)
```
##K-Means clustering
```{r reload, echo=FALSE}
data("Boston")
boston_scaled <- scale(Boston)
dist_eu <- dist(Boston)
km <-kmeans(dist_eu, centers = 15)
pairs(Boston, col = km$cluster)
```
The plot is very difficult to understand sos let us interpret it using within cluster sum of squared (wcss).
```{r wcss, echo=FALSE}
set.seed(123)
dist_eu <- dist(Boston)
k_max <- 15
twcss <- sapply(1:k_max, function(k){kmeans(dist_eu, k)$tot.withinss})
plot(1:k_max, twcss, type='b')
```
Form the above graph we can see that the highest drop of total wcss is at cluster number 2. So the best cluster number is 2.
```{r kmeans, echo=FALSE}
km <-kmeans(dist_eu, centers = 2)
pairs(Boston, col = km$cluster)
```
The new graph shows only 2 clusters.
##Bonus
The dataset is reloaded and standardized. The graph below has 4 cluster centers.
```{r bonus, echo=FALSE}
#install.packages("plotly")
library(plotly)
library(MASS)
data("Boston")
boston_scaled <- scale(Boston)
dist_eu <- dist(Boston)
km <-kmeans(dist_eu, centers = 4)
pairs(Boston, col = km$cluster)
lda.fit2 <- lda(km$cluster ~., data = Boston)
lda.arrows <- function(x, myscale = 1, arrow_heads = 0.1, color = "red", tex = 0.75, choices = c(1,2)){
heads <- coef(x)
arrows(x0 = 0, y0 = 0,
x1 = myscale * heads[,choices[1]],
y1 = myscale * heads[,choices[2]], col=color, length = arrow_heads)
text(myscale * heads[,choices], labels = row.names(heads),
cex = tex, col=color, pos=3)
}
plot(lda.fit2, dimen = 2, col= classes, pch=classes)
lda.arrows(lda.fit2, myscale = 2)
```
The graph above is LDA model with the K-means clusters as target variable and Boston data frame as data. The variable with the longest arrow is nitrogen oxides concentration (nox) and it divides the data frames observations into two separate areas. We can see more closely below.
```{r bonusclose, echo=FALSE}
plot(lda.fit2, dimen = 2, col= classes, pch=classes)
lda.arrows(lda.fit2, myscale = 10)
```
The two variables showing to be linear separators are variables index of accessibility to radial highways (rad) and Charles River dummy variable (chas) that separate the observations from cluster one from the rest of the group.
##Super bonus
```{r superbonus, echo=FALSE}
model_predictors <- dplyr::select(train, -crime)
dim(model_predictors)
dim(lda.fit$scaling)
matrix_product <- as.matrix(model_predictors) %*% lda.fit$scaling
matrix_product <- as.data.frame(matrix_product)
plot_ly(x = matrix_product$LD1, y = matrix_product$LD2, z = matrix_product$LD3, type= 'scatter3d', mode='markers', color = train$crime)
```
The three matrix produts from previous LDA is inserted into plot_ly function. We can see the above 3D graph which can be moved.