오늘은 내가 과제물 중 아주 고생했던 것을 소개해볼까 한다.
기준범주 로짓 모형을 계산할 수 있는 R 함수는 여러가지가 있는데, 여기서는 vglm() 함수를 사용해 보고자 한다. 일반선형회귀 모형의 경우 반응 변수가 1, 0으로 정형화 되어 있는데, 실제 생활에서는 여러 가지가 있을 수 있다.
여기서는 플로리다 악어의 길이와 주요 먹이에 대한 기준범주 로짓모형이다. 악어의 먹이(ftype)는 F(Fish, 물고기), I(Invertebrates, 연체류), O(Other, 기타)로 나뉘어 있고 길이(length)는 m 단위이다.
일단 아래와 같이 자료를 입력하는데, 여기서 문제점이 하나 존재한다. 이렇게 넣어 버리면 R은 length를 Factor로 인식해 버린다. 그렇다고 해서 transform에서 아래처럼 as.numeric()을 사용하면, 그냥 정수 형태로 전환해 버리고 만다. 오늘 이것 때문에 1시간을 날렸다. 대신 paste()를 사용하면 소수점 이하 숫자도 그대로 남겨진다.
> alg <- matrix(c(1.24, "I", 1.30, "I", 1.30, "I", 1.32, "F", 1.32, "F", 1.40, "F", 1.42, "I", 1.42, "F",
+ 1.45, "I", 1.45, "O", 1.47, "I", 1.47, "F", 1.50, "I", 1.52, "I", 1.55, "I", 1.60, "I",
+ 1.63, "I", 1.65, "O", 1.65, "I", 1.65, "F", 1.65, "F", 1.68, "F", 1.70, "I", 1.73, "O",
+ 1.78, "I", 1.78, "I", 1.78, "O", 1.80, "I", 1.80, "F", 1.85, "F", 1.88, "I", 1.93, "I",
+ 1.98, "I", 2.03, "F", 2.03, "F", 2.16, "F", 2.26, "F", 2.31, "F", 2.31, "F", 2.36, "F",
+ 2.36, "F", 2.39, "F", 2.41, "F", 2.44, "F", 2.46, "F", 2.56, "O", 2.67, "F", 2.72, "I",
+ 2.79, "F", 2.84, "F", 3.25, "O", 3.28, "O", 3.33, "F", 3.56, "F", 3.58, "F", 3.66, "F",
+ 3.68, "O", 3.71, "F", 3.89, "F"), ncol = 2, byrow = TRUE)
>
> alg <- data.frame(alg)
> colnames(alg) <- c("length", "ftype")
> head(alg)
length ftype
1 1.24 I
2 1.3 I
3 1.3 I
4 1.32 F
5 1.32 F
6 1.4 F
> ##
> alg <- transform(alg, length = as.numeric(paste(length)))
>
> head(alg)
length ftype
1 1.24 I
2 1.30 I
3 1.30 I
4 1.32 F
5 1.32 F
6 1.40 F
>
> class(alg$length)
[1] "numeric"
>
> table(alg$ftype)
F I O
31 20 8
>
vglm() 함수를 사용하기 위해서 VGAM 패키지를 먼저 설치해 놔야 한다.
> ##install.packages("VGAM")
> library("VGAM")
>
>
> head(alg)
length ftype
1 1.24 I
2 1.30 I
3 1.30 I
4 1.32 F
5 1.32 F
6 1.40 F
> table(alg$ftype)
F I O
31 20 8
> alg.model <- vglm(formula = ftype ~ length, family = multinomial, data = alg)
> summary(alg.model)
Call:
vglm(formula = ftype ~ length, family = multinomial, data = alg)
Pearson residuals:
Min 1Q Median 3Q Max
log(mu[,1]/mu[,3]) -2.330 -0.5075 0.5538 0.6836 1.452
log(mu[,2]/mu[,3]) -2.687 -0.4821 -0.1653 0.7093 3.439
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept):1 1.6177 1.3073 1.237 0.21591
(Intercept):2 5.6974 1.7937 3.176 0.00149 **
length:1 -0.1101 0.5171 -0.213 0.83137
length:2 -2.4654 0.8996 -2.741 0.00613 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Number of linear predictors: 2
Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])
Dispersion Parameter for multinomial family: 1
Residual deviance: 98.3412 on 114 degrees of freedom
Log-likelihood: -49.1706 on 114 degrees of freedom
Number of iterations: 5
Reference group is level 3 of the response
>
> coef(alg.model)
(Intercept):1 (Intercept):2 length:1 length:2
1.617731 5.697444 -0.110109 -2.465446
>
> fitted(alg.model)
F I O
1 0.2265307 0.721964000 0.05150528
2 0.2502564 0.692466814 0.05727683
3 0.2502564 0.692466814 0.05727683
내가 기준범주 로직 모형을 아무리 검색해도 제대로 값이 나오지 않아서, 이렇게 별도로 글을 남겨 놓는다. 다른 사람들도 나처럼 헤매지 말라고 그리고, 분명히 제대로 코딩을 다 했는데도 값이 이상하게 나오면 항상 원자료를 살펴보자.
잊지 말자, head() 다시 보자 head(). 원래는 코딩 안해도 되는데, 3시간 가까이 코딩을 하고 나니 과제물은 그냥 계산하고 모형의 계산식만 제대로 넣으면 되는 것이었다. 요즘 통계 대학원 수업을 들으면서 느끼는 것인데, SAS/R 이외에도 엑셀의 경우에는 과연 끝이 어디인지 아련할 정도로 강력하다.
강의 듣거나 과제물의 경우에는 일단 R로 돌려 놓고 엑셀에서 검증하는데, 생각보다 엑셀 강력하고 좋다. 엑셀이 없다면 이것을 다 손으로 적고 계산하고 해야 할텐데, 나처럼 대충대충 계산하는 사람들은 중간에 다 틀릴 것이다. 오늘도 중간에 )를 잘못 넣어서 죽다 살아 났다. 여하튼 과제 종료 20분 전에 간신히 마쳤다.
'SAS, R, 통계' 카테고리의 다른 글
| R *.gz file에서 컴파일 하는 방법 (0) | 2017.01.23 |
|---|---|
| How to generate quantile from SAS Proc Univariate (0) | 2016.10.24 |
| proc genmod moel y/n = x (0) | 2016.10.22 |
| R Transform (0) | 2016.10.21 |
| SAS Cards를 R에서 사용하기 (0) | 2016.09.28 |