문제

/usr/bin/x86_64-linux-gnu-ld: cannot find -lgfortran

Bioconductor를 통해 R 라이브러리를 설치할 때, 위와 같은 컴파일 오류로 인해 제대로 설치되지 않는 현상 발생

해결

https://askubuntu.com/questions/1007591/usr-bin-ld-cannot-find-lopencl

아래 위치에 gfortran 라이브러리가 없어서 문제 발생

/usr/lib/x86_64-linux-gnu/

아래 명령어로 문제 해결

sudo ln -s /usr/lib/x86_64-linux-gnu/libgfortran.so.3 /usr/lib/libgfortran.so

유명하고 많이 사용하는 R 패키지 정리

데이터 로드

RMySQL, RPostgresSQL, RSQLite - 데이터 베이스로부터 직접 데이터를 읽을 때 사용하는 패키지들이다. R[데이터베이스명] RMySQL은 MySQL의 데이터들을 직접 R로 불러올 수 있다.

XLConnect, xlsx - 이 패키지들은 Microsoft사의 엑셀을 R로부터 직접 읽어올 수 있게 한다. 물론, write도 가능하다.

foreign - SAS, SPSS 데이터셋을 읽어올 때 사용한다. 예를 들어, SAS의 경우 sas7bdat 확장자로 되어있는 파일인데, 이를 읽어올 때 foreign 패키지를 활용할 수 있다.

일반적인 텍스트 파일을 로드할 때, R에서는 별다른 패키지가 필요하지 않다. read.csv, read.table, read.fwf를 이용하면 된다. 이외의 독특한 자신만의 데이터셋을 불러오고 싶다면 CRAN guide에 데이터 import, export에 관하여 문의할 수도 있다. 

데이터 핸들링

dplyr - 데이터 subsetting, summarizing, rearranging, joining에 대한 더 쉬운 길을 제공한다. dplyr는 빠른 데이터 핸들링을 위하여 반드시 사용하는 패키지이다.(go to package) 

tidyr - 데이터셋의 레이아웃을 바꿀 때 유용한 툴이다. 데이터를 tidy format으로 바꾸기 위해 gather이나 spread 함수를 사용할 수 있다.

stringr - 문자열 다루는 것과 정규 표현식 관련 패키지이다. 

lubridate - date와 time을 더욱 다루기 쉽게 만드는 툴이다.

데이터 시각화

ggplot2 - R의 매우 유명한 시각화 패키지이다. grammar of graphics 를 활용하여 layered, customizable plot을 만들 수 있다는 장점이 있다.

ggvis - grammar of graphics을 기반으로 동작하는, 대화형, 웹베이스의 그래픽 라이브러리이다.

rgl - Interactive 3D 시각화를 위한 라이브러리이다. 

htmlwidgets - 자바스크립트 기반의 interactive 시각화를 위한 툴이다. 이 패키지는 아래와 같은 htmlwidget들을 구현하고 있다.

• leaflet (maps)
• dygraphs (time series)
• DT (tables)
• diagrammeR (diagrams)
• network3D (network graphs)
• threeJS (3D scatterplots and globes)

googleVis - R에서 데이터 시각화를 위해 구글 차트를 이용할 수 있게한다. 구글 차트 툴은 Gapminder라고 불리기도 한다. 이는 시각화 소프트웨어로, Hans Rosling의 TED 출연으로 유명해졌다.

데이터 모델링

car - car 패키지의 Anova 함수는 type2, type3 아노바 테이블을 만드는데 유명하다.

mgcv - Generalized Additive Models

lme4/nlme - Linear, Non-linear 혼합효과모형

randomForest - 머신러닝의 랜덤 포레스트

multcomp - 다중비교를 위한 툴

vcd - 범주형 변수의 시각화와 분석을 위한 툴이다.

glmnet - Lasso, elastic-net 회귀분석을 cross validation과 함께 해준다.

survival - 생존분석

caret - 회귀분석 및 분류 모델의 트레이닝을 위한 툴이다.

데이터 리포트

shiny - 인터랙티브한 웹앱을 쉽게 만들어준다. 프로그래머가 아닌 일반 사람들에게 데이터 exploring과 sharing을 쉽게 만들어준다.

R Markdown - Reproducible reporting을 위한 툴이다. R 코드를  markdown 에 쓰고, render을 하면 R Markdown은 코드를 코드 실행 결과와 함께 HTML, pdf, word 형식으로 export를 해준다. 결과를 정리하는 다른 과정 없이 자동화된 리포팅을 알 수 있게 된다. R Markdown은 RStudio와 통합된다.

xtable - xtable 함수는 R Object(예를 들어, dataframe)를 통해 latex이나 HTML 코드로 리턴해준다. 이를 통해 문서에 붙여넣기를 쉽게할 수 있다. 

공간 데이터 

sp, maptools - shapefile을 비롯한 공간 데이터를 로딩할 수 있는 툴이다.

maps -  맵에 다각형을 쉽게 그려주는 툴이다.

ggmap -  Google map으로 부터 street map을 다운로드 받고, ggplot의 background로 쓸 수 있다.

시계열, 금융 데이터

zoo - 시계열 데이터를 저장하기 위한 가장 유명한 포맷을 다룰 수 있다.

xts - 시계열 데이터를 다루기 위한 툴

quantmod - 금융 데이터를 다운로드하고, 그래프를 그리고, 분석할 수 있는 툴이다.

높은 성능을 내기 위한 R 코드 작성

Rcpp - C++을 call하는 R function을 사용한다. 

data.table - 데이터셋을 빠르게 조작하기 위한 다른 방법을 사용한다. "빅 데이터"에 유용하다.

parallel - 큰 데이터 셋을 다루기 위한 병렬 프로세싱을 사용한다.

웹으로 작업하기

XML - R을 통해 XML 문서를 읽고 만드는 패키지

jsonlite - R을 통해 JSON 데이터를 읽고 만드는 패키지

httr - HTTP Connection을 통한 작업을 위한 라이브러리

R 패키지 만들기

devtools - 코드를 R 패키지로 만들기

testthat - 프로젝트의 유닛 테스트를 위한 쉬운 방법을 제공한다.

roxygen2 - R 패키지의 도큐먼트를 만들기 위한 빠른 방법. roxygen2는 코드의 코멘트를 도큐먼트로 만들고, 패키지의 네임스페이스를 만든다.

문자 변수를 날짜 변수로 변환하기

R에서 시계열 데이터를 다룰 때, 변수를 날짜 타입으로 변환한다면 플롯을 하는데 간편한 점이 많다. R의 날짜 데이터 타입으로는 "POSIXct" 가 있으며, as.POSIXct 함수를 통해 문자 변수를 POSIXct 타입으로 변환할 수 있다. 반대로 POSIXct 타입을 문자 타입으로 바꾸는 것은 format이나 as.character를 이용한다.  

# create POSIXct variables
as.POSIXct("080406 10:11", format = "%y%m%d %H:%M")
## [1] "2008-04-06 10:11:00 CDT"

as.POSIXct("2008-04-06 10:11:01 PM", format = "%Y-%m-%d %I:%M:%S %p")
## [1] "2008-04-06 22:11:01 CDT"

as.POSIXct("08/04/06 22:11:00", format = "%m/%d/%y %H:%M:%S")
## [1] "2006-08-04 22:11:00 CDT"
# convert POSIXct variables to character strings

format(as.POSIXct("080406 10:11", format = "%y%m%d %H:%M"), "%m/%d/%Y %I:%M %p")
## [1] "04/06/2008 10:11 AM"

as.character(as.POSIXct("080406 10:11", format = "%y%m%d %H:%M"), format = "%m-%d-%y %H:%M")
## [1] "04-06-08 10:11"

만약 어떤 dataframe의 column이 POSIXct 타입이라면, 내부적으로는 정수를 저장하고 있기 때문에 Plot을 할 때는 숫자 형태로 변환하여 데이터를 순차적으로 나타내준다. 

그래프 그리기

data <- read.csv("C:\\workspace\\r\\data\\data.csv")
plot(data$dateFormat, data$value)

이렇게 기본 내장함수인 plot을 사용해도 x축을 시간으로 하여 괜찮은 모양으로 그려주는 것을 확인할 수 있다. 만약 x축 값의 포맷을 지정하거나, 적절한 간격을 주고 싶다면 ggplot2 패키지를 이용하는 것도 좋다.

GGPLOT2 라이브러리를 통한 그래프 그리기

library(scales)

ggplot(data = data, aes(x = date, y = value)) + 
    geom_line(size=2) + ylim(0,1000) + xlab("Time") + ylab("㎍/㎥") + 
    scale_x_datetime(breaks = date_breaks("1 hour"), minor_breaks=date_breaks("1 hour"), labels=date_format("%H:%M"), limits=xlim) +
    ggtitle("title") + theme_bw() + theme(text = element_text(size = 25))

ggplot 명령어 설명 

aes : (x축, y축에 사용될 변수를 지정한다.)

geom_line : line plot. size 키워드를 통해 선의 굵기 지정

ylim : y축의 범위를 지정함

xlab : x축 label

ylab : y축 label

scale_x_datetime : x축의 간격, 포맷을 지정함. (library(scales) 이후에 사용할 수 있다.)

ggtitle : title 지정

theme_bw() : 해당 테마 지정 (black & white style)

theme : 그래프의 전반적인 속성을 지정한다. 위에서는 텍스트 사이즈를 25 조정함

특정 날짜 범위에 있는 데이터만 골라내기

아래 코드는 date가 2018/01/16 인 row만 골라서 subset으로 저장한다. 이처럼 format을 통해 POSIXct 데이터 타입을 character 타입으로 변환할 수 있다.

subset <- data[(format(data$date, "%m/%d/%Y") == "01/16/2018")]

참고

http://biostat.mc.vanderbilt.edu/wiki/pub/Main/ColeBeck/datestimes.pdf

R - 특정 버전의 패키지 설치하기

우선 패키지의 버전을 알아보는 방법은 R을 실행 시킨 후, 

library(package_name)

sessionInfo()

를 입력하면 패키지의 버전을 알 수 있다.

특정 버전의 R 패키지를 설치하는 방법은 우선, Install.packages("devtools")로 devtool을 설치한다.

require(devtools)
install_version("ggplot2", version = "0.9.1", repos = "http://cran.us.r-project.org")

위와 같은 방식으로 특정 버전의 CRAN 패키지를 설치할 수 있다. 이미 패키지가 설치된 상태에서도 install_version 을 활용하면 버전이 바뀌어 다시 설치되기 때문에, 패키지를 다운그레이드 시킬 때도 유용하게 쓸 수 있다.

R에서 csv 파일을 읽어들일 때 UTF-8 포맷의 텍스트 파일을 읽어오지 못하는 문제가 있었다.

temp <- read.table(path, sep=",", stringsAsFactors=FALSE, header=FALSE, nrows=10)
head(temp)

  V1                                                        V2     V3 V4 V5 V6 V7 V8

1  4 01. 16. 18 \xec삤\xec쟾 08\xec떆 02遺\x84 25珥\x88,26.188 17.695 NA NA NA NA NA

2  5 01. 16. 18 \xec삤\xec쟾 08\xec떆 03遺\x84 25珥\x88,26.456 22.904 NA NA NA NA NA

3  6 01. 16. 18 \xec삤\xec쟾 08\xec떆 04遺\x84 25珥\x88,26.530 19.185 NA NA NA NA NA

4  7 01. 16. 18 \xec삤\xec쟾 08\xec떆 05遺\x84 25珥\x88,26.603 19.383 NA NA NA NA NA

5  8 01. 16. 18 \xec삤\xec쟾 08\xec떆 06遺\x84 25珥\x88,26.676 21.714 NA NA NA NA NA

6  9 01. 16. 18 \xec삤\xec쟾 08\xec떆 07遺\x84 25珥\x88,26.652 21.107 NA NA NA NA NA

위와 같이 파일이 깨지는 현상

temp <- read.table(path, sep=",", stringsAsFactors=FALSE, header=FALSE, nrows=10,fileEncoding = "utf-8")

Error in read.table(paste0(path, date, name, filename_hobo), skip = 5,  : 

  입력에 가능한 라인들이 없습니다

In addition: Warning messages:

1: In readLines(file, skip) :

  invalid input found on input connection 'C:\workspace\r\data\mongol\180116\1\MG2018_0116_KSY_HOBO.csv'

2: In readLines(file, skip) :

  incomplete final line found on 'C:\workspace\r\data\mongol\180116\1\MG2018_0116_KSY_HOBO.csv'

이후에 읽으니까 파일 읽기에 성공하였다.

R - 종속변수가 순서 척도인 로지스틱 회귀분석

다항 로지스틱 회귀분석에서 종속변수가 명목형이 아닌 순서형인 경우 분석하는 방법에 관한 포스팅입니다. 이 포스팅은 이 사이트의 튜토리얼을 보고 작성하였습니다.

라이브러리 로드

require(foreign)
require(ggplot2)
require(MASS)
require(Hmisc)
require(reshape2)

데이터 읽기

dat <- read.dta("https://stats.idre.ucla.edu/stat/data/ologit.dta")
head(dat)

> head(dat)
            apply pared public  gpa
1     very likely     0      0 3.26
2 somewhat likely     1      0 3.21
3        unlikely     1      1 3.94
4 somewhat likely     0      0 2.81
5 somewhat likely     0      0 2.53
6        unlikely     0      1 2.59

X : pared(부모중 한명이 대학 학위가 있는지 여부), public(졸업한 학교가 공립인지), gpa(학점)

Y : apply (“unlikely”, “somewhat likely”, “very likely”, 각각 1, 2, 3으로 코딩되어있음) => 대학원에 진학할 가능성

이 때 Y는 순서형 척도이다.

범주형 변수의 빈도 확인

lapply 함수를 이용해 범주형 변수의 빈도표를 확인

## one at a time, table apply, pared, and public
lapply(dat[, c("apply", "pared", "public")], table)

> lapply(dat[, c("apply", "pared", "public")], table)
$apply

       unlikely somewhat likely     very likely 
            220             140              40 

$pared

  0   1 
337  63 

$public

  0   1 
343  57 

범주형 독립변수와 종속변수 관계 파악 : 3-변수 빈도표

변수 빈도표를 통해 로지스틱 회귀분석을 하기에 앞서 예상되는 결과를 미리 확인한다.

## three way cross tabs (xtabs) and flatten the table
ftable(xtabs(~ public + apply + pared, data = dat))

> ftable(xtabs(~ public + apply + pared, data = dat))
                       pared   0   1
public apply                        
0      unlikely              175  14
       somewhat likely        98  26
       very likely            20  10
1      unlikely               25   6
       somewhat likely        12   4
       very likely             7   3

연속형 변수 분포 확인

summary(dat$gpa)

> summary(dat$gpa)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.900   2.720   2.990   2.999   3.270   4.000 

sd(dat$gpa)

> sd(dat$gpa)
[1] 0.3979409

연속형 독립 변수와 종속 변수와의 관계 파악

Apply에 따라 gpa가 어떻게 분포하여있는지 상자그림을 통해 확인한다. 이러한 경향을 수치화하는 것이 로지스틱 회귀분석의 목표이다.

ggplot(dat, aes(x = apply, y = gpa)) +
  geom_boxplot(size = .75) +
  geom_jitter(alpha = .5) +
  facet_grid(pared ~ public, margins = TRUE) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))

지금까지는 데이터에 대해서 살펴보았고, 본격적으로 순서 척도가 종속변수인 로지스틱 회귀분석을 어떻게 하는지 알아본다.

이 때 가능한 분석방법은 꼭 순서형 척도만은 아닐 것이다. GPA만을 독립변수로하여 ANOVA를 할 수도 있고, y를 1,2,3으로 놓고 다중회귀분석을 할 수도 있다. 또 y를 순서가 없는 명목형 척도로 보고 명목형 다항 로지스틱 분석(multinomial logistic regression)을 할 수도 있을 것이다. 하지만 해당 데이터에 대해 위 분석 방법들은 가정에 맞지 않는다. 따라서 순서형 척도를 종속변수로하는 로지스틱 회귀분석(ordinal logistic regression)이 적합하다.

순서형 척도를 종속변수로하는 로지스틱 회귀분석의 중요한 가정은 회귀계수가 같다고 가정하는 것이다. y=1,2,3인 이  데이터의 경우 어떠한 변수의 변화가 1 vs 2,3의 OR에 미치는 영향과 1,2 vs 3의 OR에 미치는 영향이 같다고 가정한다. y의 범주가 3개인 경우 2개의 식이 나오면 이 두식은 회귀계수가 같고, 단지 절편만 다르다.

MASS 패키지의 polr 함수를(proportional odds logistic regession) 이용하여 순서형 변수를 종속변수로하는 다항 로지스틱 회귀분석을 해보자.

회귀식을 적어주고 summary를 통해 적합된 결과를 본다. Hess=TRUE는 optimization을 할 때 사용하는 Hessian matrix를 의미한다. (https://en.wikipedia.org/wiki/Hessian_matrix)

## fit ordered logit model and store results 'm'
m <- polr(apply ~ pared + public + gpa, data = dat, Hess=TRUE)

## view a summary of the model
summary(m)

> summary(m)
Call:
polr(formula = apply ~ pared + public + gpa, data = dat, Hess = TRUE)

Coefficients:
          Value Std. Error t value
pared   1.04769     0.2658  3.9418
public -0.05879     0.2979 -0.1974
gpa     0.61594     0.2606  2.3632

Intercepts:
                            Value   Std. Error t value
unlikely|somewhat likely     2.2039  0.7795     2.8272
somewhat likely|very likely  4.2994  0.8043     5.3453

Residual Deviance: 717.0249 
AIC: 727.0249 

결과 해석

회귀계수에 exponential을 취하면 OR(Odds ratio)이 되며 이를 통해 결과를 해석해본다.

## odds ratios
exp(coef(m))

> exp(coef(m))
    pared    public       gpa 
2.8510579 0.9429088 1.8513972 

순서척도 로지스틱 회귀분석에서 OR은 proportional OR이라고도 부른다.

• 다른 변수가 고정되어있을 때, pared=0일 때에 비하여 pared=1일 때, “very likely” vs “somewhat likely” or “unlikely”의 OR이 2.851이다.
  
• 다른 변수가 고정되어있을 때, pared=0일 때에 비하여 pared=1일 때, “very likely” or “somewhat likely” vs “unlikely”의 OR이 2.851이다.
• 이는 즉, pared=1일 때, 대학원 진학(apply)을 할 가능성이 높다는 것을 뜻한다.
• 같은 binary 변수인 public도 위와 마찬가지로 해석한다.
• 다른 변수가 고정되어있을 때, gpa가 1 증가할 때, “very likely” vs “somewhat likely” or “unlikely” Odds가 1.85배가 된다.
  
• 다른 변수가 고정되어있을 때, gpa가 1 증가할 때, “very likely” or “somewhat likely” vs “unlikely” Odds가 1.85배가 된다.

이를 식으로 표현하면

log Odds(Y=very) = 2.2039 + 1.04pared - 0.05879public + 0.61594gpa

log Odds(Y=very or somewhat) = 4.2994 + 1.04pared - 0.05879public + 0.61594gpa

• 이 때, Odds(Y=1) = P(P1/(1-P1)), Odds(Y=1,2) = P((P1+P2)/(1-P1-P2))
• log Odds 는 logit function 이라고도 부른다.
  

위 두 방정식을 풀면 P(Y=1), P(Y=2), P(Y=3) 3개의 식을 얻을 수 있다. 이 식은 X가 주어졌을 때, Y의 예측값을 구한다.

R - Multinomial Logistic Regression

Multinomial Logistic Regression이란 y의 범주가 3개 이상(multi)이며 명목형(nomial)일 때 사용하는 로지스틱 회귀분석이다. 본 포스팅은 이 사이트의 튜토리얼을 보고 실습하고 해석한 자료이다.

데이터 설명

해당 데이터셋은 hsbdemo 데이터셋으로 총 n은 200이다.

ml <- read.dta("https://stats.idre.ucla.edu/stat/data/hsbdemo.dta")

> head(ml)
   id female    ses schtyp     prog read write math science socst       honors awards cid
1  45 female    low public vocation   34    35   41      29    26 not enrolled      0   1
2 108   male middle public  general   34    33   41      36    36 not enrolled      0   1
3  15   male   high public vocation   39    39   44      26    42 not enrolled      0   1
4  67   male    low public vocation   37    37   42      33    32 not enrolled      0   1
5 153   male middle public vocation   39    31   40      39    51 not enrolled      0   1
6  51 female   high public  general   42    36   42      31    39 not enrolled      0   1

Y : prog (프로그램 타입)

X : ses(social economic status), write (writing score)

다중회귀분석을 통해 X가 Y에 어떤 영향을 미치는지 알아본다.

빈도 테이블과 그룹 평균을 통해 X와 Y의 관계 알아보기

with(ml, table(ses, prog)) 
# 또는 table(ml$ses, ml$prog) 
# 또는 xtabs(~ses+prog, data=ml)

> with(ml, table(ses, prog))
        prog
ses      general academic vocation
  low         16       19       12
  middle      20       44       31
  high         9       42        7

with(ml, do.call(rbind, tapply(write, prog, function(x) c(M = mean(x), SD = sd(x)))))

> with(ml, do.call(rbind, tapply(write, prog, function(x) c(M = mean(x), SD = sd(x)))))
                M       SD
general  51.33333 9.397775
academic 56.25714 7.943343
vocation 46.76000 9.318754

nnet 패키지의 multinom 함수를 통해 다중 로지스틱 회귀분석 실시

ml$prog2 <- relevel(ml$prog, ref = "academic")
test <- multinom(prog2 ~ ses + write, data = ml)

multinom을 실행하기 전에 relevel이라는 함수를 통해 reference를 지정해준다. 이 reference를 기준으로 결과를 해석할 수 있다. 따라서 이 reference는 어떤 baseline이 되어야한다. 이 경우에 academic이 baseline이며 이를 기준으로 하여 분석한다.

multinom 함수를 실행하면 모델이 fitting되어 계수가 결정되고 이 정보는 test에 담긴다.

summary(test)

> summary(test)
Call:
multinom(formula = prog2 ~ ses + write, data = ml)

Coefficients:
         (Intercept)  sesmiddle    seshigh      write
general     2.852198 -0.5332810 -1.1628226 -0.0579287
vocation    5.218260  0.2913859 -0.9826649 -0.1136037

Std. Errors:
         (Intercept) sesmiddle   seshigh      write
general     1.166441 0.4437323 0.5142196 0.02141097
vocation    1.163552 0.4763739 0.5955665 0.02221996

Residual Deviance: 359.9635 
AIC: 375.9635 

이 예제의 경우 다중 로지스틱 회귀분석의 계수는 2세트가 나오게 된다.

1. academic vs general 의 log odds에 관한 계수

2. academic vs vocation의 log odds에 관한 계수

결과 해석

• write가 1단위 증가할 때, academic이 general이 될 log odds가 -0.058이 된다.
• write가 1단위 증가할 때, academic이 vocation이 될 log odds가 -0.1136037이 된다.
• ses:row에서 ses:high로 변할 때, academic이 general이 될 log odds가 -1.1628이 된다.
• ses:row에서 ses:middle로 변할 때, academic이 general이 될 log odds가 -0.5332가 된다.
• ses:row에서 ses:high로 변할 때, academic이 vocation이 될 log odds가 -0.9826이 된다.
• ses:row에서 ses:middle로 변할 때, academic이 vocation이 될 log odds가 +0.29가 된다.
  

변수의 유의성 판단

coefficient가 0라는 귀무가설 하에서 coef ~ N(0, sd^2)라는걸 이용하여 z-test 한다.

z <- summary(test)$coefficients/summary(test)$standard.errors
z

> z <- summary(test)$coefficients/summary(test)$standard.errors
> z
         (Intercept)  sesmiddle   seshigh     write
general     2.445214 -1.2018081 -2.261334 -2.705562
vocation    4.484769  0.6116747 -1.649967 -5.112689

# 2-tailed z test
p <- (1 - pnorm(abs(z), 0, 1)) * 2
p

> p <- (1 - pnorm(abs(z), 0, 1)) * 2
> p
          (Intercept) sesmiddle    seshigh        write
general  0.0144766100 0.2294379 0.02373856 6.818902e-03
vocation 0.0000072993 0.5407530 0.09894976 3.176045e-07

OR 추정값

exp(coef(test))

> exp(coef(test))
         (Intercept) sesmiddle   seshigh     write
general     17.32582 0.5866769 0.3126026 0.9437172
vocation   184.61262 1.3382809 0.3743123 0.8926116

회귀 계수들에 exponential을 취하면, 이는 OR에 대한 추정값이 된다.

Write가 평균일 때, ses에 따라 예측값이 어떻게 달라지는가?

dses <- data.frame(ses = c("low", "middle", "high"), write = mean(ml$write))
dses
predict(test, newdata = dses, "probs")

> dses <- data.frame(ses = c("low", "middle", "high"), write = mean(ml$write))
> dses
     ses  write
1    low 52.775
2 middle 52.775
3   high 52.775

> predict(test, newdata = dses, "probs")
   academic   general  vocation
1 0.4396845 0.3581917 0.2021238
2 0.4777488 0.2283353 0.2939159
3 0.7009007 0.1784939 0.1206054

Write를 변화시켜가면서 예측값의 변화 관찰하기

dwrite <- data.frame(ses = rep(c("low", "middle", "high"), each = 41), write = rep(c(30:70),
                                                                                   3))

## store the predicted probabilities for each value of ses and write
pp.write <- cbind(dwrite, predict(test, newdata = dwrite, type = "probs", se = TRUE))

## calculate the mean probabilities within each level of ses
by(pp.write[, 3:5], pp.write$ses, colMeans)

## melt data set to long for ggplot2
lpp <- melt(pp.write, id.vars = c("ses", "write"), value.name = "probability")
head(lpp)  # view first few rows

> head(lpp)  # view first few rows
  ses write variable probability
1 low    30 academic  0.09843588
2 low    31 academic  0.10716868
3 low    32 academic  0.11650390
4 low    33 academic  0.12645834
5 low    34 academic  0.13704576
6 low    35 academic  0.14827643

## plot predicted probabilities across write values for each level of ses
## facetted by program type
ggplot(lpp, aes(x = write, y = probability, colour = ses)) + geom_line() + facet_grid(variable ~
                                                                                        ., scales = "free")

• 이 그래프는 식을 전개하여 y를 sigmoid 함수 형태로 만든 것이다. 2개의 적합식에서 3개의 함수를 만들어낼 수 있다.

 Y~write의 결과랑 무엇이 다른가?

예측의 관점에서 보면 ses를 빼고 write만을 통해 Y를 예측하면 예측을 세분화하지 못하여 예측력이 떨어진다고 볼 수 있다. 이는 ses를 고려하지 않고, write와 Y의 관계만을 보기 때문이다. 따라서 ses가 추가됨으로써 write의 계수가 바뀐다는 것을 알 수 있다. 이 과정을 ses를 보정한다고도 말한다.

출처 - https://stats.idre.ucla.edu/r/dae/multinomial-logistic-regression/

R - 로지스틱 회귀분석

데이터 탐색

GRE, GPA, RANK이 입학(admission)에 어떤 영향을 주는지 로지스틱 회귀분석을 통해 분석한다.

library(aod)
library(ggplot2)

mydata <- read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv")

# view the first few rows of the data
head(mydata)

# 데이터의 대략적인 분포 확인
summary(mydata)

# 데이터 구조 확인
str(mydata)

# 변수별 표준편차 확인
sapply(mydata, sd)

# contingency table : xtabs(~ x + y, data)
xtabs(~admit+rank, data=mydata)

> head(mydata)
  admit gre  gpa rank
1     0 380 3.61    3
2     1 660 3.67    3
3     1 800 4.00    1
4     1 640 3.19    4
5     0 520 2.93    4
6     1 760 3.00    2

> summary(mydata)
     admit             gre             gpa        rank   
 Min.   :0.0000   Min.   :220.0   Min.   :2.260   1: 61  
 1st Qu.:0.0000   1st Qu.:520.0   1st Qu.:3.130   2:151  
 Median :0.0000   Median :580.0   Median :3.395   3:121  
 Mean   :0.3175   Mean   :587.7   Mean   :3.390   4: 67  
 3rd Qu.:1.0000   3rd Qu.:660.0   3rd Qu.:3.670          
 Max.   :1.0000   Max.   :800.0   Max.   :4.000  

> str(mydata)
'data.frame':	400 obs. of  4 variables:
 $ admit: int  0 1 1 1 0 1 1 0 1 0 ...
 $ gre  : int  380 660 800 640 520 760 560 400 540 700 ...
 $ gpa  : num  3.61 3.67 4 3.19 2.93 3 2.98 3.08 3.39 3.92 ...
 $ rank : Factor w/ 4 levels "1","2","3","4": 3 3 1 4 4 2 1 2 3 2 ...

> sapply(mydata, sd)
      admit         gre         gpa        rank 
  0.4660867 115.5165364   0.3805668   0.9444602 

> xtabs(~admit+rank, data=mydata)
     rank
admit  1  2  3  4
    0 28 97 93 55
    1 33 54 28 12

GLM을 통한 로지스틱 회귀모형 구축

rank 변수를 factor 타입으로 변경시킨다.

glm 을통해 로지스틱회귀모형을 fitting시킨다. family='binomial' 인자를 통해, glm으로 로지스틱 회귀모형을 쓸 수 있다. (link function이 logit function이라는 의미)

mydata$rank <- factor(mydata$rank)
mylogit <- glm(admit ~ gre + gpa + rank, data = mydata, family = "binomial")

summary(mylogit)

> summary(mylogit)

Call:
glm(formula = admit ~ gre + gpa + rank, family = "binomial", 
    data = mydata)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6268  -0.8662  -0.6388   1.1490   2.0790  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.989979   1.139951  -3.500 0.000465 ***
gre          0.002264   0.001094   2.070 0.038465 *  
gpa          0.804038   0.331819   2.423 0.015388 *  
rank2       -0.675443   0.316490  -2.134 0.032829 *  
rank3       -1.340204   0.345306  -3.881 0.000104 ***
rank4       -1.551464   0.417832  -3.713 0.000205 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 499.98  on 399  degrees of freedom
Residual deviance: 458.52  on 394  degrees of freedom
AIC: 470.52

Number of Fisher Scoring iterations: 4

로지스틱 회귀모형 결과 해석

Call : 구축한 모형에 대해 다시 상기시켜준다.

Deviance Residuals : Deviance residual에 대한 정보를 알려주는데, model fitting이 잘 되었는지에 대한 measure이다. 이를 통해 모델이 잘 적합됐는지를 평가할 수 있다.

Coefficient : 회귀계수와 그것들의 표준편차, z-statistics(wals's z-statistics), p-value를 나타낸다. 위 결과에서는 모든 변수가 유의한 것을 알 수 있다. 로지스틱 회귀모형에서는 회귀계수가 변수가 한 단위 증가했을 때 log(odds)의 증가량으로 해석할 수 있다.

example

gre가 1증가할 때, admission의 log odds(non-admission에 대한)가 0.0022 증가한다.

gpu의 경우도 마찬가지로 해석한다.

더미변수인 rank의 경우 해석이 약간 다른데, 예를 들어 rank2의 회귀계수 -0.67은 rank1에서 rank2로 바뀌었을 때, log(odds)의 변화량이다. 즉, rank1에 비해 rank2가 admission에 안좋은 영향을 준다는 것을 알 수 있다.

회귀계수들의 신뢰구간 얻기

Log-likelihood를 통해 구하는 법

confint(mylogit)

> confint(mylogit)
Waiting for profiling to be done...
                    2.5 %       97.5 %
(Intercept) -6.2716202334 -1.792547080
gre          0.0001375921  0.004435874
gpa          0.1602959439  1.464142727
rank2       -1.3008888002 -0.056745722
rank3       -2.0276713127 -0.670372346
rank4       -2.4000265384 -0.753542605

주어진 회귀계수의 표준편차를 이용해 신뢰구간을 구하는 법

## CIs using standard errors
confint.default(mylogit)

> confint.default(mylogit)
                    2.5 %       97.5 %
(Intercept) -6.2242418514 -1.755716295
gre          0.0001202298  0.004408622
gpa          0.1536836760  1.454391423
rank2       -1.2957512650 -0.055134591
rank3       -2.0169920597 -0.663415773
rank4       -2.3703986294 -0.732528724

이 신뢰구간은 이렇게도 계산 가능하다.예를 들어, gre의 95% 신뢰구간은 아래와 같이 구한다.

> 0.002264-1.96*0.001094
[1] 0.00011976
> 0.002264+1.96*0.001094
[1] 0.00440824

Wald Test를 통해 범주형 변수의 overall effect 파악하기

rank 변수를 dummy 변수로 만들어, reference인 rank1과 비교해 어떤 effect가 있는지 분석하였다. 하지만 aod 패키지의 walt.test를 이용하면, rank의 overall effect를 파악할 수 있다. [b:회귀계수, Sigma:error term의 공분산행렬, Terms:rank변수가 있는 열] wald.test는 순서를 통해 해당 범주형 변수의 위치를 파악하기 때문에 순서를 잘 신경 써야 한다.

wald.test(b = coef(mylogit), Sigma = vcov(mylogit), Terms = 4:6)

> wald.test(b = coef(mylogit), Sigma = vcov(mylogit), Terms = 4:6)
Wald test:
----------

Chi-squared test:
X2 = 20.9, df = 3, P(> X2) = 0.00011

카이스퀘어 테스트의 p-value가 0.00011이므로 rank 변수는 유의하다.

변수간의 회귀계수가 같은지 검정하기

wald test를 통해 rank가 admission에 유의한 효과가 있다는 것을 파악하였고, rank1이라는 reference에 비해 rank2와 rank2가 유의한 영향을 준다는 것도 확인하였다. 한가지 의문은 rank2, rank3, rank4 의 회귀계수가 동일한지에 대한 것이다. 회귀 계수가 동일하다면, rank1일 때에 비해 rank2, rank3, rank4의 효과가 동일하다는 것을 의미한다.

이 구문은 rank2의 효과와 rank3의 효과가 동일한지를 검정한다.

l <- cbind(0, 0, 0, 1, -1, 0)
wald.test(b = coef(mylogit), Sigma = vcov(mylogit), L = l)

> l <- cbind(0, 0, 0, 1, -1, 0)
> wald.test(b = coef(mylogit), Sigma = vcov(mylogit), L = l)
Wald test:
----------

Chi-squared test:
X2 = 5.5, df = 1, P(> X2) = 0.019

검정 결과, rank2, rank3의 효과는 유의하다 다른 것으라 나타났다.

OR과 회귀계수의 관계

앞서 로지스틱 회귀분석에서의 회귀계수는 log odds의 증가량이라고 언급하였다. 그러면 회귀계수에 exponential을 취하면, OR의 증분이 된다.

exp(coef(mylogit))

> exp(coef(mylogit))
(Intercept)         gre         gpa       rank2       rank3       rank4 
  0.0185001   1.0022670   2.2345448   0.5089310   0.2617923   0.2119375 

해석 ex) gpa가 1증가하면 admission의 non-admission에 대한 OR이 2.23배 증가한다.

예측하기

1. 다른 변수들을 고정하고 rank가 변할 때 예측값의 변화를 보기

newdata1 <- with(mydata, data.frame(gre = mean(gre), gpa = mean(gpa), rank = factor(1:4)))

> newdata1
    gre    gpa rank
1 587.7 3.3899    1
2 587.7 3.3899    2
3 587.7 3.3899    3
4 587.7 3.3899    4

newdata1$rankP <- predict(mylogit, newdata = newdata1, type = "response") # probability 점추정값
newdata1

> newdata1
    gre    gpa rank     rankP
1 587.7 3.3899    1 0.5166016
2 587.7 3.3899    2 0.3522846
3 587.7 3.3899    3 0.2186120
4 587.7 3.3899    4 0.1846684

2. rank, gpa를 고정한 후, gre의 효과 보기

newdata2 <- with(mydata, data.frame(gre = rep(seq(from = 200, to = 800, length.out = 100),
                                              4), gpa = mean(gpa), rank = factor(rep(1:4, each = 100))))
newdata2

> newdata2
         gre    gpa rank
1   200.0000 3.3899    1
2   206.0606 3.3899    1
3   212.1212 3.3899    1
4   218.1818 3.3899    1
5   224.2424 3.3899    1
6   230.3030 3.3899    1
7   236.3636 3.3899    1
8   242.4242 3.3899    1
9   248.4848 3.3899    1
10  254.5455 3.3899    1

newdata3 <- cbind(newdata2, predict(mylogit, newdata = newdata2, type = "link", 
                                    se = TRUE))
newdata3 <- within(newdata3, {
  PredictedProb <- plogis(fit)
  LL <- plogis(fit - (1.96 * se.fit))
  UL <- plogis(fit + (1.96 * se.fit))
})

> head(newdata3)
       gre    gpa rank        fit    se.fit residual.scale        UL        LL PredictedProb
1 200.0000 3.3899    1 -0.8114870 0.5147714              1 0.5492064 0.1393812     0.3075737
2 206.0606 3.3899    1 -0.7977632 0.5090986              1 0.5498513 0.1423880     0.3105042
3 212.1212 3.3899    1 -0.7840394 0.5034491              1 0.5505074 0.1454429     0.3134499
4 218.1818 3.3899    1 -0.7703156 0.4978239              1 0.5511750 0.1485460     0.3164108
5 224.2424 3.3899    1 -0.7565919 0.4922237              1 0.5518545 0.1516973     0.3193867
6 230.3030 3.3899    1 -0.7428681 0.4866494              1 0.5525464 0.1548966     0.3223773

이 때 type="link" 옵션을 주면, link scale로 예측값을 내준다. 이는 logit 변환을 하기전의 예측값이다. 이를 통해 예측한 확률의 신뢰구간을 구할 수 있다. 왜냐하면 이 link scale의 예측값 (fit 변수)의 표준편차를 쉽게 알 수 있기 때문이다. 따라서 link scale의 예측값의 신뢰구간을 구한후 이를 다시 logit 변환하여 예측 확률의 신뢰구간을 얻을 수 있다. 이를 구현한 것이 바로 위의 코드이다.

plogis 함수를 link scale의 예측값에서 확률 추정값 p를 구해준다. 즉 log(p/1-p) = Xb 의 방정식을 풀어서 p를 구한다. 예를 들어 첫번째 라인의 경우, log(p/(1-p)) = -0.811 의 방정식을 양변에 exponential을 취해 풀면, 대략 p=0.308이 나온다.

newdata3의 예측값 시각화하기

rank를 고정시켜놓고 gre가 admission 확률에 어떤 영향을 주는지를 시각화

ggplot(newdata3, aes(x = gre, y = PredictedProb)) + geom_ribbon(aes(ymin = LL,
                                                                    ymax = UL, fill = rank), alpha = 0.2) + geom_line(aes(colour = rank),
                                                                                                                      size = 1)

출처

https://stats.idre.ucla.edu/r/dae/logit-regression/

원래 데이터

해당 지역의 10000명당 기관의 수를 나타낸 데이터입니다. (출처-통계청)

           sgg_kor facility_type 2003fac 2004fac 2005fac 2006fac 2007fac
 4               계  상급종합병원  0.0086  0.0086 0.00852  0.0087  0.0087
6               계      종합병원  0.0494  0.0494 0.05054  0.0510  0.0530
8               계          병원  0.1645  0.1755 0.18450  0.1937  0.2127
10              계      요양병원  0.0139  0.0231 0.04120  0.0727  0.1200
14              계      치과병원  0.0211  0.0221 0.02517  0.0274  0.0311
18              계      한방병원  0.0311  0.0320 0.03024  0.0292  0.0288
22              계        조산원  0.0143  0.0129 0.01055  0.0103  0.0104
24              계        보건소  0.0473  0.0475 0.04750  0.0478  0.0481
26              계      보건지소  0.2601  0.2605 0.25859  0.2571  0.2602
28              계    보건진료소  0.3865  0.3881 0.38504  0.3843  0.3875
30              계    보건의료원  0.0035  0.0035 0.00345  0.0034  0.0035
36            서울  상급종합병원  0.0195  0.0195 0.01942  0.0193  0.0196
38            서울      종합병원  0.0418  0.0409 0.04079  0.0406  0.0412
40            서울          병원  0.1158  0.1168 0.12431  0.1304  0.1432
42            서울      요양병원  0.0029  0.0058 0.01068  0.0251  0.0461

Melt

Melt는 Column을 Row로 바꾸는 함수입니다. (Key=id에 지정, 나머지값들은 Row로 변환)

facility <- melt(facility, id=c("sgg_kor", 'facility_type'))

sgg_kor facility_type variable  value
 1      계  상급종합병원  2003fac 0.0086
 2      계      종합병원  2003fac 0.0494
 3      계          병원  2003fac 0.1645
 4      계      요양병원  2003fac 0.0139
 5      계      치과병원  2003fac 0.0211
 6      계      한방병원  2003fac 0.0311

reshape

Reshape함수를 통해 Row를 Column으로 바꿀 수 있습니다. (Key를 idvar로 놓고, 컬럼으로 바꾸고 싶은 변수를 timevar에 지정)

facility <- reshape(facility, idvar=c("sgg_kor", "variable"), timevar="facility_type", direction="wide")

  sgg_kor variable value.상급종합병원 value.종합병원 value.병원 value.요양병원 value.치과병원      
계         2003fac             0.0086         0.0494     0.1645         0.0139         0.02      
    서울  2003fac             0.0195         0.0418     0.1158         0.0029         0.0350         
  강남구  2003fac             0.0373         0.0560     0.2425             NA         0.1679         
  강동구  2003fac             0.0209         0.0417     0.1043             NA             NA         
  강서구  2003fac                 NA         0.0185     0.1668         0.0185         0.0185         
  관악구  2003fac                 NA         0.0190     0.0949             NA             NA             

R 랜덤포레스트 사용하기

머신러닝 분야에서 많이 쓰이는 예제 데이터셋 Pima Diabetes 자료를 통해 간단한 Random Forest를 사용하는 방법을 알아보겠습니다.

라이브러리 로드

library(MASS)
library(randomForest)
library(caret)

MASS : 예제 데이터셋이 있는 라이브러리

randomForest : randomForest 이용을 위한 라이브러리

caret : 성능 평가시 confusion matrix를 그리기 위한 라이브러리

데이터 불러오기

data("Pima.tr")
data("Pima.te")

R에서 데이터 불러오는 방법은 라이브러리가 로드된 상태에서 data("데이터이름")이라고 실행하면, 작업공간상에 데이터가 생성됩니다.

> head(Pima.tr)
  npreg glu bp skin  bmi   ped age type
1     5  86 68   28 30.2 0.364  24   No
2     7 195 70   33 25.1 0.163  55  Yes
3     5  77 82   41 35.8 0.156  35   No
4     0 165 76   43 47.9 0.259  26   No
5     0 107 60   25 26.4 0.133  23   No
6     5  97 76   27 35.6 0.378  52  Yes

RandomForest 모형 만들기

set.seed(10)
rf.fit = randomForest(type ~ npreg + glu + bp + skin + bmi + ped + age
                      , data=Pima.tr, mtry = floor(sqrt(7)), ntree = 500, importance = T)
rf.fit

randomForest 모형을 만드는 문법은 위와 같습니다. hyper parameter로 이 모형에서는 mtry와 ntree를 사용하였습니다. mtry는 각각의 tree마다 몇 개의 feature를 사용할 것인지를 정하는 것입니다. 골드스탠다다는 보통, regression의 경우 변수갯수/3, classification의 경우 sqrt(변수갯수)를 사용합니다. 이 데이터셋은 classification 문제이므로, sqrt를 사용하였습니다. ntree는 tree의 총 갯수를 의미합니다.

테스트 데이터 생성

test_x = Pima.te[c("npreg", "glu", "bp", "skin", "bmi", "ped", "age")]
test_y = Pima.te$type

test_y를 불러올 때는 $로 불러와야지 vector로 불러와지므로 아래 confusionMatrix가 잘 실행됩니다.

성능 평가

y_pred = predict(rf.fit, test_x)

confusionMatrix(y_pred, test_y)

Confusion Matrix and Statistics

          Reference
Prediction  No Yes
       No  192  46
       Yes  31  63
                                          
               Accuracy : 0.7681          
                 95% CI : (0.7189, 0.8124)
    No Information Rate : 0.6717          
    P-Value [Acc > NIR] : 7.678e-05       
                                          
                  Kappa : 0.455           
 Mcnemar's Test P-Value : 0.1106          
                                          
            Sensitivity : 0.8610          
            Specificity : 0.5780          
         Pos Pred Value : 0.8067          
         Neg Pred Value : 0.6702          
             Prevalence : 0.6717          
         Detection Rate : 0.5783          
   Detection Prevalence : 0.7169          
      Balanced Accuracy : 0.7195          
                                          
       'Positive' Class : No    

변수 중요도

> importance(rf.fit)
             No       Yes MeanDecreaseAccuracy MeanDecreaseGini
npreg  9.467923  1.762747             9.098049         9.308814
glu   14.543787 16.760373            20.263701        21.771324
bp     4.210273 -8.185459            -1.740227         7.920807
skin   1.370477  1.662921             2.029630         9.350372
bmi    3.431359  7.500978             7.200368        13.322145
ped    4.957824  4.666885             6.778956        13.479470
age   12.162707  6.341079            13.414806        14.528669

변수 중요도 플롯

17페이지에 random forest hyper-parameter 관련한 설명들이 나온다.

https://cran.r-project.org/web/packages/randomForest/randomForest.pdf

Random Forest 관련 튜토리얼

https://machinelearningmastery.com/tune-machine-learning-algorithms-in-r/