정보
-
업무명 : 서울대 통계 연구소 R을 이용한 빅데이터 분석 교육 연수
-
작성자 : 이상호
-
작성일 : 2020-05-12
-
설 명 :
-
수정이력 :
내용
[개요]
-
안녕하세요? 기상 연구 및 웹 개발을 담당하고 있는 해솔입니다.
-
오늘은 통계 연구소 R을 이용한 빅데이터 분석 교육 연수 내용을 다루고자 합니다.
-
이 교육은 5일에 걸쳐 진행하였고 각 시간마다 이론 및 실습을 병행하여 진행했습니다.
-
특히 배운 지식을 활용하여 학술 논문 또는 보고서 작성 시 많은 도움이 되었습니다.
-
그에 따라 "교육 일정, 실습 자료, 관련 자료"순으로 소개해 드리겠습니다.
[특징]
-
R을 이용한 빅데이터 교육에 대한 이해를 돕기위해 작성
[기능]
-
교육 일정
-
실습 자료
-
관련 자료
[사용 OS]
-
Windows 10
[사용 언어]
-
R v4.0.0
세부 내용
[교육 일정]
-
강좌 구성
-
1일 : R tidyverse 패키지를 사용한 데이터 가공
-
2일 : 데이터 시각화 (Data Visualization)
-
3일 : 회귀 분석 (Regression Analysis)
-
4일 : 분류, 군집 (Classification, Clustering)
-
5일 : 추천 알고리즘 (Recommendation System)
-
-
강의 시간표
-
09:30-12:30 : 강의
-
13:30-15:30 : 실습
-
실습 자료
-
소스 코드에 대한 설명은 주석을 참조하시기 바랍니다.
1일 : R tidyverse 패키지를 사용한 데이터 가공
- 필수 패키지 설치 및 읽기
#################################################################################################
### ### ###
### Description : This is for the special course on "Big Data" of SNU-SRI. ###
### This concentrates on "Data Wrangling" with R package "dplyr" and "tidyr". ###
### ###
### Reference : Jaimie (Jaimyoung) Kwon's github (and lecture note) ###
### ###
#################################################################################################
# 필수 패키지 설치
install.packages("ggplot2")
install.packages("dplyr")
install.packages("tidyr")
install.packages("data.table")
# 학습을 위해 필요한 패키지 설치
install.packages("lubridate")
install.packages("Lahman")
# 데이터 설치
install.packages("nycflights13")
install.packages("babynames")
# 패키지 로딩
library(ggplot2)
library(dplyr)
library(tidyr)
library(data.table)
-
실습 자료
#################################################################################################
### ###
### Description : This is for the special course on "Big Data" of SNU-SRI. ###
### This concentrates on "Data Wrangling" with R package "dplyr" and "tidyr". ###
### ###
### Reference : Jaimie (Jaimyoung) Kwon's github (and lecture note) ###
### ###
#################################################################################################
#################################################################################################
########## dplyr ##########
#################################################################################################
library(nycflights13)
help(package="nycflights13")
# 데이터를 transpose 취해서 보여준다. 즉, 각 설명변수(열)들이 어떤 값들이 나왔는지 순서대로 보여주는 것이다.
# 이것을 이용하면 데이터프레임 안에 있는 모든 열들을 콘솔에서 편하게 볼 수 있다.
# str 과 다른 점은 콘솔 상에서 가능한 최대한의 값들을 보여준다.
glimpse(flights)
str(flights)
# 연습문제 12/43
glimpse(airlines)
glimpse(airports)
glimpse(planes)
glimpse(weather)
flights
print(flights)
class(flights)
# x <- data.frame(c(1,1))
# class(x)
# class(tbl_df(x))
################################# 고급 연습문제(13/43) #################################
head(as.data.frame(flights))
##################################### 1. filter() #####################################
# 설명 : 주어진 조건을 만족하는 행들을 추려낸다.
# 기초적인 사용법 : filter(주어진 자료(데이터프레임), 필터링하고자 하는 조건)
# 1월 1일에 출발한 항공편 자료를 찾는 명령
filter(flights, month == 1, day == 1)
filter(flights, month == 1 & day == 1)
# 1월과 2월에 출발한 항공편 자료를 찾는 명령
filter(flights, month == 1 | month == 2)
### slice() : 행을 행 번호로 추려낸다.
# flights[1:10,] ==
slice(flights, 1:10)
##################################### 2. arrange() #####################################
# 설명 : 행을 변수들의 오름차순으로 정렬한다. 즉 지정한 열을 기준으로 가장 작은 값에서 시작하여
# 가장 큰 값으로 끝나는 순서대로 데이터를 정렬한다.
# 기초적인 사용법 : arrange(주어진 자료(데이터프레임), 오름차순으로 정렬하고자 하는 열(기준))
# flights 자료를 year, month, day 순으로 정렬하는 명령
# flights[order(flights$year, flights$month, flights$day),] ==
arrange(flights, year, month, day)
### desc() : 내림차순으로 정렬하도록 바꾸는 함수이다.
# arr_delay가 가장 큰 것부터 정렬하는 명령
arrange(flights, desc(arr_delay))
##################################### 3. select() #####################################
# 설명 : 필요한 열을 선택한다. 열 이름을 써주는 연산이 가장 흔히 쓰인다.
# 기초적인 사용법 : select(주어진 자료(데이터프레임), 추출하고자 하는 열의 이름들)
# year, month, day라는 이름을 가진 열을 추출하는 명령
select(flights, year, month, day)
# 열의 범위를 지정하여 열을 추출하는 명령
select(flights, year:arr_time)
# year, month, day라는 이름을 가진 열을 제외하고 나머지를 추출하는 명령
select(flights, -year, -month, -day)
# 열의 범위를 지정하여 그 열들을 제외한 나머지 열들을 추출하는 명령
select(flights, -(year:arr_time))
### + starts_with() : 정해진 문자열로 시작하는 열들을 모두 지정할 때 사용
# arr로 시작하는 이름을 가진 열들을 모두 제외하는 명령
select(flights, -starts_with("arr"))
### + ends_with() : 정해진 문자열로 끝나는 열들을 모두 지정할 때 사용
# delay로 시작하는 이름을 가진 열들을 모두 추출하는 명령
select(flights, ends_with("delay"))
### + matches() : 정해진 형식으로 구성된 이름을 가진 열들을 모두 지정할 때 사용
# 무엇_무엇 꼴의 이름을 가진 열들을 모두 추출하는 명령
select(flights, matches("._."))
select(flights, matches(".r_t."))
### + contains() : 정해진 문자열이 포함된 이름을 가진 열들을 모두 지정할 때 사용
# _이 들어가는 이름을 가진 열들을 모두 추출하는 명령
select(flights, contains("_"))
select(flights, contains("_."))
## ?select 로 추가적으로 확인할 수 있다.
### rename() : 열의 이름을 바꿔준다.
rename(flights, tail_num = tailnum)
##################################### 4. mutate() #####################################
# 설명 : 기존의 열(들)의 값을 이용해 새로운 열을 추가한다.
# 함수에서 생성한 열을 바로 사용할 수 있다.
# 기초적인 사용법 : mutate(주어진 자료(데이터프레임), 생성하고자 하는 열 계산)
mutate(flights, gain = arr_delay - dep_delay, gain_per_hour = gain / (air_time / 60), speed = distance / air_time * 60)
# transform() 함수와 비교하면?
mutate(flights, gain = arr_delay - dep_delay, gain_per_hour = gain / (air_time / 60))
transform(flights, gain = arr_delay - dep_delay, gain_per_hour = gain / (air_time / 60)) #error
##################################### 5. summarize() #####################################
# 설명 : 주어진 자료(데이터프레임)을 한 줄(행)으로 요약한다. 즉 요약통계량을 계산한다.
# summarise() 함수도 같은 역할을 한다.
# 기초적인 사용법 : summarize(주어진 자료(데이터프레임), 원하는 요약통계량)
# dep_delay의 평균을 계산하는 명령
summarize(flights, delay = mean(dep_delay, na.rm = TRUE))
# 요약통계함수를 사용한다. 요약통계함수는 벡터값을 입력으로 받아 하나의 숫자를 출력한다.
# 기본 R 패키지의 요약통계함수는 min(), max(), mean(), sum(), sd(), median(), IQR() 등이 있다.
# dplyr 패키지에서 제공하는 요약통계함수는 n(), n_distinct(), first(), last(), nth() 이 있다.
# group_by() 함수와 함께 이용하면 매우 편리하다.
############################## 6. sample_n(), sample_frac() ##############################
# 설명 : sample_n() 함수는 정해진 숫자의 열을 랜덤샘플한다.
# sample_frac() 함수는 정해진 비율의 열을 랜덤샘플한다.
# 기초적인 사용법 : sample_n(주어진 자료(데이터프레임), 추출하고자 하는 열의 개수)
# sample_frac(주어진 자료(데이터프레임),
# 추출하고자 하는 열의 개수의 전체 열의 개수에 대한 비율)
# 10줄의 열을 랜덤샘플하는 명령
sample_n(flights, 10)
# 1%의 열을 랜덤샘플하는 명령
sample_frac(flights, 0.01)
# 비복원추출이 default이지만 replace 옵션으로 복원추출할 수도 있다.
# weight 옵션으로 가중치를 정할 수 있다.
# 재현가능한 연구를 위해서는 set.seed해줄 수 있다.
##################################### 7. distinct() #####################################
# 설명 : 주어진 자료에서 고유한 행을 찾아낸다.
# 기초적인 사용법 : distinct(주어진 자료(데이터프레임) 혹은 자료에서 선택한 열)
# 서로 다른 기체번호를 찾는 명령
distinct(select(flights, tailnum))
# 서로 다른 노선을 찾는 명령
distinct(select(flights, origin, dest))
############################# Language of data manipulation #############################
# (1) 행정렬 arrange()
# (2) 행선택 filter()
# (3) 열선택 select()
# (4) 새로운 변수 계산 mutate()
# (5) 요약통계량 계산 summarise()
########################### group_by() 함수를 이용한 그룹연산 ###########################
# 설명 : 데이터셋을 그룹으로 나눈다.
# 기초적인 사용법 : group_by(주어진 자료(데이터프레임), 그룹으로 나누고자 하는 기준이 되는 열(들))
flights
groups(flights)
(by_tailnum = group_by(flights, tailnum))
groups(by_tailnum)
### 연동사용법
## select() 는 그룹변수를 항상 포함한다.
select(by_tailnum, year, month, day)
## arrange() 는 그룹변수로 우선 정렬한다.
arrange(by_tailnum, day)
arrange(flights, day)
## mutate() 와 filter()는 윈도우 함수(window function)과 더불어 사용한다.
# vignette("window-functions")
# 각 비행기 별로 가장 큰 air_time을 가진 항공편을 추출하는 명령
filter(by_tailnum, min_rank(desc(air_time)) == 1)
# 각 비행기 별로 air_time의 크기를 작은 순서대로 순위화하는 명령
mutate(by_tailnum, rank = min_rank(air_time))
## sample_n() 과 sample_frac() 은 그룹별로 랜덤샘플한다.
sample_n(by_tailnum, 1)
## slice() 는 각 그룹별로 행을 추출한다.
#slice는 그룹별 원소의 개수보다 많이 숫자를 설정하면 개수만큼 나온다.
slice(by_tailnum, 5)
# summarise() 는 각 그룹별로 요약통계량을 계산한다.
# 아래 참고
#(by_tailnum = group_by(flights, tailnum))
delay = summarise(by_tailnum,
count = n(),
dist = mean(distance, na.rm = TRUE),
delay = mean(arr_delay, na.rm = TRUE))
delay = filter(delay, count > 20, dist < 2000)
delay
### exercise 26/43
groups(by_tailnum)
by_tailnum = ungroup(by_tailnum)
groups(by_tailnum)
### exercise 26/43
rowwise(by_tailnum)
delay = summarise(rowwise(by_tailnum),
count = n(),
dist = mean(distance, na.rm = TRUE),
delay = mean(arr_delay, na.rm = TRUE))
delay = filter(delay, dist < 2000)
delay
########################### ggplot2 패키지를 이용한 시각화 ###########################
by_tailnum = group_by(flights, tailnum)
delay = summarise(by_tailnum,
count = n(),
dist = mean(distance, na.rm = TRUE),
delay = mean(arr_delay, na.rm = TRUE))
delay = filter(delay, count > 20, dist < 2000)
delay
ggplot(delay, aes(dist, delay)) +
geom_point(aes(size = count), alpha = 1/2) +
geom_smooth() +
scale_size_area()
sum(select(delay, count)>=1000)
sum(is.na(select(delay, delay)))
########################### %>% 연산자와 체이닝(chaining) ###########################
# The downside of the functional nature of dplyr is that when you combine multiple data manipulation
# operations, you have to read from the inside out and the arguments may be very distant to the function
# call. These functions providing an alternative way of calling dplyr (and other data manipulation)
# functions that you read can from left to right. (dyplyr 패키지 설명의 chain 부분)
# 임시변수 방법
a1 = group_by(flights, year, month, day)
a2 = select(a1, arr_delay, dep_delay)
a3 = summarise(a2,
arr = mean(arr_delay, na.rm = TRUE),
dep = mean(dep_delay, na.rm = TRUE))
(a4 = filter(a3, arr > 30 | dep > 30))
# 중첩 방법
filter(
summarise(
select(
group_by(flights, year, month, day),
arr_delay, dep_delay
),
arr = mean(arr_delay, na.rm = TRUE),
dep = mean(dep_delay, na.rm = TRUE)
),
arr > 30 | dep > 30
)
# 두 방법 모두 가독성이 떨어진다.
### %>% 연산자
# simple1
head(iris)
iris %>% head
# simple2
head(iris, 10)
iris %>% head(10)
# 위의 임시변수 방법과 중첩 방법과 결과가 같다.
# 또한 가독성이 좋다.
flights %>%
group_by(year, month, day) %>%
select(arr_delay, dep_delay) %>%
summarise(
arr = mean(arr_delay, na.rm = TRUE),
dep = mean(dep_delay, na.rm = TRUE)
) %>%
filter(arr > 30 | dep > 30)
################################## 연습문제(31/43) ##################################
# 날짜 형식을 다루는 패키지
library(lubridate)
# C 언어 중에서 지역화를 지정하는 기능.
# 구체적으로는 setlocale 함수나 localconv 함수로 통화(通貨)의 표현,
# 10진수 표현(소수 점이 점인가 콤마인가 등), 일시(日時) 표현의 차이 등 각국 특유의 표현을 지정한다.
# [네이버 지식백과] 로컬 [locale] (컴퓨터인터넷IT용어대사전, 2011. 1. 20., 일진사)
Sys.setlocale("LC_TIME", "usa")
(per_day = flights %>%
group_by(year, month, day) %>%
summarise(flights = n(),
distance = mean(distance),
arr = mean(arr_delay, na.rm = TRUE),
dep = mean(dep_delay, na.rm = TRUE)
) %>%
mutate(dt_str = sprintf("%04d-%02d-%02d", year, month, day),
dt = parse_date_time(dt_str, "ymd",tz = "US/Eastern"),
dow = wday(dt, label=TRUE)
)
)
# 각각의 도착지에 대해 비행기 개수
flights %>% group_by(tailnum, dest) %>% distinct() %>% group_by(dest) %>% summarise(count = n())
# 각각의 도착지에 대해 다른 항공편 편수
flights %>% group_by(dest) %>% summarise(count = n())
# 각 날짜별 비행기 편수
select(per_day, year:flights)
flights %>% group_by(year,month, day) %>% summarise(count = n())
# 각 요일별 비행기 편수
ungroup(per_day) %>% group_by(dow) %>% summarise(count = sum(flights))
ungroup(per_day) %>% group_by(dow) %>% summarise(count = mean(flights))
per_day %>% ggplot(aes(dt, flights)) +
geom_line() + geom_point() + geom_smooth()
per_day %>% ggplot(aes(dow, flights)) +
geom_boxplot()
per_day %>% ggplot(aes(dt, dep)) +
geom_line() + geom_point() + geom_smooth()
per_day %>% ggplot(aes(dow, dep)) +
geom_boxplot()
per_day %>% ggplot(aes(dt, distance)) +
geom_line() + geom_point() + geom_smooth()
per_day %>% ggplot(aes(dow, distance)) +
geom_boxplot()
########################## join 함수를 이용한 테이블 결합 ##########################
(df1 <- data_frame(x = c(1, 2), y = 2:1))
(df2 <- data_frame(x = c(1, 3), a = 10, b = "a"))
# inner_join() 함수는 x와 y에 모두 매칭되는 행만 포함한다. 교집합이다.
df1 %>% inner_join(df2)
# left_join() 함수는 x 의 모든 행을 포함한다. 매칭되지 않은 y 테이블 변수들을 NA가 된다. 차집합이다.
df1 %>% left_join(df2)
# right_join() 함수는 y 테이블의 모든 행을 포함한다.
# left_join(y,x) 의 경우와 행의 순서만 다를 뿐 동일한 결과를 준다.
df1 %>% right_join(df2)
# full_join() 함수는 x와 y의 모든 행을 포함한다. 합집합이다.
df1 %>% full_join(df2)
################################## 연습문제(35/43) ##################################
(flights2 = flights %>% select(year:day, hour, origin, dest, tailnum, carrier))
airlines
flights2 %>% left_join(airlines)
# flights2 %>% left_join(airlines, by="carrier")
# rm(list = ls())
################################# 연습문제(37/43)-1 #################################
ds <- tbl_df(mtcars)
ds
as.data.frame(ds)
if (require("Lahman") && packageVersion("Lahman") >= "3.0.1") {
batting <- tbl_df(Batting)
dim(batting)
colnames(batting)
head(batting)
# Data manipulation verbs ---------------------------------------------------
filter(batting, yearID > 2005, G > 130)
select(batting, playerID:lgID)
arrange(batting, playerID, desc(yearID))
summarise(batting, G = mean(G), n = n())
mutate(batting, rbi2 = if(is.null(AB)) 1.0 * R / AB else 0)
# Group by operations -------------------------------------------------------
# To perform operations by group, create a grouped object with group_by
players <- group_by(batting, playerID)
head(group_size(players), 100)
summarise(players, mean_g = mean(G), best_ab = max(AB))
best_year <- filter(players, AB == max(AB) | G == max(G))
progress <- mutate(players, cyear = yearID - min(yearID) + 1,
rank(desc(AB)), cumsum(AB))
# When you group by multiple level, each summarise peels off one level
per_year <- group_by(batting, playerID, yearID)
stints <- summarise(per_year, stints = max(stint))
filter(stints, stints > 3)
summarise(stints, max(stints))
mutate(stints, cumsum(stints))
# Joins ---------------------------------------------------------------------
player_info <- select(tbl_df(Master), playerID, birthYear)
hof <- select(filter(tbl_df(HallOfFame), inducted == "Y"),
playerID, votedBy, category)
# Match players and their hall of fame data
inner_join(player_info, hof)
# Keep all players, match hof data where available
left_join(player_info, hof)
# Find only players in hof
semi_join(player_info, hof)
# Find players not in hof
anti_join(player_info, hof)
}
# rm(list = ls())
################################# 연습문제(37/43)-2 #################################
library(babynames)
babynames
class(babynames)
babynames %>% group_by(year) %>% filter(min_rank(desc(n)) == 1)
babynames %>% group_by(year) %>% filter(min_rank(desc(n)) <= 10) %>% summarise(rate = sum(prop))
# rm(list = ls())
################################# 연습문제(37/43)-3 #################################
# 주소로 지정할 수도 있지만...
batting = read.csv(file.choose())
master = read.csv(file.choose())
head(batting)
head(master)
batting = tbl_df(batting)
master = tbl_df(master)
class(batting)
class(master)
# 각 연도별 타자들의 타율을 계산하라. ( AVG(타율) = H(안타) / AB(타수) )
mutate(batting, AVG = H / AB) %>% select(playerID:lgID, AVG)
full.batting = mutate(batting, AVG = H / AB)
# 각 연도별로 선수들의 타율의 분포를 시각화하라.
### plot
full.batting %>% mutate(year = factor(yearID)) %>%
group_by(year) %>%
ggplot(aes(year, AVG)) + geom_point() + geom_smooth(aes(group=1))
# 자료를 선별하여 분포를 그린 것이므로 참고만 하자.
full.batting %>% mutate(year = factor(yearID)) %>%
filter(AB >= 100, yearID >= 1985) %>%
group_by(year) %>%
ggplot(aes(year, AVG)) + geom_point() + geom_smooth(aes(group=1))
full.batting %>% filter(AB >= 100) %>%
group_by(yearID) %>%
summarise( n = n()) %>%
ggplot(aes(yearID, n))+ geom_point()
### box-plot
full.batting %>% mutate(year = factor(yearID)) %>%
group_by(year) %>%
ggplot(aes(year, AVG)) + geom_boxplot()
# 자료를 선별하여 분포를 그린 것이므로 참고만 하자.
full.batting %>% mutate(year = factor(yearID)) %>%
filter(AB >= 100, yearID >= 1985) %>%
group_by(year) %>%
ggplot(aes(year, AVG)) + geom_boxplot()
# 1980년 이후 최고타율을 기록한 선수들의 리스트를 구하라.
(AVGtable = full.batting %>% select(playerID:lgID, AVG))
# min_rank() 함수의 성질
cw <- c(1,2,3,4,NA, 2913)
min_rank(cw)
temp = AVGtable %>% group_by(playerID) %>% filter(min_rank(desc(AVG)) == 1 & yearID >= 1980)
AVGtable %>% filter(playerID == "blackti01")
# 결과 테이블을 Master 테이블과 join 해서 선수의 이름을 표시하라.
(name = master %>% select(playerID, nameFirst, nameLast, nameGiven, retroID, bbrefID))
# select() 순서를 바꾸기 위해 사용해보자.
list = temp %>% left_join(name, by="playerID") %>% select(playerID, nameFirst:bbrefID, yearID:AVG)
list
#################################################################################################
########## tidyr ##########
#################################################################################################
# Key concepts : Long format <-> Wide format, Pivot table
# Core contents : spread() and gather() function
# library(tidyr)
##################################### 1. spread() #####################################
# 설명 : Long format -> Wide format
# 기초적인 사용법 : spread(주어진 자료(데이터프레임), 피벗이 되는 열, 대상이 되는 열)
digest = flights %>% group_by(carrier, origin) %>%
summarise(air_time_mean = mean(air_time, na.rm=TRUE))
as.data.frame(digest)
# 피벗 테이블
spread(digest, origin, air_time_mean)
flights %>% group_by(carrier, origin) %>%
summarise(air_time_mean = mean(air_time, na.rm=TRUE)) %>%
spread(origin, air_time_mean)
# NA를 0으로 채우기
spread(digest, origin, air_time_mean, fill = 0)
##################################### 2. gather() #####################################
# 설명 : Wide format -> Long format
# 기초적인 사용법 : gather(주어진 자료(데이터프레임), 변수명을 포함한 열, 대상이 되는 열,
# long format 으로 전환할 열)
series = flights %>% group_by(month, day) %>%
summarise(air_time_mean = mean(air_time, na.rm=TRUE), dep_delay_mean = mean(dep_delay, na.rm=TRUE))
series
series %>% gather(indicators, value, air_time_mean, dep_delay_mean)
series %>% gather(indicators, value, -month, -day)
2일 : 데이터 시각화 (Data Visualization)
# Install Packages
install.packages('ggmap')
install.packages('ggplot2')
install.packages('tm')
install.packages('SnowballC')
install.packages('wordcloud')
install.packages('networkD3')
install.packages('curl')
install.packages('GGally')
install.packages('ggvis')
install.packages('shiny')
install.packages('plotly')
# 1차원 자료용 그래픽
# Data: 초미세먼지 자료
setwd("~/Desktop/개인/Big Data 시각화 (2016.02)")
pollution <- read.csv("avgpm25.csv", colClasses = c("numeric", "character", "factor", "numeric", "numeric"))
head(pollution)
## pm25 fips region longitude latitude
## 1 9.771185 01003 east -87.74826 30.59278
## 2 9.993817 01027 east -85.84286 33.26581
## 3 10.688618 01033 east -87.72596 34.73148
## 4 11.337424 01049 east -85.79892 34.45913
## 5 12.119764 01055 east -86.03212 34.01860
## 6 10.827805 01069 east -85.35039 31.18973
# Five Number Summary
summary(pollution$pm25)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.383 8.549 10.050 9.836 11.360 18.440
# Boxplot
boxplot(pollution$pm25, col="blue")
boxplot(pollution$pm25, col="blue")
abline(h=12)
# Histogram
hist(pollution$pm25, col="green")
rug(pollution$pm25)
hist(pollution$pm25, col="green",breaks=100)
rug(pollution$pm25)
hist(pollution$pm25, col='green', breaks=c(2,5,8,10,14,18,20))
# Histogram with bar
hist(pollution$pm25, col="green")
abline(v=12,lwd=2)
abline(v=median(pollution$pm25), col="magenta", lwd=4)
abline(h=60, col="skyblue", lwd=4)
# Barplot
barplot(table(pollution$region), col = "wheat", main = "Number of Counties in Each Region")
# n차원 자료용 그래픽
# Multiple Boxplot
boxplot(pm25 ~ region, data = pollution, col = "red") #formula in R; y~x1+x2+...
# Multiple Histograms
par(mfrow = c(2, 1), mar = c(4, 4, 2, 1)) ### par(mfrow = c(1, 1)) or dev.off() ; 초기화
hist(subset(pollution, region == "east")$pm25, col = "green")
hist(subset(pollution, region == "west")$pm25, col = "green")
# Scatterplot
with(pollution, plot(latitude, pm25)) ## plot(pm25 ~ latitude, data=pollution) 과 같은 표현!
abline(h=12,lwd=2,lty=2)
with(pollution, plot(latitude, pm25, col = region))
abline(h=12,lwd=2,lty=2)
# Multiple Scatterplots
par(mfrow = c(1, 2), mar = c(5, 4, 2, 1)) ## margin : c(bottom, left, top, right)
with(subset(pollution, region == "west"), plot(latitude, pm25, main = "West"))
with(subset(pollution, region == "east"), plot(latitude, pm25, main = "East"))
# R Plot System
# Base plotting system
# Data loading : airquality
library(datasets)
with(airquality, plot(Wind, Ozone))
# Base Plot with Annotation
with(airquality, plot(Wind, Ozone))
title(main="Ozone and Wind in New Yorm City") ## Add a title
# Base Plot with Regression Line
with(airquality, plot(Wind, Ozone, main = "Ozone and Wind in New York City", pch = 20))
?points
model <- lm(Ozone ~ Wind, airquality)
abline(model, lwd=2)
abline(1,2, lwd=2) ### intercept(y절편), slope(기울기), thickness(두께)
# Multiple Base Plots
par(mfrow = c(1, 3), mar = c(4, 4, 2, 1), oma = c(0, 0, 2, 0))
with (airquality, {
plot(Wind, Ozone, main = "Ozone and Wind")
plot(Solar.R, Ozone, main = "Ozone and Solar Radiation")
plot(Temp, Ozone, main = "Ozone and Temperature")
mtext("Ozone and Weather in New York City", outer = TRUE) #main title;;; 통합 제목
}
)
# Lattice System in R
# Lattice system in R
library(lattice)
str(airquality)
airquality <- transform(airquality, Month = factor(Month)) # Convert 'Month' to a factor variable
levels(airquality$Month)
xyplot(Ozone ~ Wind | Month, data = airquality, layout = c(5, 1))
xyplot(Ozone ~ Wind, data = airquality)
# Lattice, Why?
set.seed(10)
x <- rnorm(100)
f <- rep(0:1, each = 50)
y<-x+f-f*x+rnorm(100,sd=0.5)
f <- factor(f, labels = c("Group 1", "Group 2")) # group1(f=0; y=x+norm)=1:50, group2(f=1; y=1+norm)=51:100
xyplot(y ~ x)
xyplot(y ~ x | f, layout = c(2, 1))
# Lattice Panel Functions: Regression line
xyplot(y ~ x | f, panel = function(x, y, ...) {
panel.xyplot(x, y, ...) # First call default panel function
panel.lmline(x, y, col = 2) # Overlay a simple linear regression line
})
ftn = function(x, y, ...) {
panel.xyplot(x, y, ...) # First call default panel function
panel.lmline(x, y, col = 2) # Overlay a simple linear regression line
}
xyplot(y~x | f, panel=ftn)
# ggplot2
# Data ; MAACS
library("ggplot2")
load("maacs.Rda")
# basic ggplot ; using qplot (quick plot)
qplot(logpm25, NocturnalSympt, data = maacs, facets = . ~ bmicat,
geom = c("point"))+stat_smooth(method="lm")
# Function ggplot
g<-ggplot(maacs, aes(logpm25, NocturnalSympt))
g
g+geom_point()
g+geom_point()+geom_smooth(method = "lm")
g+geom_point()+geom_smooth(method = "lm")+facet_grid(.~bmicat)
# Modifying Aesthetics
g + geom_point(color = "steelblue", size=4, alpha = 1/2)
g + geom_point(aes(color = bmicat), size=4, alpha=1/2)
# Modifying Labels
g + geom_point(aes(color = bmicat), size=4, alpha=1/2) +
labs(title = "MAACS Cohort") + labs(x=expression("log"*PM[2.5]),y="Nocturnal Symptoms")
# Customizing the Smooth
g + geom_point(aes(color = bmicat), size=4, alpha=1/2) +
geom_smooth(size = 4, linetype = 3, method = "lm", se = FALSE)
# Parallel Coordinate Plot
require(GGally)
head(mtcars)
ggparcoord(mtcars, columns = c(1,5:10))+geom_line()
# ggmap의 활용(지도 기반 시각화)
# ggmap 패키지
library(ggmap)
geocode("Seoul National University")
## lon lat
## 1 126.978 37.56654
korea <- c(left = 124, bottom = 32, right = 132, top = 43.5)
map <- get_stamenmap(korea, zoom = 5, maptype = "toner-lite")
ggmap(map)
# 휴스턴(Houston) 도심의 범죄 자료
# Basic plot on the map
# crime Data from Google Maps (included in ggmap packages)
str(crime)
## 'data.frame': 86314 obs. of 17 variables:
## $ time : POSIXt, format: "2010-01-01 15:00:00" "2010-01-01 15:00:00" ...
## $ date : chr "1/1/2010" "1/1/2010" "1/1/2010" "1/1/2010" ...
## $ hour : int 0 0 0 0 0 0 0 0 0 0 ...
## $ premise : chr "18A" "13R" "20R" "20R" ...
## $ offense : Factor w/ 7 levels "aggravated assault",..: 4 6 1 1 1 3 3 3 3 3 ...
## $ beat : chr "15E30" "13D10" "16E20" "2A30" ...
## $ block : chr "9600-9699" "4700-4799" "5000-5099" "1000-1099" ...
## $ street : chr "marlive" "telephone" "wickview" "ashland" ...
## $ type : chr "ln" "rd" "ln" "st" ...
## $ suffix : chr "-" "-" "-" "-" ...
## $ number : int 1 1 1 1 1 1 1 1 1 1 ...
## $ month : Ord.factor w/ 8 levels "january"<"february"<..: 1 1 1 1 1 1 1 1 1 1 ...
## $ day : Ord.factor w/ 7 levels "monday"<"tuesday"<..: 5 5 5 5 5 5 5 5 5 5 ...
## $ location: chr "apartment parking lot" "road / street / sidewalk" "residence / house" "residence / house" ...
## $ address : chr "9650 marlive ln" "4750 telephone rd" "5050 wickview ln" "1050 ashland st" ...
## $ lon : num -95.4 -95.3 -95.5 -95.4 -95.4 ...
## $ lat : num 29.7 29.7 29.6 29.8 29.7 ...
downtown<-subset(crime,
-95.39681 <= lon & lon <= -95.34188 &
29.73631 <= lat & lat <= 29.78400)
qmplot(lon, lat, data = downtown, maptype = "toner-background", color = I("red"))
Different Source / Type of map
배경에 표시될 지도의 출처(source option)와 종류(maptype option)를 다양하게 바꿈으로써,
상황에 어울리는 시각화가 가능하다.
source = “stamen”(기본설정)일 때, 가능한 maptype :
“terrain”, “terrain-background”, “terrain-labels”, “terrain-lines”, “toner”,
“toner-2010”, “toner-2011”, “toner-background”, “toner-hybrid”, “toner-labels”,
“toner-lines”, “toner-lite”, “watercolor”
source = “google”일 때, 가능한 maptype :
“terrain”, “satellite”, “roadmap”, and “hybrid”
qmplot(lon, lat, data = downtown, maptype = "watercolor", color = I("red"))
qmplot(lon, lat, data = downtown, maptype = "terrain-background", color = I("red"))
qmplot(lon, lat, data = downtown, maptype = "toner-hybrid", color = I("red"))
qmplot(lon, lat, data = downtown, source = "google", maptype = "terrain", color = I("red"))
qmplot(lon, lat, data = downtown, source = "google", maptype = "satellite", color = I("red"))
qmplot(lon, lat, data = downtown, source = "google", maptype = "hybrid", color = I("red"))
# 범죄율의 밀도(Density) 표현
robberies <- subset(downtown, offense == "robbery")
qmplot(lon, lat, data = robberies, geom = "blank", zoom = 15, maptype = "toner-background", darken = .7) +
stat_density_2d(aes(fill = ..level..), geom = "polygon", alpha = .3, color = NA)+
scale_fill_gradient2("Robbery\nPropensity", low = "white", mid = "yellow", high = "red", midpoint = 700)
# 범죄율의 밀도(Density) 표현 - ggmap 함수를 이용한 표현
downtownmap<-get_stamenmap(c(left=-95.39681,bottom= 29.73631,right=-95.34188,top=29.78400),zoom = 15, maptype = "toner-background", darken= .7)
ggmap(downtownmap, darken=.7)+
stat_density_2d(aes(fill = ..level..), geom = "polygon", alpha = .3, color = NA,data=robberies)+
scale_fill_gradient2("Robbery\nPropensity", low = "white", mid = "yellow", high = "red", midpoint = 700)
# 범죄의 종류별 시각화
qmplot(lon, lat, data = downtown, maptype = "toner-background", color = offense) +
facet_wrap(~ offense)
# Quiz1
# baseball.csv 파일에는 한국 프로야구 10개의 팀의 연고지(구장)의 위치의 경도(lon)와 위도(lat) 값과 2015년 프로야구 관중수(crowd)가 담겨있다.
# 파일을 불러들여서(header = TRUE), 다음과 같이 지도에 각 구단의 연고지를 구단 별로 다른 색으로 위치를 표시해 보자.
# Hint :
# qmplot(lon, lat, data = downtown, maptype = “toner-background”, color = I(“red”)) 이 코드를 알맞게 변형해보세요.
# Quiz1 - Advanced
# Quiz의 결과는 다소 눈에 익지 않다. 가장 큰 이유중 하나는 대한민국 지도가 한꺼번에 보이지 않기 때문일 것이며, 10개의 팀을 모두 각기 다른 색으로 표현하는 것은 10개의 팀을 쉽게 구분할 수 있는 방법은 아닐 수 있다. 또한, baseball.csv 파일안에는 프로야구 관중수(crowd)라는 정보도 포함되어 있다.
# ggmap과 ggplot의 함수들을 함께 사용하면 다음과 같이 지도 영역의 자유로운 설정과 다양한 정보를 담을 수 있다.
# 텍스트 자료의 시각화
library(tm)
library(SnowballC)
library(wordcloud)
jeopQ <- read.csv('JEOPARDY_CSV.csv', stringsAsFactors = FALSE)
###Data Cleaning
jeopCorpus <- Corpus(VectorSource(jeopQ$Question))
jeopCorpus <- tm_map(jeopCorpus, removePunctuation) ## . , ? ! ()
jeopCorpus <- tm_map(jeopCorpus, removeWords, c("the", "this", "The", "This", stopwords('english'))) ## stopwords("en")
jeopCorpus <- tm_map(jeopCorpus, stemDocument) ## went -> go
jeopCorpus <- tm_map(jeopCorpus, PlainTextDocument)
wordcloud(jeopCorpus, max.words = 100, random.order = FALSE)
# 네트워크 자료의 시각화
# Basic Network Diagram
library(networkD3)
src <- c("A", "A", "A", "A",
"B", "B", "C", "C", "D")
target <- c("B", "C", "D", "J",
"E", "F", "G", "H", "I")
networkData <- data.frame(src, target)
simpleNetwork(networkData)
# Sankey Diagram
library(networkD3)
URL <- paste0(
"https://cdn.rawgit.com/christophergandrud/networkD3/",
"master/JSONdata/energy.json")
Energy <- jsonlite::fromJSON(URL)
head(Energy$nodes)
## name
## 1 Agricultural 'waste'
## 2 Bio-conversion
## 3 Liquid
## 4 Losses
## 5 Solid
## 6 Gas
head(Energy$links)
## source target value
## 1 0 1 124.729
## 2 1 2 0.597
## 3 1 3 26.862
## 4 1 4 280.322
## 5 1 5 81.144
## 6 6 2 35.000
sankeyNetwork(Links = Energy$links, Nodes = Energy$nodes, Source = "source",
Target = "target", Value = "value", NodeID = "name",
units = "TWh", fontSize = 12, nodeWidth = 30)
# Interactive Graphics
# ggvis system
# Data : mtcars data
# ggvis
# The %>% operator passes the left hand object to the first argument of the right hand expression
library(ggvis)
library(shiny)
p <- ggvis(mtcars, x = ~wt, y = ~mpg)
layer_points(p)
# ggvis using pipe function “%>%”
mtcars %>% ## DATA
ggvis(x = ~wt, y = ~mpg) %>% ## model
layer_points() ## plotting
# Different Layers
mtcars %>%
ggvis(x = ~wt, y = ~mpg) %>%
layer_points()
mtcars %>%
ggvis(x = ~wt, y = ~mpg) %>%
layer_lines()
mtcars %>%
ggvis(x = ~wt, y = ~mpg) %>%
layer_bars()
mtcars %>%
ggvis(x = ~wt, y = ~mpg) %>%
layer_smooths()
# Interactive Controls
mtcars %>%
ggvis(~wt, ~mpg,
size := input_slider(10, 100), # := set a property to a value
#size = ~wt, # = set a property w.r.t a variable with ~variable format
opacity := input_slider(0, 1)
) %>%
layer_points() ##Press Escape/Ctrl + C to stop.
mtcars %>%
ggvis(~wt, ~mpg,
size := input_numeric(30),
opacity := input_radiobuttons(choices=c(0.2,0.4,0.6,0.8,1))
) %>%
layer_points() ##Press Escape/Ctrl + C to stop.
# Examples
mtcars %>% ggvis(x=~wt) %>%
layer_densities(
adjust = input_slider(.1,2, value =1, step=.1, label="Bandwidth adjustment"),
kernel = input_select(
c("Gaussian"="gaussian","Epanechnikov"="epanechnikov","Rectangular"="rectangular", "Triangular"="triangular", "Biweight"="biweight","Cosine"="cosine", "Optcosine"="optcosine"), label="Kernel"))
# Quiz2
# Basic R plot 을 통해 “pollution” 데이터의 scatter plot을 다음과 같이 그려볼 수 있었다.
with(pollution, plot(latitude, pm25))
# 이번에는 같은 그래프를 “ggvis”패키지와 pipe function “%>%”을 사용하여
# 다음과 같이 투명도(opacity), 크기(size), 색깔(fill)을 interactive하게 바꿀 수 있는 그래프를 그려보자.
# Hint
# 투명도(“opacity” 옵션) : input_slider 이용
# 크기(“size” 옵션) : input_numeric 이용
# 색깔(“fill” 옵션) : input_select 이용
# plotly system
# plotly
library(plotly)
set.seed(100)
d<-diamonds[sample(nrow(diamonds),1000), ]
plot_ly(d, x=carat, y=price, text=
paste("Clarity: ", clarity),
mode = "markers", color = carat, size= carat)
plot_ly(ggplot2::diamonds, y=price, color=cut, type="box")
plot_ly(z=volcano, type="surface")
3일 : 회귀 분석 (Regression Analysis)
install.packages("ISLR")
# 선형회귀모형의 적합
Advertising = read.csv("~/Desktop/Advertising.csv", header = TRUE)
lm.fit = lm(Sales ~ TV, data = Advertising)
summary(lm.fit)
##
## Call:
## lm(formula = Sales ~ TV, data = Advertising)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.3860 -1.9545 -0.1913 2.0671 7.2124
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.032594 0.457843 15.36 <2e-16 ***
## TV 0.047537 0.002691 17.67 <2e-16 ***
## ---
## Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
##
## Residual standard error: 3.259 on 198 degrees of freedom
## Multiple R-squared: 0.6119,Adjusted R-squared: 0.6099
## F-statistic: 312.1 on 1 and 198 DF, p-value: < 2.2e-16
confint(lm.fit, level=0.95)
## 2.5 % 97.5 %
## (Intercept) 6.12971927 7.93546783
## TV 0.04223072 0.05284256
# 중회귀모형의 적합도
lm.fit = lm(Sales ~ TV + Radio + Newspaper, data = Advertising)
summary(lm.fit)
##
## Call:
## lm(formula = Sales ~ TV + Radio + Newspaper, data = Advertising)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.8277 -0.8908 0.2418 1.1893 2.8292
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.938889 0.311908 9.422 <2e-16 ***
## TV 0.045765 0.001395 32.809 <2e-16 ***
## Radio 0.188530 0.008611 21.893 <2e-16 ***
## Newspaper -0.001037 0.005871 -0.177 0.86
## ---
## Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
##
## Residual standard error: 1.686 on 196 degrees of freedom
## Multiple R-squared: 0.8972,Adjusted R-squared: 0.8956
## F-statistic: 570.3 on 3 and 196 DF, p-value: < 2.2e-16
cor(Advertising)
## X TV Radio Newspaper Sales
## X 1.00000000 0.01771469 -0.11068044 -0.15494414 -0.05161625
## TV 0.01771469 1.00000000 0.05480866 0.05664787 0.78222442
## Radio -0.11068044 0.05480866 1.00000000 0.35410375 0.57622257
## Newspaper -0.15494414 0.05664787 0.35410375 1.00000000 0.22829903
## Sales -0.05161625 0.78222442 0.57622257 0.22829903 1.00000000
# 가변수 (Dummy variable)
Credit = read.csv("~/Desktop/Credit.csv", header = TRUE)
lm.fit = lm(Balance ~ Gender, data = Credit)
summary(lm.fit)
##
## Call:
## lm(formula = Balance ~ Gender, data = Credit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -529.54 -455.35 -60.17 334.71 1489.20
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 509.80 33.13 15.389 <2e-16 ***
## GenderFemale 19.73 46.05 0.429 0.669
## ---
## Residual standard error: 460.2 on 398 degrees of freedom
## Multiple R-squared: 0.0004611,Adjusted R-squared: -0.00205
## F-statistic: 0.1836 on 1 and 398 DF, p-value: 0.6685
# 가변수를 포함한 회귀모형
lm.fit = lm(Balance~ Ethnicity, data=Credit)
summary(lm.fit)
##
## Call:
## lm(formula = Balance ~ Ethnicity, data = Credit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -531.00 -457.08 -63.25 339.25 1480.50
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 531.00 46.32 11.464 <2e-16 ***
## EthnicityAsian -18.69 65.02 -0.287 0.774
## EthnicityCaucasian -12.50 56.68 -0.221 0.826
## ---
## Residual standard error: 460.9 on 397 degrees of freedom
## Multiple R-squared: 0.0002188,Adjusted R-squared: -0.004818
## F-statistic: 0.04344 on 2 and 397 DF, p-value: 0.9575
# 가법성 가정 완화 (교호작용)
lm.fit = lm(Sales ~ TV*Radio, data = Advertising)
summary(lm.fit)
##
## Call:
## lm(formula = Sales ~ TV * Radio, data = Advertising)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.3366 -0.4028 0.1831 0.5948 1.5246
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.750e+00 2.479e-01 27.233 <2e-16 ***
## TV 1.910e-02 1.504e-03 12.699 <2e-16 ***
## Radio 2.886e-02 8.905e-03 3.241 0.0014 **
## TV:Radio 1.086e-03 5.242e-05 20.727 <2e-16 ***
## ---
## Residual standard error: 0.9435 on 196 degrees of freedom
## Multiple R-squared: 0.9678,Adjusted R-squared: 0.9673
## F-statistic: 1963 on 3 and 196 DF, p-value: < 2.2e-16
# 변수선택
library(ISLR)
names(Hitters)
## [1] "AtBat" "Hits" "HmRun" "Runs" "RBI"
## [6] "Walks" "Years" "CAtBat" "CHits" "CHmRun"
## [11] "CRuns" "CRBI" "CWalks" "League" "Division"
## [16] "PutOuts" "Assists" "Errors" "Salary" "NewLeague"
dim(Hitters)
## [1] 322 20
sum(is.na(Hitters$Salary))
## [1] 59
Hitters = na.omit(Hitters)
dim(Hitters)
## [1] 263 20
sum(is.na(Hitters))
## [1] 0
# RSS를 포함한 변수선택을 위한 통계량 변화
par(mfrow = c(2,2))
plot(reg.summary$rss, xlab = "Number of Variables", ylab = "RSS", type = "l")
plot(reg.summary$adjr2, xlab = "Number of Variables", ylab = "Adjusted RSq", type = "l")
which.max(reg.summary$adjr2)
## [1] 11
points(11, reg.summary$adjr2[11], col = "red", cex = 2, pch = 20)
plot(reg.summary$cp, xlab = "Number of Variables", ylab = "Cp", type = "l")
which.min(reg.summary$cp)
## [1] 10
points(10, reg.summary$cp[10], col = "red", cex = 2, pch = 20)
plot(reg.summary$bic, xlab = "Number of Variables", ylab = "BIC", type = "l")
which.min(reg.summary$bic)
## [1] 6
points(6, reg.summary$bic[6], col = "red", pch = 20)
# 능선회귀와 라쏘 (Ridge Regression and LASSO)
x = model.matrix(Salary ~ . , data = Hitters)[,-1]
y = Hitters$Salary
library(glmnet)
grid = 10^seq(10, -2, length=100)
ridge.mod = glmnet(x, y, alpha=0, lambda = grid)
dim(coef(ridge.mod))
## [1] 20 100
ridge.mod$lambda[50]
## [1] 11497.57
coef(ridge.mod)[,50]
## (Intercept) AtBat Hits HmRun Runs
## 407.356050200 0.036957182 0.138180344 0.524629976 0.230701523
## RBI Walks Years CAtBat CHits
## 0.239841459 0.289618741 1.107702929 0.003131815 0.011653637
## CHmRun CRuns CRBI CWalks LeagueN
## 0.087545670 0.023379882 0.024138320 0.025015421 0.085028114
## DivisionW PutOuts Assists Errors NewLeagueN
## -6.215440973 0.016482577 0.002612988 -0.020502690 0.301433531
sqrt(sum(coef(ridge.mod)[-1,50]^2))
## [1] 6.360612
ridge.mod$lambda[60]
## [1] 705.4802
coef(ridge.mod)[,60]
## (Intercept) AtBat Hits HmRun Runs
## 54.32519950 0.11211115 0.65622409 1.17980910 0.93769713
## RBI Walks Years CAtBat CHits
## 0.84718546 1.31987948 2.59640425 0.01083413 0.04674557
## CHmRun CRuns CRBI CWalks LeagueN
## 0.33777318 0.09355528 0.09780402 0.07189612 13.68370191
## DivisionW PutOuts Assists Errors NewLeagueN
## -54.65877750 0.11852289 0.01606037 -0.70358655 8.61181213
sqrt(sum(coef(ridge.mod)[-1,60]^2))
## [1] 57.11001
predict(ridge.mod, s=50, type = "coefficients")[1:20,]
## (Intercept) AtBat Hits HmRun Runs
## 4.876610e+01 -3.580999e-01 1.969359e+00 -1.278248e+00 1.145892e+00
## RBI Walks Years CAtBat CHits
## 8.038292e-01 2.716186e+00 -6.218319e+00 5.447837e-03 1.064895e-01
## CHmRun CRuns CRBI CWalks LeagueN
## 6.244860e-01 2.214985e-01 2.186914e-01 -1.500245e-01 4.592589e+01
## DivisionW PutOuts Assists Errors NewLeagueN
## -1.182011e+02 2.502322e-01 1.215665e-01 -3.278600e+00 -9.496680e+00
set.seed(1)
train = sample(1:nrow(x), nrow(x)/2)
test = (-train)
y.test = y[test]
ridge.mod = glmnet(x[train, ], y[train], alpha = 0, lambda = grid, thresh = 1e-12)
ridge.pred = predict(ridge.mod, s=4, newx = x[test,])
mean((ridge.pred - y.test)^2)
## [1] 101036.8
mean((mean(y[train]) - y.test)^2)
## [1] 193253.1
ridge.pred = predict(ridge.mod, s=1e10, newx = x[test,])
mean((ridge.pred - y.test)^2)
## [1] 193253.1
set.seed(1)
cv.out = cv.glmnet(x[train,], y[train], alpha = 0)
plot(cv.out)
bestlam = cv.out$lambda.min
bestlam
## [1] 211.7416
ridge.pred = predict(ridge.mod, s=bestlam, newx = x[test,])
mean((ridge.pred - y.test)^2)
## [1] 96015.51
out = glmnet(x, y, alpha = 0)
predict(out, type = "coefficients", s = bestlam)[1:20,]
## (Intercept) AtBat Hits HmRun Runs
## 9.88487157 0.03143991 1.00882875 0.13927624 1.11320781
## RBI Walks Years CAtBat CHits
## 0.87318990 1.80410229 0.13074381 0.01113978 0.06489843
## CHmRun CRuns CRBI CWalks LeagueN
## 0.45158546 0.12900049 0.13737712 0.02908572 27.18227535
## DivisionW PutOuts Assists Errors NewLeagueN
## -91.63411299 0.19149252 0.04254536 -1.81244470 7.21208390
# 라쏘의 경우
lasso.mod=glmnet(x[train ,],y[train],alpha =1, lambda =grid)
plot(lasso.mod)
set.seed(1)
cv.out=cv.glmnet(x[train ,],y[train],alpha =1)
plot(cv.out)
bestlam=cv.out$lambda.min
bestlam
## [1] 16.78016
lasso.pred=predict(lasso.mod,s=bestlam,newx=x[test ,])
mean((lasso.pred-y.test)^2)
## [1] 100743.4
out=glmnet(x,y,alpha=1,lambda=grid)
lasso.coef=predict(out,type="coefficients",s=bestlam )[1:20 ,]
lasso.coef
## (Intercept) AtBat Hits HmRun Runs
## 18.5394844 0.0000000 1.8735390 0.0000000 0.0000000
## RBI Walks Years CAtBat CHits
## 0.0000000 2.2178444 0.0000000 0.0000000 0.0000000
## CHmRun CRuns CRBI CWalks LeagueN
## 0.0000000 0.2071252 0.4130132 0.0000000 3.2666677
## DivisionW PutOuts Assists Errors NewLeagueN
## -103.4845458 0.2204284 0.0000000 0.0000000 0.0000000
# 다항회귀모형의 적합
fit = lm(wage ~ poly(age, 4), data = Wage)
coef(summary(fit))
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 111.70361 0.7287409 153.283015 0.000000e+00
## poly(age, 4)1 447.06785 39.9147851 11.200558 1.484604e-28
## poly(age, 4)2 -478.31581 39.9147851 -11.983424 2.355831e-32
## poly(age, 4)3 125.52169 39.9147851 3.144742 1.678622e-03
## poly(age, 4)4 -77.91118 39.9147851 -1.951938 5.103865e-02
fit2 = lm(wage ~ poly(age, 4, raw = TRUE), data = Wage)
coef(summary(fit2))
fit2a = lm(wage ~ age + I(age^2) + I(age^3) + I(age^4), data = Wage)
coef(summary(fit2a))
fit2b = lm(wage ~ cbind(age, age^2, age^3, age^4), data = Wage)
coef(summary(fit2b))
agelims = range(Wage$age)
age.grid = seq(from = agelims[1], to = agelims[2])
preds = predict(fit, newdata = list(age = age.grid), se = TRUE)
se.bands = cbind(preds$fit + 2*preds$se.fit, preds$fit - 2*preds$se.fit)
attach(Wage)
plot(age, wage, xlim = agelims, cex = .5, col = "darkgrey", main = "Degree-4 Polynomial")
lines(age.grid, preds$fit, lwd = 2, col = "blue")
matlines(age.grid, se.bands, lwd=1, col = "blue", lty = 3)
# 회귀스플라인 (Regression Spline)
library(splines)
fit = lm(wage~bs(age, knots=c(25,40,60)), data = Wage)
pred = predict(fit, newdata = list(age = age.grid), se = TRUE)
plot(age, wage, col = "gray")
lines(age.grid, pred$fit, lwd = 2)
lines(age.grid, pred$fit + 2*pred$se, lty = "dashed")
lines(age.grid, pred$fit - 2*pred$se, lty = "dashed")
# 자연스플라인 (Natural Spline)
fit2 = lm(wage ~ ns(age, df=4), data = Wage)
pred2 = predict(fit2, newdata = list(age = age.grid), se = TRUE)
plot(age, wage, col = "gray")
lines(age.grid, pred2$fit, col = "red", lwd = 2)
# 평활스플라인 (Smoothing Spline)
plot(age, wage, xlim=agelims, cex = .5, col = "darkgrey", main = "Smoothing Spline")
fit = smooth.spline(age, wage, df = 16)
fit2 = smooth.spline(age, wage, cv = TRUE)
fit2$df
## [1] 6.794596
lines(fit, col = "red", lwd = 2)
lines(fit2, col = "blue", lwd = 2)
legend("topright",legend=c("16 DF","6.8 DF"),col=c("red "," blue "),lty =1, lwd =2, cex =.8)
# 국소회귀 (Local Regression)
plot(age,wage,xlim=agelims,cex=.5,col="darkgrey")
title("Local Regression")
fit=loess(wage~age,span=.2,data=Wage)
fit2=loess(wage~age,span=.5,data=Wage)
lines(age.grid,predict(fit,data.frame(age=age.grid)),col="red",lwd=2)
lines(age.grid,predict(fit2,data.frame(age=age.grid)),col="blue",lwd=2)
legend("topright",legend=c("Span=0.2","Span=0.5"),col=c("red","blue"),lty=1,lwd=2,cex=.8)
# 일반화가법모형 (Generalized Additive Models)
gam1 = lm(wage~ns(year, 4)+ns(age, 5) + education, data=Wage)
library(gam)
gam.m3 = gam(wage~s(year, 4) + s(age, 5) + education, data=Wage)
summary(gam.m3)
##
## Call: gam(formula = wage ~ s(year, 4) + s(age, 5) + education, data = Wage)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -119.43 -19.70 -3.33 14.17 213.48
##
## (Dispersion Parameter for gaussian family taken to be 1235.69)
##
## Null Deviance: 5222086 on 2999 degrees of freedom
## Residual Deviance: 3689770 on 2986 degrees of freedom
## AIC: 29887.75
##
## Number of Local Scoring Iterations: 2
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## s(year, 4) 1 27162 27162 21.981 2.877e-06 ***
## s(age, 5) 1 195338 195338 158.081 < 2.2e-16 ***
## education 4 1069726 267432 216.423 < 2.2e-16 ***
## Residuals 2986 3689770 1236
## ---
## Signif. codes: 0
***
0.001
**
0.01
*
0.05
.
0.1
1
##
## Anova for Nonparametric Effects
## Npar Df Npar F Pr(F)
## (Intercept)
## s(year, 4) 3 1.086 0.3537
## s(age, 5) 4 32.380 <2e-16 ***
## education
## ---
preds = predict(gam.m3, newdata = Wage)
par(mfrow = c(1,3))
plot(gam.m3, se=TRUE, col = "blue")
plot.gam(gam1, se=TRUE, col = "red")
gam.m1=gam(wage~s(age,5)+education,data=Wage)
gam.m2=gam(wage~year+s(age,5)+education,data=Wage)
anova(gam.m1,gam.m2,gam.m3,test="F")
## Analysis of Deviance Table
##
## Model 1: wage ~ s(age, 5) + education
## Model 2: wage ~ year + s(age, 5) + education
## Model 3: wage ~ s(year, 4) + s(age, 5) + education
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 2990 3711731
## 2 2989 3693842 1 17889.2 14.4771 0.0001447 ***
## 3 2986 3689770 3 4071.1 1.0982 0.3485661
## ---
gam.lo=gam(wage~s(year,df=4)+lo(age,span=0.7)+education,data=Wage)
par(mfrow=c(1,3))
plot.gam(gam.lo, se=TRUE, col="green")
gam.lo.i=gam(wage~lo(year,age,span=0.5)+education,data=Wage)
library(akima)
par(mfrow=c(1,2))
plot(gam.lo.i)
table(education,I(wage>250))
##
## education FALSE TRUE
## 1. < HS Grad 268 0
## 2. HS Grad 966 5
## 3. Some College 643 7
## 4. College Grad 663 22
## 5. Advanced Degree 381 45
par(mfrow=c(1,3))
gam.lr.s=gam(I(wage>250)~year+s(age, df=5)+education, family=binomial
, data=Wage, subset =( education !="1. < HS Grad"))
plot(gam.lr.s, se=T, col="green")
4일 : 분류, 군집 (Classification, Clustering)
[실습]
-
분류
## Data for An Introduction to Statistical Learning with Applications in R
library(ISLR)
## for lda() and qda()
library(MASS)
## for knn()
library(class)
## Plotting
library(ggplot2)
library(gridExtra)
data(package = "ISLR")
## Load Smarket (S&P Stock Market Data)
data(Smarket)
head(Smarket)
names(Smarket)
dim(Smarket)
summary(Smarket)
pairs(Smarket)
attach(Smarket)
# Logistic Regression
glm.fit=glm(Direction~Lag1+Lag2+Lag3+Lag4+Lag5+Volume,data=Smarket,family=binomial)
summary(glm.fit)
coef(glm.fit)
summary(glm.fit)$coef
summary(glm.fit)$coef[,4]
glm.probs=predict(glm.fit,type="response")
glm.probs[1:10]
contrasts(Direction)
glm.pred=rep("Down",1250)
glm.pred[glm.probs>.5]="Up"
table(glm.pred,Direction)
(507+145)/1250
mean(glm.pred==Direction)
train=(Year<2005)
Smarket.2005=Smarket[!train,]
dim(Smarket.2005)
Direction.2005=Direction[!train]
glm.fit=glm(Direction~Lag1+Lag2+Lag3+Lag4+Lag5+Volume,data=Smarket,family=binomial,subset=train)
glm.probs=predict(glm.fit,Smarket.2005,type="response")
glm.pred=rep("Down",252)
glm.pred[glm.probs>.5]="Up"
table(glm.pred,Direction.2005)
mean(glm.pred==Direction.2005)
mean(glm.pred!=Direction.2005)
## Plot
plotData <- ggplot(data = Smarket,
mapping = aes(x = Lag1, y = Lag2, color = Direction, shape = factor(Year))) +
geom_point(stat = "identity", alpha = 0.8) +
scale_color_manual(values = c("Down" = "yellow", "Up" = "red")) +
theme_bw() +
theme(legend.key = element_blank()) +
labs(title = "Original data")
plotData
glm.fit=glm(Direction~Lag1+Lag2,data=Smarket,family=binomial(link = "logit"),subset=train)
glm.probs=predict(glm.fit,Smarket.2005,type="response")
glm.pred=rep("Down",252)
glm.pred[glm.probs>.5]="Up"
table(glm.pred,Direction.2005)
mean(glm.pred==Direction.2005)
106/(106+76)
predict(glm.fit,newdata=data.frame(Lag1=c(1.2,1.5),Lag2=c(1.1,-0.8)),type="response")
## Plot
plotData <- ggplot(data = Smarket.2005,
mapping = aes(x = Lag1, y = Lag2, color = Direction, shape = factor(Year))) +
geom_point(stat = "identity", alpha = 0.8) +
scale_color_manual(values = c("Down" = "yellow", "Up" = "red")) +
theme_bw() +
theme(legend.key = element_blank()) +
labs(title = "Test data")
plotData
## Put the predicted probability
predProbLogit <- predict(glm.fit,Smarket.2005,type="response")
predClassLogit <- factor(predict(glm.fit,Smarket.2005,type="response") > 0.5,
levels = c(FALSE,TRUE), labels = c("Down","Up"))
## Plot (probability)
plotLogit <- ggplot(data = Smarket.2005,
mapping = aes(x = Lag1, y = Lag2, color = predProbLogit, shape = factor(Year))) +
geom_point(alpha = 0.8) +
scale_color_gradient(low = "yellow", high = "red") +
theme_bw() +
theme(legend.key = element_blank()) +
labs(title = "Predicted probability of outcome (Logistic)")
grid.arrange(plotData, plotLogit, ncol = 2)
plotLdaClass <- ggplot(data = Smarket.2005,
mapping = aes(x = Lag1, y = Lag2, color = predClassLogit, shape = factor(Year))) +
geom_point(alpha = 0.7) +
scale_color_manual(values = c("Down" = "yellow", "Up" = "red")) +
theme_bw() +
theme(legend.key = element_blank()) +
labs(title = "Predicted outcome (Logistic; p>0.5)")
grid.arrange(plotData, plotLdaClass, ncol = 2)
#Linear Discriminant Analysis
library(MASS)
lda.fit=lda(Direction~Lag1+Lag2,data=Smarket,subset=train)
lda.fit
plot(lda.fit)
lda.pred=predict(lda.fit, Smarket.2005)
names(lda.pred)
lda.class=lda.pred$class
table(lda.class,Direction.2005)
mean(lda.class==Direction.2005)
sum(lda.pred$posterior[,1]>=.5)
sum(lda.pred$posterior[,1]<.5)
lda.pred$posterior[1:20,1]
lda.class[1:20]
sum(lda.pred$posterior[,1]>.9)
## Put into the dataset
predProbLda <- lda.pred$posterior[,"Up"]
predClassLda <- lda.pred$class
## Plot (probability)
plotLdaProb <- ggplot(Smarket.2005,
mapping = aes(x = Lag1, y = Lag2, color = predProbLda, shape = factor(Year))) +
geom_point(alpha = 0.5) +
scale_color_gradient(low = "yellow", high = "red") +
theme_bw() +
theme(legend.key = element_blank()) +
labs(title = "Predicted probability of outcome (LDA)")
grid.arrange(plotData, plotLdaProb, ncol = 2)
## Plot (classification)
plotLdaClass <- ggplot(Smarket.2005,
mapping = aes(x = Lag1, y = Lag2, color = predClassLda, shape = factor(Year))) +
geom_point(alpha = 0.5) +
scale_color_manual(values = c("Down" = "yellow", "Up" = "red")) +
theme_bw() +
theme(legend.key = element_blank()) +
labs(title = "Predicted outcome (LDA)")
grid.arrange(plotData, plotLdaClass, ncol = 2)
# Quadratic Discriminant Analysis
qda.fit=qda(Direction~Lag1+Lag2,data=Smarket,subset=train)
qda.fit
qda.pred=predict(qda.fit,Smarket.2005)
qda.class=qda.pred$class
table(qda.class,Direction.2005)
mean(qda.class==Direction.2005)
## Put into the dataset
predProbQda <- qda.pred$posterior[,"Up"]
predClassQda <- qda.pred$class
## Plot (probability)
plotQdaProb <- ggplot(Smarket.2005,
mapping = aes(x = Lag1, y = Lag2, color = predProbQda, shape = factor(Year))) +
geom_point(alpha = 0.5) +
scale_color_gradient(low = "yellow", high = "red") +
theme_bw() +
theme(legend.key = element_blank()) +
labs(title = "Predicted probability of outcome (QDA)")
grid.arrange(plotData, plotQdaProb, ncol = 2)
## Plot (classification)
plotQdaClass <- ggplot(Smarket.2005,
mapping = aes(x = Lag1, y = Lag2, color = predClassQda, shape = factor(Year))) +
geom_point(alpha = 0.5) +
scale_color_manual(values = c("Down" = "yellow", "Up" = "red")) +
theme_bw() +
theme(legend.key = element_blank()) +
labs(title = "Predicted outcome (QDA)")
grid.arrange(plotData, plotQdaClass, ncol = 2)
# K-Nearest Neighbors
library(class)
train.X=cbind(Lag1,Lag2)[train,]
test.X=cbind(Lag1,Lag2)[!train,]
train.Direction=Direction[train]
set.seed(1)
knn.pred=knn(train.X,test.X,train.Direction,k=1)
table(knn.pred,Direction.2005)
(83+43)/252
knn.pred=knn(train.X,test.X,train.Direction,k=3)
table(knn.pred,Direction.2005)
mean(knn.pred==Direction.2005)
## Fit KNN (the output is a vector of prediction)
knn.pred <- knn(train = train.X,
test = test.X,
cl = train.Direction, k = 3)
predClassKnn <- knn.pred
## Plot (classification)
plotKnnClass <- ggplot(Smarket.2005,
mapping = aes(x = Lag1, y = Lag2, color = predClassKnn, shape = factor(Year))) +
geom_point(alpha = 0.8) +
scale_color_manual(values = c("Down" = "yellow", "Up" = "red")) +
theme_bw() +
theme(legend.key = element_blank()) +
labs(title = "Predicted outcome (KNN)")
grid.arrange(plotData, plotKnnClass, ncol = 2)
-
군집분석, 요인분석
# K-Means Clustering
library(ggplot2)
set.seed(2)
nd=matrix(rnorm(50*2), ncol=2)
nd[1:25,1]=nd[1:25,1]+3
nd[1:25,2]=nd[1:25,2]-4
colnames(nd)=c("x1","x2")
km.out=kmeans(nd,2,nstart=20)
km.out$cluster
cluster=factor(km.out$cluster)
nd1=data.frame(nd,cluster)
c1=ggplot(nd1, aes(x1, x2 , color = cluster)) +
geom_point(alpha = 0.4, size = 3.5) + labs(title = "K-Means Clustering Results with K=2")
print(c1)
dev.off()
set.seed(4)
km.out=kmeans(nd,3,nstart=20)
km.out
cluster=factor(km.out$cluster)
nd2=data.frame(nd,cluster)
c2=ggplot(nd2, aes(x1, x2 , color = cluster)) +
geom_point(size = 3.5) + labs(title = "K-Means Clustering Results with K=3")
print(c2)
centers=as.data.frame(km.out$centers)
ggplot(data=nd2, aes(x=x1, y=x2,color=cluster )) +
geom_point() +
geom_point(data=centers, aes(x=x1,y=x2, color='Center')) +
geom_point(data=centers, aes(x=x1,y=x2, color='Center'), size=52, alpha=.3, show.legend=FALSE)
set.seed(3)
km.out=kmeans(nd,3,nstart=1)
km.out$tot.withinss
km.out=kmeans(nd,3,nstart=20)
km.out$tot.withinss
#cf)
library(ggfortify)
autoplot(km.out, data=nd1,frame = TRUE)
autoplot(km.out, data=nd1,frame = TRUE, frame.type = 'norm')
autoplot(km.out, data=nd2, frame.type = "convex", frame.level = 0.95,
frame.alpha = 0.2, main="K-Means Clustering Results with K=3")
# Hierarchical Clustering
hc.complete=hclust(dist(nd), method="complete")
hc.average=hclust(dist(nd), method="average")
hc.single=hclust(dist(nd), method="single")
#your own labels (now rownames) are supplied in geom_text() and label=label
dendr= dendro_data(hc.complete, type="rectangle")
ggplot() +
geom_segment(data=segment(dendr), aes(x=x, y=y, xend=xend, yend=yend)) +
geom_text(data=label(dendr), aes(x=x, y=y, label=label, hjust=0), size=4) +
coord_flip() + scale_y_reverse(expand=c(0.2, 0)) +
theme(axis.line.y=element_blank(),
axis.ticks.y=element_blank(),
axis.text.y=element_blank(),
axis.title.y=element_blank(),
panel.background=element_rect(fill="white"),
panel.grid=element_blank())+
labs(title="Complete Linkage")
dendr = dendro_data(hc.average, type="rectangle")
ggplot() +
geom_segment(data=segment(dendr), aes(x=x, y=y, xend=xend, yend=yend)) +
geom_text(data=label(dendr), aes(x=x, y=y, label=label, hjust=0), size=4) +
coord_flip() + scale_y_reverse(expand=c(0.2, 0)) +
theme(axis.line.y=element_blank(),
axis.ticks.y=element_blank(),
axis.text.y=element_blank(),
axis.title.y=element_blank(),
panel.background=element_rect(fill="white"),
panel.grid=element_blank())+
labs(title="Average Linkage")
dendr = dendro_data(hc.single, type="rectangle")
ggplot() +
geom_segment(data=segment(dendr), aes(x=x, y=y, xend=xend, yend=yend)) +
geom_text(data=label(dendr), aes(x=x, y=y, label=label, hjust=0), size=4) +
coord_flip() + scale_y_reverse(expand=c(0.2, 0)) +
theme(axis.line.y=element_blank(),
axis.ticks.y=element_blank(),
axis.text.y=element_blank(),
axis.title.y=element_blank(),
panel.background=element_rect(fill="white"),
panel.grid=element_blank())+
labs(title="Single Linkage")
ndsc=scale(nd)
hcsc.complete=hclust(dist(nd), method="complete")
dendr = dendro_data(hcsc.complete, type="rectangle")
ggplot() +
geom_segment(data=segment(dendr), aes(x=x, y=y, xend=xend, yend=yend)) +
geom_text(data=label(dendr), aes(x=x, y=y, label=label, hjust=0), size=4) +
coord_flip() + scale_y_reverse(expand=c(0.2, 0)) +
theme(axis.line.y=element_blank(),
axis.ticks.y=element_blank(),
axis.text.y=element_blank(),
axis.title.y=element_blank(),
panel.background=element_rect(fill="white"),
panel.grid=element_blank())+
labs(title="Hierarchical Clustering with Scaled Features")
nd3=matrix(rnorm(30*3), ncol=3)
ndd=as.dist(1-cor(t(nd3)))
hcc.complete=hclust(dist(ndd), method="complete")
dendr = dendro_data(hcc.complete, type="rectangle")
ggplot() +
geom_segment(data=segment(dendr), aes(x=x, y=y, xend=xend, yend=yend)) +
geom_text(data=label(dendr), aes(x=x, y=y, label=label, hjust=0), size=4) +
coord_flip() + scale_y_reverse(expand=c(0.2, 0)) +
theme(axis.line.y=element_blank(),
axis.ticks.y=element_blank(),
axis.text.y=element_blank(),
axis.title.y=element_blank(),
panel.background=element_rect(fill="white"),
panel.grid=element_blank())+
labs(title="Complete Linkage with Correlation-Based Distance")
cutree(hc.complete, 2)
cutree(hc.average, 2)
cutree(hc.single, 2)
cutree(hc.single, 4)
hc.average=hclust(dist(nd), method="average") # heirarchal clustering
dendr = dendro_data(hc.average, type="rectangle") # convert for ggplot
clust = cutree(hc.average,k=3) # find 2 clusters
names(clust)=c(1:50)
clust.df <- data.frame(label=names(clust), cluster=factor(clust))
# dendr[["labels"]] has the labels, merge with clust.df based on label column
dendr[["labels"]] <- merge(dendr[["labels"]],clust.df, by="label")
# plot the dendrogram; note use of color=cluster in geom_text(...)
ggplot() +
geom_segment(data=segment(dendr), aes(x=x, y=y, xend=xend, yend=yend)) +
geom_text(data=label(dendr), aes(x, y, label=label, hjust=0, color=cluster),
size=3) +
coord_flip() + scale_y_reverse(expand=c(0.2, 0)) +
theme(axis.line.y=element_blank(),
axis.ticks.y=element_blank(),
axis.text.y=element_blank(),
axis.title.y=element_blank(),
panel.background=element_rect(fill="white"),
panel.grid=element_blank())
#(cf)
# demonstrate plotting directly from object class hclust
ggdendrogram(hc.complete)
ggdendrogram(hc.complete, rotate = TRUE, size = 5, theme_dendro = FALSE)
# demonstrate converting hclust to dendro using dendro_data first
hcdata <- dendro_data(hc.complete)
ggdendrogram(hcdata, rotate=TRUE) +
labs(title="Dendrogram in ggplot2")
# Triangular lines
ddata <- dendro_data(as.dendrogram(hc.complete), type = "triangle")
ggplot(segment(ddata)) + geom_segment(aes(x = x, y = y, xend = xend,
yend = yend)) +
geom_text(data = label(ddata), aes(x = x,y = y, label = label), angle = 90, lineheight = 0)
# factor analysis
life<-read.table("life_expectation.txt",header=TRUE)
cor.life<-cor(life)
library(corrplot)
corrplot(cor.life, type="upper", order="hclust", tl.col="black", tl.srt=45)
d.factanal=factanal(life,factors=1, ,method="mle",rotation="varimax")
d.factanal=factanal(life,factors=2, method="mle",rotation="varimax")
d.factanal=factanal(life,factors=3, method="mle",rotation="varimax")
# How could we interpret the three factors?
sapply(1:3, function(f)
factanal(life, factors = f, method="mle")$PVAL)
# The communalities:
life.exp <- factanal(life,factors=3, rotation="varimax")
1-life.exp$uniquenesses
### Plot loadings against one another
# Factor 1 vs. Factor 2
text=rownames(d.factanal$loadings)
load = d.factanal$loadings
ggplot(data.frame(load[,1:2]),aes(Factor1, Factor2,label = text))+
geom_point() +
geom_text(angle = 20,size=5)+
labs(title="Factor 1 vs Factor 2", x="Factor1 loadings", y="Factor 2 loadings")
# Factor 1 vs. Factor 3
text=rownames(d.factanal$loadings)
load = d.factanal$loadings
ggplot(data.frame(load[,c(1,3)]),aes(Factor1, Factor3,label = text))+
geom_point() +
geom_text(angle = 20,size=5)+
labs(title="Factor 1 vs Factor 2", x="Factor1 loadings", y="Factor 3 loadings")
# Factor 2 vs. Factor 3:
ggplot(data.frame(load[, 2:3]),aes(Factor2, Factor3,label = text))+
geom_point() +
geom_text(angle = 20,size=5)+
labs(title="Factor 2 vs Factor 3", x="Factor2 loadings", y="Factor 3 loadings")
# Estimating factor scores for the life expectancy data set:
scores <- factanal(life,factors=3, rotation="varimax",scores="regression")$scores
### Doing separate 2-D scatterplots of the factor scores:
# Factor 1 vs. Factor 2
scores[,c(1,2,3)]
text=rownames(scores)
ggplot(data.frame(scores[, 1:2]),aes(Factor1, Factor2,label = text))+
geom_point() +
geom_text(angle = 20)+
labs(title="Factor 1 vs Factor 2", x="Factor1 scores", y="Factor 2 scores")
# Factor 1 vs. Factor 3:
ggplot(data.frame(scores[, c(1,3)]),aes(Factor1, Factor3,label = text))+
geom_point() +
geom_text(angle = 20)+
labs(title="Factor 1 vs Factor 3", x="Factor1 scores", y="Factor 3 scores")
# Factor 2 vs. Factor 3:
ggplot(data.frame(scores[, 2:3]),aes(Factor2, Factor3,label = text))+
geom_point() +
geom_text(angle = 20)+
labs(title="Factor 2 vs Factor 3", x="Factor2 scores", y="Factor 3 scores")
#(cf)1.
d.factanal = factanal(life, factors = 3, scores = 'regression')
autoplot(d.factanal, label = TRUE, label.size = 3,
loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3)
#(cf)2.
#Using fa() from package psych with rotation
install.packages("psych")
library(psych)
d.fa = fa(r=cor.life, nfactors=3, n.obs=31, rotate="varimax")
d.fa
#Factor scores
scores <- factor.scores(x=life, f=d.fa, method="regression")
head(scores)
#Visualize loadings
factor.plot(d.fa, cut=0.5)
-
나무모형
### Decision Trees
# Fitting Classification Trees
Carseats<-read.table("Carseats.txt")
library(tree)
library(ggplot2)
library(ggdendro)
attach(Carseats)
High=ifelse(Sales<=8,"No","Yes")
Carseats=data.frame(Carseats,High)
tree.carseats=tree(High~.-Sales, Carseats)
plot(tree.carseats)
text(tree.carseats)
summary(tree.carseats)
dev.off()
# Confusion matrix
set.seed(2)
train=sample(1:nrow(Carseats), 200)
Carseats.test=Carseats[-train,]
High.test=High[-train]
tree.carseats=tree(High~.-Sales,Carseats,subset=train)
tree.pred=predict(tree.carseats,Carseats.test,type="class")
table(tree.pred,High.test)
(86+57)/200
# perform cross-validation in order to determine the optimal level of tree complexity
set.seed(3)
cv.carseats=cv.tree(tree.carseats,FUN=prune.misclass)
# prune.tree() enacts optimal deviance pruning
# prune.misclass() enacts optimal misclassification rate pruning
names(cv.carseats)
cv.carseats
# plot the error rate as a function of both size and k
library(gridExtra)
cv=data.frame(cv.carseats$size,cv.carseats$k,cv.carseats$dev)
colnames(cv) = c("size","k","dev")
min=subset(cv, dev==min(dev))
cv1=ggplot(cv,aes(x=size,y=dev))+
geom_point()+
geom_point(data=min, size=3, colour=alpha("red", 0.7))+
geom_line()
cv2=ggplot(cv,aes(x=k,y=dev))+
geom_point()+
geom_point(data=min, size=3, colour=alpha("red", 0.7))+
geom_line()
grid.arrange(cv1, cv2, nrow=1, ncol=2)
prune.carseats=prune.misclass(tree.carseats,best=9)
plot(prune.carseats)
text(prune.carseats,pretty=0)
# Confusion matrix
tree.pred=predict(prune.carseats,Carseats.test,type="class")
table(tree.pred,High.test)
(94+60)/200
prune.carseats=prune.misclass(tree.carseats,best=15)
plot(prune.carseats)
text(prune.carseats,pretty=0)
tree.pred=predict(prune.carseats,Carseats.test,type="class")
table(tree.pred,High.test)
(86+62)/200
# (cf)
library(rpart)
tree.carseats = rpart(High~.-Sales, method = "class", data = Carseats)
# display the results
printcp(tree.carseats)
# Classification tree diagrams
ddata <- dendro_data(tree.carseats)
ggplot() +
geom_segment(data = ddata$segments,
aes(x = x, y = y, xend = xend, yend = yend)) +
geom_text(data = ddata$labels,
aes(x = x, y = y, label = label), size = 4, vjust = 0) +
geom_text(data = ddata$leaf_labels,
aes(x = x, y = y, label = label), size = 4, vjust = 1) +
theme_dendro()
# visualize cross-validation results
plotcp(tree.carseats)
# plot tree
plot(tree.carseats, uniform=TRUE,
main="Classification Tree for Carseats")
text(tree.carseats, use.n=TRUE, all=TRUE, cex=.8)
# prune the tree
pfit<- prune(tree.carseats, cp=tree.carseats$cptable[which.min(tree.carseats$cptable[,"xerror"]),"CP"])
# plot the pruned tree
plot(pfit, uniform=TRUE,
main="Pruned Classification Tree for Carseats")
text(pfit, use.n=F, all=TRUE, cex=.8)
# perform cross-validation in order to determine the optimal level of tree complexity
ddata <- dendro_data(pfit)
ggplot() +
geom_segment(data = ddata$segments,
aes(x = x, y = y, xend = xend, yend = yend)) +
geom_text(data = ddata$labels,
aes(x = x, y = y, label = label), size = 4, vjust = 0) +
geom_text(data = ddata$leaf_labels,
aes(x = x, y = y, label = label), size = 4, vjust = 1) +
theme_dendro()
install.packages("rpart.plot")
library(rpart.plot)
prp(pfit, faclen = 0, cex = 0.8, extra = 1)
only_count <- function(x, labs, digits, varlen)
{
paste(x$frame$n)
}
boxcols <- c("pink", "lightblue")[pfit$frame$yval]
par(xpd=TRUE)
prp(pfit, faclen = 0, cex = 0.8, node.fun=only_count, box.col = boxcols)
legend("bottomright", legend = c("No","Yes"), fill = c("pink", "lightblue"),
title = "Group")
# Fitting Regression Trees
library(MASS)
set.seed(1)
train = sample(1:nrow(Boston), nrow(Boston)/2)
tree.boston=tree(medv~.,Boston,subset=train)
summary(tree.boston)
# Regression tree diagrams
tree_data <- dendro_data(tree.boston)
ggplot(segment(tree_data)) +
geom_segment(aes(x = x, y = y, xend = xend, yend = yend, size = n),
colour = "blue", alpha = 0.5) +
scale_size("n") +
geom_text(data = label(tree_data),
aes(x = x, y = y, label = label), vjust = -0.5, size = 4) +
geom_text(data = leaf_label(tree_data),
aes(x = x, y = y, label = label), vjust = 0.5, size = 3) +
theme_dendro()+
labs(title="Regression tree diagrams")
# perform cross-validation in order to determine the optimal level of tree complexity
cv.boston=cv.tree(tree.boston)
# plot the error rate as a function of size
cv=data.frame(cv.boston$size,cv.boston$k,cv.boston$dev)
colnames(cv) = c("size","k","dev")
min=subset(cv, dev==min(dev))
ggplot(cv,aes(x=size,y=dev))+
geom_point()+
geom_point(data=min, size=3, colour=alpha("red", 0.7))+
geom_line()
# plot the pruned tree
prune.boston=prune.tree(tree.boston,best=5)
tree_data <- dendro_data(prune.boston)
ggplot(segment(tree_data)) +
geom_segment(aes(x = x, y = y, xend = xend, yend = yend, size = n),
colour = "blue", alpha = 0.5) +
scale_size("n") +
geom_text(data = label(tree_data),
aes(x = x, y = y, label = label), vjust = -0.5, size = 4) +
geom_text(data = leaf_label(tree_data),
aes(x = x, y = y, label = label), vjust = 0.5, size = 3) +
theme_dendro()+
labs(title="Regression tree diagrams")
# the test set MSE associated with the regression tree
yhat=predict(tree.boston,newdata=Boston[-train,])
boston.test=Boston[-train,"medv"]
mean((yhat-boston.test)^2)
ggplot()+
geom_point(aes(x=yhat, y=boston.test))+
geom_abline(intercept = 0, slope = 1, color="red",
linetype="dashed", size=1.5)+
xlab(expression(hat(y)))+
ylab(expression(y[0]))
# Bagging and Random Forests
# Bagging
library(randomForest)
set.seed(1)
bag.boston=randomForest(medv~.,data=Boston,subset=train,mtry=13,importance=TRUE)
bag.boston
# the test MSE for bagging
yhat.bag = predict(bag.boston,newdata=Boston[-train,])
mean((yhat.bag-boston.test)^2)
ggplot()+
geom_point(aes(x=yhat.bag, y=boston.test))+
geom_abline(intercept = 0, slope = 1, color="red",
linetype="dashed", size=1.5)+
xlab(expression(hat(y)))+
ylab(expression(y[0]))
bag.boston=randomForest(medv~.,data=Boston,subset=train,mtry=13,ntree=25)
yhat.bag = predict(bag.boston,newdata=Boston[-train,])
mean((yhat.bag-boston.test)^2)
# Random Forests
set.seed(1)
rf.boston=randomForest(medv~.,data=Boston,subset=train,mtry=6,importance=TRUE)
# the test MSE for random forests
yhat.rf = predict(rf.boston,newdata=Boston[-train,])
mean((yhat.rf-boston.test)^2)
# Variable Importance Plot
var_importance <- data.frame(variable=setdiff(colnames(Boston), "medv"),
importance=as.vector(rf.boston$importance[,1]))
var_importance <- var_importance[order(var_importance[,2]),]
var_importance$variable <- factor(var_importance$variable, levels=var_importance$variable)
vi=ggplot(var_importance,aes(x=variable, y=importance, fill=variable)) +
geom_bar(stat='identity', position='identity')+
scale_fill_discrete(name="Variable Name",guide = guide_legend(reverse=TRUE)) +
coord_flip()+
theme(axis.text.x=element_blank(),
axis.text.y=element_text(size=12),
axis.title=element_text(size=16),
plot.title=element_text(size=18),
legend.title=element_text(size=16),
legend.text=element_text(size=12))+
ylab("Mean square error")
var_purity <- data.frame(variable=setdiff(colnames(Boston), "medv"),
purity=as.vector(rf.boston$importance[,2]))
var_purity = var_purity[order(var_purity[,2]),]
var_purity$variable <- factor(var_purity$variable, levels=var_purity$variable)
np=ggplot(var_purity,aes(x=variable, y=purity, fill=variable)) +
geom_bar(stat='identity', position='identity') +
scale_fill_discrete(name="Variable Name",guide = guide_legend(reverse=TRUE)) +
coord_flip()+
theme(axis.text.x=element_blank(),
axis.text.y=element_text(size=12),
axis.title=element_text(size=16),
plot.title=element_text(size=18),
legend.title=element_text(size=16),
legend.text=element_text(size=12))+
ylab("Node impurity")
grid.arrange(vi, np, nrow=1, ncol=2,top="Variable Importance from Random Forest Fit")
varImpPlot(rf.boston)
# Boosting
library(gbm)
set.seed(1)
boost.boston=gbm(medv~.,data=Boston[train,],distribution="gaussian",n.trees=5000,interaction.depth=4)
summary(boost.boston)
# partial dependence plots for rm and lstat.
par(mfrow=c(1,2))
plot(boost.boston,i="rm",col="red")
plot(boost.boston,i="lstat",col="blue")
# the test MSE for boosting
yhat.boost=predict(boost.boston,newdata=Boston[-train,],n.trees=5000)
mean((yhat.boost-boston.test)^2)
boost.boston=gbm(medv~.,data=Boston[train,],distribution="gaussian",n.trees=5000,interaction.depth=4,shrinkage=0.2,verbose=F)
yhat.boost=predict(boost.boston,newdata=Boston[-train,],n.trees=5000)
mean((yhat.boost-boston.test)^2)
[연습문제]
-
분류
# Logistic Regression
# Exercise1
#(1)
titanic = read.csv("titanic3.csv")
dim(titanic)
str(titanic)
summary(titanic)
sum(is.na(titanic))
titanic = titanic[,c(1,2,4:7)]
titanic=na.omit(titanic)
glm.fit= glm(survived ~ pclass + sex + age + sibsp,
family = binomial, data = titanic)
summary(glm.fit)
coef(glm.fit)
summary(glm.fit)$coef
summary(glm.fit)$coef[,4]
#(2)
attach(titanic)
train = (parch < 1)
test = !train
train.titanic = titanic[train, ]
test.titanic = titanic[test, ]
test.survived = survived[test]
#(3)
glm.fit=glm(survived ~ pclass + sex + age,
family = binomial,subset=train)
glm.probs=predict(glm.fit,test.titanic,type="response")
glm.pred=rep(0,length(glm.probs))
glm.pred[glm.probs>.5]=1
mean(glm.pred!=test.survived)
# Linear Discriminant Analysis
# Exercise 2
#(1)
library(MASS)
lda.fit=lda(survived ~ pclass + sex + age + sibsp,data=titanic, subset=train)
lda.fit
#(2)
lda.pred=predict(lda.fit,test.titanic)
names(lda.pred)
lda.class=lda.pred$class
table(lda.class,test.survived)
mean(lda.class!==test.survived)
# Quadratic Discriminant Analysis
# Exercise 3
#(1)
qda.fit=qda(survived ~ pclass + sex + age + sibsp,data=titanic, subset=train)
qda.fit
#(2)
qda.class=predict(qda.fit ,test.titanic)$class
table(qda.class,test.survived)
mean(qda.class!=test.survived)
# K-Nearest Neighbors
# Exercise 4
#(1)
library(class)
train.x = cbind(pclass, sex, age, sibsp)[train, ]
rownames(train.x)=NULL
test.x = cbind(pclass, sex, age, sibsp)[test, ]
rownames(test.x)=NULL
train.survived = survived[train]
# K=1
set.seed(0413)
knn.pred=knn(train.x,test.x,train.survived,k=1)
mean(knn.pred != test.survived)
# K=10
knn.pred=knn(train.x,test.x,train.survived,k=3)
mean(knn.pred != test.survived)
# K=100
knn.pred = knn(train.x,test.x,train.survived,k = 100)
mean(knn.pred != test.survived)
-
군집분석, 요인분석
# K-Means Clustering
# Exercise 10
library(ggfortify)
geneData = read.csv("Ch10Ex11.csv", header=F)
dim(geneData)
geneData = t(geneData)
rownames(geneData)=c(paste("H", 1:20, sep = ""),paste("D", 1:20, sep = ""))
colnames(geneData)=c(paste("g",1:1000,sep=""))
set.seed(0413)
km.out=kmeans(geneData,2,nstart=20)
cluster = km.out$cluster
cluster
autoplot(km.out,data=geneData,label = TRUE, label.size = 3, frame=TRUE)
# Hierarchical Clustering
# Exercise 11
geneData = read.csv("Ch10Ex11.csv", header=F)
dim(geneData)
geneData = t(geneData)
rownames(geneData)=c(paste("H", 1:20, sep = ""),paste("D", 1:20, sep = ""))
colnames(geneData)=c(paste("g",1:1000,sep=""))
dd = as.dist(1 - cor(t(geneData)))
hc.complete=hclust(dd, method="complete")
hc.average=hclust(dd, method="average")
hc.single=hclust(dd, method="single")
library(ggdendro)
library(gridExtra)
dendr = dendro_data(hc.complete, type="rectangle")
h1=ggplot() +
geom_segment(data=segment(dendr), aes(x=x, y=y, xend=xend, yend=yend)) +
geom_text(data=label(dendr), aes(x=x, y=y, label=label, hjust=0), size=4) +
coord_flip() + scale_y_reverse(expand=c(0.2, 0)) +
theme(axis.line.y=element_blank(),
axis.ticks.y=element_blank(),
axis.text.y=element_blank(),
axis.title.y=element_blank(),
panel.background=element_rect(fill="white"),
panel.grid=element_blank())+
labs(title="Complete Linkage with Correlation-Based Distance")
dendr = dendro_data(hc.average, type="rectangle")
h2=ggplot() +
geom_segment(data=segment(dendr), aes(x=x, y=y, xend=xend, yend=yend)) +
geom_text(data=label(dendr), aes(x=x, y=y, label=label, hjust=0), size=4) +
coord_flip() + scale_y_reverse(expand=c(0.2, 0)) +
theme(axis.line.y=element_blank(),
axis.ticks.y=element_blank(),
axis.text.y=element_blank(),
axis.title.y=element_blank(),
panel.background=element_rect(fill="white"),
panel.grid=element_blank())+
labs(title="Average Linkage with Correlation-Based Distance")
dendr = dendro_data(hc.single, type="rectangle")
h3= ggplot() +
geom_segment(data=segment(dendr), aes(x=x, y=y, xend=xend, yend=yend)) +
geom_text(data=label(dendr), aes(x=x, y=y, label=label, hjust=0), size=4) +
coord_flip() + scale_y_reverse(expand=c(0.2, 0)) +
theme(axis.line.y=element_blank(),
axis.ticks.y=element_blank(),
axis.text.y=element_blank(),
axis.title.y=element_blank(),
panel.background=element_rect(fill="white"),
panel.grid=element_blank())+
labs(title="Single Linkage with Correlation-Based Distance")
grid.arrange(h1, h2, h3,nrow=1, ncol=3)
# factor analysis
# Exercise 12
#(1)
pers<-read.csv("personality.csv")
library(corrplot)
corrplot(cor(pers), order = "hclust", tl.col="black", tl.srt=90)
stand.pers = as.data.frame(scale(pers))
#(2)
fac.pers = factanal(stand.pers, factors = 10, rotation = "none", na.action = na.omit)
fac.pers
#3)
1-fac.pers$uniquenesses
# Plot loadings against one another
text=rownames(fac.pers$loadings)
load = fac.pers$loadings
ggplot(data.frame(load[,1:2]),aes(Factor1, Factor2,label = text))+
geom_point() +
geom_text(angle = 20,size=5)+
labs(title="Factor 1 vs Factor 2", x="Factor1 loadings", y="Factor 2 loadings")
#(4)
# Factor analysis with rotation
fac.pers2 = factanal(stand.pers, factors = 10, rotation = "varimax", na.action = na.omit)
fac.pers2
load2 = fac.pers2$loadings[,1:2]
### Plot loadings against one another
text=rownames(fac.pers2$loadings)
load = fac.pers2$loadings
ggplot(data.frame(load[,1:2]),aes(Factor1, Factor2,label = text))+
geom_point() +
geom_text(angle = 20,size=5)+
labs(title="Factor 1 vs Factor 2", x="Factor1 loadings", y="Factor 2 loadings")
-
나무모형
### Decision Trees
# classification Trees
# Exercise 5
#(1)
titanic = read.csv("titanic3.csv")
titanic = titanic[,c(1,2,4:7)]
titanic=na.omit(titanic)
attach(titanic)
train = (parch < 1)
test = !train
train.titanic = titanic[train, ]
test.titanic = titanic[test, ]
test.survived = survived[test]
tree.titanic=tree(as.factor(survived)~ pclass+age+sex+sibsp, titanic,subset=train)
summary(tree.titanic)
#(2)
plot(tree.titanic)
text(tree.titanic,pretty = 0)
tree.pred=predict(tree.titanic,test.titanic, type="class")
table(tree.pred,test.survived)
(70+4)/length(tree.pred)
#(3)
set.seed(0413)
cv.titanic = cv.tree(tree.titanic, FUN = prune.misclass)
names(cv.titanic)
cv.titanic
# plot the error rate as a function of both size and k
library(gridExtra)
cv=data.frame(cv.titanic$size,cv.titanic$k,cv.titanic$dev)
colnames(cv) = c("size","k","dev")
min=subset(cv, dev==min(dev))
cv1=ggplot(cv,aes(x=size,y=dev))+
geom_point()+
geom_point(data=min, size=3, colour=alpha("red", 0.7))+
geom_line()
cv2=ggplot(cv,aes(x=k,y=dev))+
geom_point()+
geom_point(data=min, size=3, colour=alpha("red", 0.7))+
geom_line()
grid.arrange(cv1, cv2, nrow=1, ncol=2)
prune.titanic=prune.misclass(tree.titanic,best=3)
# Confusion matrix
tree.pred=predict(prune.titanic,test.titanic,type="class")
table(tree.pred,test.survived)
(70+4)/length(tree.pred)
# Regression Trees
# Exercise 6
#(1)
library(ISLR)
dim(Hitters)
sum(is.na(Hitters$Salary))
Hitters =na.omit(Hitters)
dim(Hitters)
sum(is.na(Hitters))
train=1:200
train.hitters=Hitters[train,]
test.hitters=Hitters[-train,]
attach(Hitters)
tree.hitters=tree(log(Salary) ~., data=Hitters,subset=train)
tree.hitters
#(2)
tree_data <- dendro_data(tree.hitters)
ggplot(segment(tree_data)) +
geom_segment(aes(x = x, y = y, xend = xend, yend = yend, size = n),
colour = "blue", alpha = 0.5) +
scale_size("n") +
geom_text(data = label(tree_data),
aes(x = x, y = y, label = label), vjust = -0.5, size = 4) +
geom_text(data = leaf_label(tree_data),
aes(x = x, y = y, label = label), vjust = 0.5, size = 3) +
theme_dendro()+
labs(title="Regression tree diagrams")
pred.hitters = predict(tree.hitters, test.hitters)
mean((log(test.hitters$Salary) - pred.hitters)^2)
#(3)
set.seed(0413)
cv.hitters = cv.tree(tree.hitters, FUN = prune.tree)
cv=data.frame(cv.hitters$size,cv.hitters$k,cv.hitters$dev)
colnames(cv) = c("size","k","dev")
min=subset(cv, dev==min(dev))
ggplot(cv,aes(x=size,y=dev))+
geom_point()+
geom_point(data=min, size=3, colour=alpha("red", 0.7))+
geom_line()
pruned.hitters = prune.tree(tree.hitters, best = 5)
pred.pruned = predict(pruned.hitters, test.hitters)
mean((log(test.hitters$Salary) - pred.pruned)^2)
# Bagging
# Exercise 7
#(1)
train=1:200
train.hitters=Hitters[train,]
test.hitters=Hitters[-train,]
attach(Hitters)
library(randomForest)
set.seed(0413)
bag.hitters = randomForest(log(Salary) ~ ., data = train.hitters, mtry = 19, ntree = 500,
importance = T)
#(2)
bag.pred = predict(bag.hitters, test.hitters)
mean((log(test.hitters$Salary) - bag.pred)^2)
bag.hitters = randomForest(log(Salary) ~ ., data = train.hitters, mtry = 19, ntree = 25,
importance = T)
bag.pred = predict(bag.hitters, test.hitters)
mean((log(test.hitters$Salary) - bag.pred)^2)
# Random Forest
# Exercise 8
#(1)
set.seed(0413)
rf.hitters = randomForest(log(Salary) ~ ., data = train.hitters, mtry = 6, ntree = 500,
importance = T)
rf.pred = predict(rf.hitters, test.hitters)
mean((log(test.hitters$Salary) - rf.pred)^2)
#(2)
# Variable Importance Plot
var_importance <- data.frame(variable=setdiff(colnames(Hitters), "Salary"),
importance=as.vector(rf.hitters$importance[,1]))
var_importance <- var_importance[order(var_importance[,2]),]
var_importance$variable <- factor(var_importance$variable, levels=var_importance$variable)
ggplot(var_importance,aes(x=variable, y=importance, fill=variable)) +
geom_bar(stat='identity', position='identity')+
scale_fill_discrete(name="Variable Name",guide = guide_legend(reverse=TRUE)) +
coord_flip()+
theme(axis.text.x=element_blank(),
axis.text.y=element_text(size=12),
axis.title=element_text(size=16),
plot.title=element_text(size=18),
legend.title=element_text(size=16),
legend.text=element_text(size=12))+
ylab("Mean square error")
# Boosting
# Exercise 9
#(1)
library(gbm)
set.seed(0413)
boost.hitters=gbm(log(Salary)~.,data=train.hitters,distribution="gaussian",
n.trees=5000,interaction.depth=4,shrinkage = 0.01)
summary(boost.hitters)
#(2)
yhat.boost=predict(boost.hitters,newdata=test.hitters,n.trees=5000)
mean((yhat.boost-log(test.hitters$Salary))^2)
5일 : 추천 알고리즘 (Recommendation System)
########################
#3. Colaborate filtering
#3.1 Jester5k data ####
#data loading
install.packages("recommenderlab")
library(recommenderlab)
data(Jester5k)
head(as(Jester5k,"matrix"))
image(Jester5k[1:80]) #Jester5k user1~80까지만 추출해서 ratingmatrix를 확인
#colaborate filtering
r=Recommender(Jester5k[1:1000],method="UBCF") #method를 IBCF로하면 item-based CF
pr=predict(r, Jester5k[1:2], type="ratings")
as(pr,"matrix") #이미 평가한 자료의 경우 예측값을 주지 않는다
ptype=predict(r, Jester5k[1:2], n=5) #n을 통해 추천품목수를 조정 가능
as(ptype,"list")
as(pr[1],"matrix")[,order(as(pr[1],"matrix"),decreasing = T)] #첫번째 user에 대해 예상별점을 정렬
#3.2 scale 보정 #######
#making dummy variable
a = as(Jester5k[1:1000],"matrix")
dgmat=cbind(a[1:100000],as.data.frame(cbind(rep(rownames(a),ncol(Jester5k[1:1000])),rep(colnames(a),each=nrow(Jester5k[1:1000]))))) #각 행은 하나의 평점에 대한 정보를 가짐
colnames(dgmat)=c("rating","user","item")
user=unique(dgmat$user); item=unique(dgmat$item)
head(dgmat)
dgmat=dgmat[is.na(dgmat$rating)==F,] #평점이 있는 정보만 사용
dummy=model.matrix(rating~user+item,dgmat) #user와 item에 대해 dummy variable 생성
head(dummy)
#estimate mu by ridge regression
library(glmnet)
#cv.lm=cv.glmnet(dummy, dgmat$rating, type.measure="deviance", alpha=0) #시간이 오래 걸림.. 집에서 해보세요..
lm = glmnet(dummy,dgmat$rating,family="gaussian",lambda=0.2, alpha=0) #cv돌려본 결과 최적의 lambda와 대충 일치. 실제로는 cv.lm$lambda.min을 넣어준다.
#기존 형식으로 자료 변환
dgmat$rating=dgmat$rating-dummy%*%lm$beta
user.index=match(dgmat$user,user); item.index=match(dgmat$item,item) #각 row의 user(item)가 몇번째 user(item) 인지 확인
mat=sparseMatrix(i=user.index,j=item.index,x=dgmat$rating)
mat=as.matrix(mat) ; sum(dgmat$rating==0,na.rm=T) ;mat[mat==0]=NA #dgmat$rating 값 중 0이 없음을 확인하고 NA를 넣음
colnames(mat)=item ; rownames(mat)=user
mat = as(mat,"realRatingMatrix") #Recommender 함수를 위해 다시 rating matrix로..
#예상 평점 계산
r1=Recommender(mat,method="UBCF")
pr1=predict(r1, mat[1:2], type="ratings")
pr1 = as(pr1,"matrix")
rownames(lm$beta)=gsub('user','',rownames(lm$beta)) ;rownames(lm$beta)=gsub('item','',rownames(lm$beta)) #lm$beta의 이름 정제(왜 필요한지는 하기 전에 확인)
item=as.character(item) ; user=as.character(user)
user1_dummy=c()
for(i in 1:length(item)){ #user1과 user2에 대해 앞에서 추정한 mu들을 더해주기 위한 과정
b=rep(0,length(lm$beta))
b[match(c(user[1],item[i]),rownames(lm$beta))]=1
b[1]=1
user1_dummy=cbind(user1_dummy,b)
}
scale.value=rbind(t(lm$beta) %*% user1_dummy ,t(lm$beta) %*% user2_dummy)
pr1.final= scale.value+pr1 ; colnames(pr1.final)=item
pr1.final[,order(pr1.final[1,],decreasing = T)] #user1에 대해 예상평점이 높은 순으로 정렬
#예제 - MovieLense data
data("MovieLense")
#위 과정을 따라서 해보세요!
########################
#4. matrix factorization
#4.1 Jester5k data ####
#install library
library(matfact)
mf = matfact(as(Jester5k,"matrix")[1:1000,])
pred = mf$P %*% t(mf$Q)
colnames(pred)=item ; rownames(pred)=user
head(pred)
#prediction
i1 = is.na(as(Jester5k,"matrix")[1,]) #사용자 1,2가 평점을 내리지 않은 상품들을 추출
i2 = is.na(as(Jester5k,"matrix")[2,])
pred[1,i1==T] ; sort(pred[1,i1==T],decreasing=T) #평점을 내리지 않은 상품들에 대해서만 예상평점을 관찰
pred[2,i2==T] ; sort(pred[2,i2==T],decreasing=T)
#4.2 scale 보정 ####
mat1=as(mat,"matrix") #앞에서 scale 보정에 사용했던 mat 변수를 그대로 사용하면 된다.
mf1 = matfact(mat1)
pred1 = mf1$P %*% t(mf1$Q)
colnames(pred1)=item ; rownames(pred1)=user
pred1.final=pred1[1:2,]+scale.value #앞에서와 같은 보정값을 더해주면 된다
pred1.final[1,i1==T] ; pred1.final[2,i2==T] #matrix factorization을 통한 user1,2의 예상별점
sort(pred1.final[1,i1==T],decreasing=T); sort(pred1.final[2,i2==T],decreasing=T) #matrix factorization을 통한 user1,2의 예상별점(별점 순으로 정렬)
#matrix factorization
#예제 - MovieLense data
#위 과정을 따라서 해보세요!
##########################
#7. 구매자료를 이용한 추천
#7.1 연관성분석 ####
#tot data 불러오기
install.packages("arules")
library(arules)
shop=as.matrix(shop)
rules=apriori(shop, parameter=list(supp=0.1,conf=0.9)) #연관규칙 형성
rules.sorted=sort(rules, by="lift")
inspect(rules.sorted)
#7.2 conditional regression approach ####
#너무 오래 걸림....
#pred=c()
#for(i in 1:ncol(tot)){
# cv.fit=cv.glmnet(tot[,-i],tot[,i],family="binomial",alpha=0,nfolds=3)
# pred.new=predict(cv.fit,newx=tot[,-i],s="lambda.min",type="response")
# pred=cbind(pred,pred.new)
#}
#미리 해놓은 결과값 불러오기
pred=read.csv("F:\\Stat\\job\\recommender_공개강좌\\실습자료\\prediction.csv")
rownames(pred)=pred[,1];pred=pred[,-1]
pred[which(tot==1)]=0 #구매내역이 있으면 기대확률을 볼 필요가 없으므로 0으로 처리
sort(pred[1,],decreasing=T)[1:10] #첫번째 user에게 추천하는 상품들 상위 10개와 그 기대 확률ㄷ
#인터넷 쇼핑몰 자료
#shop 자료 불러오기
#구매자료를 이용한 추천 다시 하기
관련 자료
-
강의안 및 교육 자료 참조
03. [서울대 통계연구소] R을 이용한 빅데이터 분석.zip
drive.google.com
참고 문헌
[논문]
- 없음
[보고서]
- 없음
[URL]
- 없음
문의사항
[기상학/프로그래밍 언어]
- sangho.lee.1990@gmail.com
[해양학/천문학/빅데이터]
- saimang0804@gmail.com
'프로그래밍 언어 > R' 카테고리의 다른 글
| [R] R 및 Python을 이용한 '디시인사이드' MBTI 갤러리 웹 크롤링 및 키워드 분석을 통한 워드 클라우드 생성 (1) | 2020.05.23 |
|---|---|
| [R] R을 이용한 수치해석 : 2020년 대학수학능력시험 (수능) 가형 기출문제 (0) | 2020.05.18 |
| [R] KLAPS 수치예측 모델 자료를 이용하여 내삽 방법 (Multi level B-Spline Approximation, Kriging)에 따른 전처리 및 가시화 (0) | 2020.04.14 |
| [R] R을 이용한 통계 분석 및 데이터 시각화 : dplyr (2) | 2020.04.09 |
| [R] R을 이용한 통계 분석 및 데이터 시각화 : readr (0) | 2020.04.09 |