USairpollution Data Analysis
이진규
2018-10-01
1 Load libraries
library(tidyverse)
library(gridExtra)
library(corrplot)
library(MVA)2 Read Data
HSAUR2에 내장되어 있는 USairpollution 데이터로 간단한 시각화와 분석을 해보자.
data <- as.data.frame(HSAUR2::USairpollution)
# View first 6 rows
head(data)## SO2 temp manu popul wind precip predays
## Albany 46 47.6 44 116 8.8 33.36 135
## Albuquerque 11 56.8 46 244 8.9 7.77 58
## Atlanta 24 61.5 368 497 9.1 48.34 115
## Baltimore 47 55.0 625 905 9.6 41.31 111
## Buffalo 11 47.1 391 463 12.4 36.11 166
## Charleston 31 55.2 35 71 6.5 40.75 148# Structure of data set
str(data)## 'data.frame': 41 obs. of 7 variables:
## $ SO2 : int 46 11 24 47 11 31 110 23 65 26 ...
## $ temp : num 47.6 56.8 61.5 55 47.1 55.2 50.6 54 49.7 51.5 ...
## $ manu : int 44 46 368 625 391 35 3344 462 1007 266 ...
## $ popul : int 116 244 497 905 463 71 3369 453 751 540 ...
## $ wind : num 8.8 8.9 9.1 9.6 12.4 6.5 10.4 7.1 10.9 8.6 ...
## $ precip : num 33.36 7.77 48.34 41.31 36.11 ...
## $ predays: int 135 58 115 111 166 148 122 132 155 134 ...# Data Variables
colnames(data)## [1] "SO2" "temp" "manu" "popul" "wind" "precip" "predays"# Missing value
table(is.na(data))##
## FALSE
## 287# Basic statistics
summary(data)## SO2 temp manu popul
## Min. : 8.00 Min. :43.50 Min. : 35.0 Min. : 71.0
## 1st Qu.: 13.00 1st Qu.:50.60 1st Qu.: 181.0 1st Qu.: 299.0
## Median : 26.00 Median :54.60 Median : 347.0 Median : 515.0
## Mean : 30.05 Mean :55.76 Mean : 463.1 Mean : 608.6
## 3rd Qu.: 35.00 3rd Qu.:59.30 3rd Qu.: 462.0 3rd Qu.: 717.0
## Max. :110.00 Max. :75.50 Max. :3344.0 Max. :3369.0
## wind precip predays
## Min. : 6.000 Min. : 7.05 Min. : 36.0
## 1st Qu.: 8.700 1st Qu.:30.96 1st Qu.:103.0
## Median : 9.300 Median :38.74 Median :115.0
## Mean : 9.444 Mean :36.77 Mean :113.9
## 3rd Qu.:10.600 3rd Qu.:43.11 3rd Qu.:128.0
## Max. :12.700 Max. :59.80 Max. :166.0USairpollution data는 총 41개의 관측치와 7개의 변수들로 구성되었고 변수들은 모두 숫자형(정수 및 실수)의 타입을 가졌다.
또한 결측값은 존재하지 않았고 데이터 변수명과 기초 통계량을 확인하였다.
3 Data Visualization
3.1 Histogram
총 7개의 연속형 변수에 대해서 Histogram을 확인해보자.
p1 <- ggplot(data, aes(x=SO2)) +
geom_histogram(binwidth=4, fill="palegreen", colour="black") +
labs(x="SO2")
p2 <- ggplot(data, aes(x=temp)) +
geom_histogram(binwidth=1, fill="palegreen", colour="black") +
labs(x="temp")
p3 <- ggplot(data, aes(x=manu)) +
geom_histogram(binwidth=100, fill="palegreen", colour="black") +
labs(x="manu")
p4 <- ggplot(data, aes(x=popul)) +
geom_histogram(binwidth=100, fill="palegreen", colour="black") +
labs(x="popul")
p5 <- ggplot(data, aes(x=wind)) +
geom_histogram(binwidth=0.5, fill="palegreen", colour="black") +
labs(x="wind")
p6 <- ggplot(data, aes(x=precip)) +
geom_histogram(binwidth=4, fill="palegreen", colour="black") +
labs(x="precip")
p7 <- ggplot(data, aes(x=predays)) +
geom_histogram(binwidth=4, fill="palegreen", colour="black") +
labs(x="predays")
grid.arrange(p1, p2, p3, p4, p5, p6, p7,
layout_matrix = cbind(c(1,4,7), c(2,5,7), c(3,6,7)))
3.2 Kernel Density Curve
Histogram은 빈도를 가지고 그리지만 Kernel Density Curve는 확률을 가지고 그린다.
p8 <- ggplot(data, aes(x=SO2)) +
geom_density(fill="darkslategray", alpha=0.5) +
labs(x="SO2")
p9 <- ggplot(data, aes(x=temp)) +
geom_density(fill="darkslategray", alpha=0.5) +
labs(x="temp")
p10 <- ggplot(data, aes(x=manu)) +
geom_density(fill="darkslategray", alpha=0.5) +
labs(x="manu")
p11 <- ggplot(data, aes(x=popul)) +
geom_density(fill="darkslategray", alpha=0.5) +
labs(x="popul")
p12 <- ggplot(data, aes(x=wind)) +
geom_density(fill="darkslategray", alpha=0.5) +
labs(x="wind")
p13 <- ggplot(data, aes(x=precip)) +
geom_density(fill="darkslategray", alpha=0.5) +
labs(x="precip")
p14 <- ggplot(data, aes(x=predays)) +
geom_density(fill="darkslategray", alpha=0.5) +
labs(x="predays")
grid.arrange(p8, p9, p10, p11, p12, p13, p14,
layout_matrix=cbind(c(1,4,7), c(2,5,7), c(3,6,7)))
3.3 Boxplot
연속형 변수에 대해서 다양한 정보를 알 수 있는 Boxplot
p15 <- ggplot(data, aes(x=1, y=SO2)) + # x=1(임의의 값) 지정해줘야함.
geom_boxplot(fill="paleturquoise", width=0.3,
outlier.color = "red", outlier.shape=16) +
coord_flip() +
scale_x_continuous(breaks = NULL) + # x축 이름 생략
theme(axis.title.y= element_blank()) + # y축 구분자 생략(coord_flip())
labs(y="SO2")
p16 <- ggplot(data, aes(x=1, temp)) +
geom_boxplot(fill="paleturquoise", width=0.3,
outlier.color = "red", outlier.shape=16) +
coord_flip() +
scale_x_continuous(breaks = NULL) +
theme(axis.title.y= element_blank()) +
labs(y="temp")
p17 <- ggplot(data, aes(x=1, y=manu)) +
geom_boxplot(fill="paleturquoise", width=0.3,
outlier.color = "red", outlier.shape=16) +
coord_flip() +
scale_x_continuous(breaks = NULL) +
theme(axis.title.y= element_blank()) +
labs(y="manu")
p18 <- ggplot(data, aes(x=1, y=popul)) +
geom_boxplot(fill="paleturquoise", width=0.3,
outlier.color = "red", outlier.shape=16) +
coord_flip() +
scale_x_continuous(breaks = NULL) +
theme(axis.title.y= element_blank()) +
labs(y="popul")
p19 <- ggplot(data, aes(x=1, y=wind)) +
geom_boxplot(fill="paleturquoise", width=0.3,
outlier.color = "red", outlier.shape=16) +
coord_flip() +
scale_x_continuous(breaks = NULL) +
theme(axis.title.y= element_blank()) +
labs(y="wind")
p20 <- ggplot(data, aes(x=1, y=precip)) +
geom_boxplot(fill="paleturquoise", width=0.3,
outlier.color = "red", outlier.shape=16) +
coord_flip() +
scale_x_continuous(breaks = NULL) +
theme(axis.title.y= element_blank()) +
labs(y="precip")
p21 <- ggplot(data, aes(x=1, y=predays)) +
geom_boxplot(fill="paleturquoise", width=0.3,
outlier.color = "red", outlier.shape=16) +
coord_flip() +
scale_x_continuous(breaks = NULL) +
theme(axis.title.y= element_blank()) +
labs(y="predays")
grid.arrange(p15, p16, p17, p18, p19, p20, p21,
layout_matrix=cbind(c(1,4,7), c(2,5,7), c(3,6,7)))
4 SO2 vs manu
USairpollution의 변수들의 상관관계를 알아보자.
data_cor <- cor(data)
data_cor## SO2 temp manu popul wind
## SO2 1.00000000 -0.43360020 0.64476873 0.49377958 0.09469045
## temp -0.43360020 1.00000000 -0.19004216 -0.06267813 -0.34973963
## manu 0.64476873 -0.19004216 1.00000000 0.95526935 0.23794683
## popul 0.49377958 -0.06267813 0.95526935 1.00000000 0.21264375
## wind 0.09469045 -0.34973963 0.23794683 0.21264375 1.00000000
## precip 0.05429434 0.38625342 -0.03241688 -0.02611873 -0.01299438
## predays 0.36956363 -0.43024212 0.13182930 0.04208319 0.16410559
## precip predays
## SO2 0.05429434 0.36956363
## temp 0.38625342 -0.43024212
## manu -0.03241688 0.13182930
## popul -0.02611873 0.04208319
## wind -0.01299438 0.16410559
## precip 1.00000000 0.49609671
## predays 0.49609671 1.000000004.1 Correlation matrix
모든 변수들의 상관계수를 한 눈에 알아볼 수 있는 상관계수 행렬.
corrplot(data_cor)
다음과 같이 더 예쁘게 그릴 수 있다.
corrplot(data_cor,
method = "color", # 색깔로 표현
type = "lower", # 왼쪽 아래 행렬만 표시
order = "hclust", # 유사한 상관계수끼리 군집화
addCoef.col = "black", # 상관계수 색깔
tl.col = "black", # 변수명 색깔
tl.srt = 45, # 변수명 45도 기울임
diag = F) # 대각 행렬
manu와 popul이 0.96으로 상관계수가 가장 높다.
위의 관계는 수업자료에 다룬 내용이므로 상관계수가 그 다음으로 높은 manu와 SO2를 다뤄보자.
4.2 Scatter plot
ggplot(data, aes(x=SO2, y=manu)) +
geom_point()
popul 변수를 원의 크기로 하면 Bubble chart도 그릴 수 있다.
ggplot(data, aes(x=SO2, y=manu)) +
geom_point(aes(size=popul),shape=21, colour="black", fill="yellow", alpha=0.5) +
scale_size_area()
산점도에 label로 rowname을 추가
ggplot(data, aes(x=SO2, y=manu)) +
geom_point(shape=20, size=3) +
geom_text(aes(label=row.names(data), size=1, vjust=-1, hjust=0)) +
theme(legend.position = "none")
Chicago와 Philadelphia 같은 도시들이 이상치로 의심된다.
4.3 Bivariate boxplot
data1 <- data %>%
select(SO2, manu)
bvbox(data1)
text(data1, labels=row.names(data1), cex=0.7)
이상치를 제거해보자.
outcity <-
match(c("Chicago", "Detroit","Cleveland", "Philadelphia", "Providence"), rownames(USairpollution)) # 7, 14, 9, 30, 33
# 이상치를 제거한 후 상관계수
round(with(data, cor(manu[-outcity], popul[-outcity])), 2)## [1] 0.810.64 -> 0.81로 많이 커졌음을 알 수 있다.
5 linear regression
manu를 독립변수, SO2를 종속변수로 간단하게 회귀분석을 해보자.
data_lm <- lm(manu ~ SO2, data=data)
summary(data_lm)##
## Call:
## lm(formula = manu ~ SO2, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1109.95 -157.12 -12.86 184.22 1643.40
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.006 111.517 -0.018 0.986
## SO2 15.478 2.938 5.268 5.36e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 436.2 on 39 degrees of freedom
## Multiple R-squared: 0.4157, Adjusted R-squared: 0.4007
## F-statistic: 27.75 on 1 and 39 DF, p-value: 5.363e-06manu = -2.006 + 15.478 * SO2
추정값은 15.478, 표준편차는 2.938, t값은 5.268, p값은 5.36e-06이므로 0으로 봐도 무방하다.
따라서 manu에 대한 SO2의 선형 효과는 통계적으로 유의하다.
5.1 Scatter plot + lm
산점도에 회귀식을 표시하면 다음과 같다.
# Scatter plot + lm
ggplot(data, aes(x=SO2, y=manu)) +
geom_point() +
geom_smooth(method="lm", se=TRUE)
6 Summary
간단하게 데이터 분석과 시각화를 해봤는데 데이터의 수가 많지 않았고 데이터 변수들의 타입도 단순 연속형이었기 때문에 다양한 분석과 시각화를 할 수 없었다.
따라서 복잡한 분석을 통해 유의한 결과를 얻기 보다는 단순한 분석과 시각화를 통해 접근하는 것이 좋아보였다.