Python 중고급 - 정규표현식, 클래스

정규 표현식

파이썬에서 정규표현식을 이용하는 것은 re 라이브러리를 이용하면 된다. import re 구문으로 re 라이브러리를 임포트할 수 있다. 

print(all([not re.match("a", "cat"), # cat은 a로 시작하지 않는다.
          re.search("a", "cat"), # cat에 a가 존재하기 때문에 
          not re.search("c", "dog"), # c가 dog에 존재하지 않기 때문에
          3 == len(re.split("[ab]", "carbs")), # a,b로 나누면 'c', 'r', 's' 를 원소로 하는 list가 생성된다.
          "R-D-" == re.sub("[0-9]", "-", "R2D2")])) # sub를 통해 숫자를 -로 대체

객체 지향 프로그래밍

Stack 자료구조를 파이썬의 객체 지향을 이용해서 간단하게 구현하였다. add 함수는 stack에 데이터를 저장하고, pop 함수는 stack의 가장 위에 있는 데이터를 삭제한다. stack은 파이썬의 기본 자료구조인 list를 이용하여 구현하였다. 

class stack : 
    def __init__(self) : 
        "이것은 생성자이다."
        self.list = []
        self.length = 0
        
    def add(self, value) :
        self.list.append(value)
        self.length += 1
    
    def pop(self) : 
        del self.list[self.length-1]
        
    def __repr__(self) : 
        """ 파이썬 프롬프트에서 이 함수를 입력하거나 str()으로 보내주면
            Set 객체를 문자열로 표현해줌 """
        return "Stack : " + str(self.list)

my_stack = stack()

my_stack.add(value=3)
print(my_stack)

Stack : [3]

my_stack.add(value=5)

print(my_stack)

Stack : [3, 5]

my_stack.pop()

print(my_stack)

Stack : [3]

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

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

Python “SSL: CERTIFICATE_VERIFY_FAILED” Error 해결법

Keras에서 ResNet의 weight를 다운 받다가 이러한 에러가 났다. 환경은 리눅스, 파이썬 3 버전이다.

base_model = ResNet50(weights='imagenet', pooling=max, include_top = False) 

이 코드를 실행시키다가 에러가 났는데 h5 파일을 다운 받아야하는데, SSL certificate 문제가 있다면서 다운을 받지 못하는 문제이다. 

http://blog.pengyifan.com/how-to-fix-python-ssl-certificate_verify_failed/

위 링크를 보고 해결하였는데, 

첫번째 해결법은 PYTHONHTTPSVERIFY 환경변수를 0으로 설정하는 방법이다.

두 번재 해결법은 대안적인 방법으로 문제가 있는 코드 앞에 아래 코드를 넣는 방법이다. 

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 을 활용하면 버전이 바뀌어 다시 설치되기 때문에, 패키지를 다운그레이드 시킬 때도 유용하게 쓸 수 있다.

Python - Pandas isin 구문

Python에서 테이블 형식의 데이터를 읽고 처리할 때 가장 많이 쓰이는 pandas 라이브러리에서는 다양한 데이터 처리 기능을 구현하고 있다. 이 중에 isin 구문은 열이 list의 값들을 포함하고 있는 모든 행들을 골라낼 때 주로 쓰인다. 

예를 들어, 아래 예제를 보면

df = DataFrame({'A': [1, 2, 3], 'B': ['a', 'b', 'f']})
df.isin([1, 3, 12, 'a'])

이와 같이 이진값을 반환한다. 

       A      B
0   True   True
1  False  False
2   True  False

이를 그대로 쓰는 경우보다 Dataframe의 컬럼에서 어떤 list의 값을 포함하고 있는것만 걸러낼 때 isin 구문이 유용하다.

이러한 데이터프레임이 있을 때

df = pd.DataFrame({'A': [1, 2, 3], 'B': ['a', 'b', 'f']})

A 컬럼의 값이 [1,3,12]를 포함하는 것만 골라낸다.

df[df['A'].isin([1, 3, 12])]

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'

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

SAS를 통한 표본 추출

sample.csv

임의추출

코드

FILENAME REFFILE 'C:\Users\\sample.csv';

PROC IMPORT DATAFILE=REFFILE
    DBMS=CSV
    OUT=WORK.population;
    GETNAMES=YES;
RUN;

proc surveyselect data=population method=srs n=200
    out=Work.sample;
run;

결과

표본 평균 계산

코드

proc surveymeans data=sample total=1459;
    var MSSubclass;
run;

결과

• 이 때, 표본 평균의 분산 추정량은 (N-n)/N * s^2/n 으로 계산된다. (유한 모집단이기 때문에 유한모집단 수정계수를 곱한다.)
• 모분산(sigma^2)을 아는 경우에는 (N-n)/(N-1) * sigma^2/n이다.

표본 비율의 추정

코드

data new;
    set sample;
    bin=(MSZoning='RL');
run;

proc surveymeans data=new total=1459;
    var bin;
run;

결과

층화 추출법

* LotShape라는 변수를 기준으로 층화추출한다. ;

proc sort data=population;
    by MasVnrType;
run;

proc freq data=population;
    tables MasVnrType;
run;

proc surveyselect data=population method=srs n=(5,5,5,5,5)
    out=sample2;
    strata MasVnrType;
run;

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/