1. 데이터 살펴보기

데이터가 어떤 변수로 이루어져있는지 보자.
데이터를 살펴보면 변수가 16개로 이루어져있다는 것을 알 수 있다.
모두 연속형 변수이다.

'data.frame':   252 obs. of  16 variables:
 $ brozek : num  12.6 6.9 24.6 10.9 27.8 20.6 19 12.8 5.1 12 ...
 $ age    : int  23 22 22 26 24 24 26 25 25 23 ...
 $ weight : num  154 173 154 185 184 ...
 $ height : num  67.8 72.2 66.2 72.2 71.2 ...
 $ adipos : num  23.7 23.4 24.7 24.9 25.6 26.5 26.2 23.6 24.6 25.8 ...
 $ free   : num  135 161 116 165 133 ...
 $ neck   : num  36.2 38.5 34 37.4 34.4 39 36.4 37.8 38.1 42.1 ...
 $ chest  : num  93.1 93.6 95.8 101.8 97.3 ...
 $ abdom  : num  85.2 83 87.9 86.4 100 94.4 90.7 88.5 82.5 88.6 ...
 $ hip    : num  94.5 98.7 99.2 101.2 101.9 ...
 $ thigh  : num  59 58.7 59.6 60.1 63.2 66 58.4 60 62.9 63.1 ...
 $ knee   : num  37.3 37.3 38.9 37.3 42.2 42 38.3 39.4 38.3 41.7 ...
 $ ankle  : num  21.9 23.4 24 22.8 24 25.6 22.9 23.2 23.8 25 ...
 $ biceps : num  32 30.5 28.8 32.4 32.2 35.7 31.9 30.5 35.9 35.6 ...
 $ forearm: num  27.4 28.9 25.2 29.4 27.7 30.6 27.8 29 31.1 30 ...
 $ wrist  : num  17.1 18.2 16.6 18.2 17.7 18.8 17.7 18.8 18.2 19.2 ...

> head(fat, n=3)

  brozek age weight height adipos  free neck chest abdom  hip thigh knee
1   12.6  23 154.25  67.75   23.7 134.9 36.2  93.1  85.2 94.5  59.0 37.3
2    6.9  22 173.25  72.25   23.4 161.3 38.5  93.6  83.0 98.7  58.7 37.3
3   24.6  22 154.00  66.25   24.7 116.0 34.0  95.8  87.9 99.2  59.6 38.9
  ankle biceps forearm wrist
1  21.9   32.0    27.4  17.1
2  23.4   30.5    28.9  18.2
3  24.0   28.8    25.2  16.6

결측값이 없다는 것을 확인할 수 있다.

> sum(is.na(fat)) # 결측값이 없다.

[1] 0

2. best subset approach based on adjusted R^2, Mallows Cp, BIC

성별을 변수로 추가했더라면 body fat을 더 잘 설명할 수 있었을 것이다.
하지만 현재 데이터에는 성별이 없으므로, 성별없이 분석을 진행해보도록 한다.
새로이 데이터를 수집할 기회가 있다면 성별을 포함하는 것이 더 좋을 것 같다.

2.1선형회귀분석을 돌려보자.

아래는 변수개수(p)에 따라 구한 최적의 모형들이다.
예를 들어, 변수를 한개만 포함시킬 때는(p=1) 변수 abdom을 사용하는 것이 좋으며, 변수를 두개 포함시킬때는(p=2) 변수 free와 weight을 포함시켜야 모형이 좋다고 한다.

          age weight height adipos free neck chest abdom hip thigh knee
1  ( 1 )  " " " "    " "    " "    " "  " "  " "   "*"   " " " "   " " 
2  ( 1 )  " " "*"    " "    " "    "*"  " "  " "   " "   " " " "   " " 
3  ( 1 )  " " "*"    " "    " "    "*"  " "  " "   " "   " " " "   " " 
4  ( 1 )  " " "*"    " "    " "    "*"  " "  " "   "*"   " " " "   " " 
5  ( 1 )  " " "*"    " "    "*"    "*"  " "  " "   "*"   " " " "   " " 
6  ( 1 )  " " "*"    " "    "*"    "*"  " "  " "   "*"   " " "*"   " " 
7  ( 1 )  " " "*"    " "    "*"    "*"  " "  "*"   "*"   " " "*"   " " 
8  ( 1 )  " " "*"    " "    "*"    "*"  " "  "*"   "*"   " " "*"   " " 
9  ( 1 )  " " "*"    " "    "*"    "*"  " "  "*"   "*"   " " "*"   " " 
10  ( 1 ) " " "*"    " "    "*"    "*"  " "  "*"   "*"   " " "*"   " " 
11  ( 1 ) " " "*"    "*"    "*"    "*"  " "  "*"   "*"   " " "*"   " " 
12  ( 1 ) " " "*"    "*"    "*"    "*"  " "  "*"   "*"   " " "*"   "*" 
13  ( 1 ) "*" "*"    "*"    "*"    "*"  " "  "*"   "*"   " " "*"   "*" 
14  ( 1 ) "*" "*"    "*"    "*"    "*"  "*"  "*"   "*"   " " "*"   "*" 
15  ( 1 ) "*" "*"    "*"    "*"    "*"  "*"  "*"   "*"   "*" "*"   "*" 
          ankle biceps forearm wrist
1  ( 1 )  " "   " "    " "     " "  
2  ( 1 )  " "   " "    " "     " "  
3  ( 1 )  " "   " "    "*"     " "  
4  ( 1 )  " "   " "    "*"     " "  
5  ( 1 )  " "   " "    "*"     " "  
6  ( 1 )  " "   " "    "*"     " "  
7  ( 1 )  " "   " "    "*"     " "  
8  ( 1 )  "*"   " "    "*"     " "  
9  ( 1 )  "*"   "*"    "*"     " "  
10  ( 1 ) "*"   "*"    "*"     "*"  
11  ( 1 ) "*"   "*"    "*"     "*"  
12  ( 1 ) "*"   "*"    "*"     "*"  
13  ( 1 ) "*"   "*"    "*"     "*"  
14  ( 1 ) "*"   "*"    "*"     "*"  
15  ( 1 ) "*"   "*"    "*"     "*"

2.2 몇개 변수를 추가하는 것 좋을지 살펴보자

p에 따라 \[R^2\], \(R_{adj}^2\), Mallow’s Cp와 BIC 통계량이 어떻게 변화하는지 살펴보자.
각각 변수 15개, 10개, 12개,7개 추가할 때 최적의 모형을 얻을 수 있다.
하지만 개인적으로 판단했을 때, 변수 2개 추가했을 때가 가장 적절해보인다. 변수 3개~15개 추가했을 때는 통계량 값은 크게 차이가 나지 않기 때문이다.

2.3 위 결과를 이용해서 모수를 추정하여 최종모형식을 만들어보자.

변수를 2개 사용한다고 했을 때, 모형식을 도출해보자. \(Y= 19.53 + 0.422Weight - 0.53free + \epsilon\)

> coef(regfit.full,2)

(Intercept)      weight        free 
 19.6535898   0.4229089  -0.5314993

2.4 회귀모형을 진단해보자.

QQplot의 꼬리 부분이 45도 각도 직성에서 많이 벗어나있다. 그래서 정규성이 위배될 수도 있을 것으로 보인다.
잔차&예측치를 확인해보았더니, 체계적인 관계가 있는 것으로 보인다. 종속변수와 독립변수 사이에 선형적이기 위해서는 잔차와 예측치 사이에 체계적인 관계가 있으면 안된다. 다항식으로 바꾸면 문제를 해결할 수 있을 것으로 보인다.

직교다항식으로 바꾸었을 때 회귀진단결과가 더 좋아보인다.
일단 잔차&예측치plot을 보면 체계적인 관계가 안 보인다. 종속변수와 독립변수사이의 선형적인 관계가 만족되었다고 볼 수 있겠다.
QQplot에 찍힌 점들도 45도 각도 직선으로부터 많이 벗어나지 않는 것으로 보아 정규성도 만족된다고 판단한다.
따라서 분석을 진행하게 된다면 polynomial하게 바꾸는 것이 좋겠다.

3 validation set method

변수 2개를 포함하는 것이 가장 안정적으로 보인다.
변수를 3개 포함하는 순간 값이 불안해진다.
따라서 최적의 모형은 변수 2개를 넣었을 때이고, 그때의 회귀계수는 아래와 같다. \(Brozek= 18.5 + 0.43Weight - 0.53free\)

> coef(regfit.best,2)

(Intercept)      weight        free 
 18.0518456   0.4321000  -0.5308647

Validation Set approach로 수행한 결과도 회귀진단을 거쳐야하지만 생략하도록 한다.

4. LOOCV

변수를 2개 포함시킬 때 MSE가 가장 작다.
따라서 최적의 모형은 변수 2개를 넣었을 때이고, 그때의 회귀계수는 아래와 같다. \(Brozek= 18.5 + 0.43Weight - 0.53free\)

> coef(regfit.best,2)

(Intercept)      weight        free 
 18.0518456   0.4321000  -0.5308647

LOOCV로 수행한 결과도 회귀진단을 거쳐야하지만 생략하도록 한다.

5.10-fold CV —————————–

변수가 2개일 때 MSE값이 가장 작다.
따라서 최적의 모형은 변수 2개를 넣었을 때이고, 그때의 회귀계수는 아래와 같다. \(Brozek= 18.5 + 0.43Weight - 0.53free\)

> coef(regfit.best,2)

(Intercept)      weight        free 
 18.0518456   0.4321000  -0.5308647

-k-fold로 수행한 결과도 회귀진단을 거쳐야하지만 생략하도록 한다.

6. Appendix_Code

> # 데이터 살펴보기
> Sys.setlocale("LC_ALL", "eng")
> setwd("E:/Dropbox/00.2018/01.2018_1_semester/01.DataMining/04.HW/HW2")
> fat <- read.csv("./fat.csv", header = T)
> str(fat)
> head(fat, n=3)
> sum(is.na(fat)) # 결측값이 없다. 
> 
> # Best Subset Approach
> library(leaps)
> 
> regfit.full=regsubsets(brozek~.,fat)
> regfit.full=regsubsets(brozek~.,data=fat,nvmax=15)
> reg.summary=summary(regfit.full)
> reg.summary$outmat
> 
> # 결과를 보자. 몇 개의 변수를 추가하는 것이 좋은가. 
> par(mfrow=c(2,2))
> 
> plot(reg.summary$rss, xlab="Number of Variables",
+ ylab = "RSS",
+ type="l")
> points(x = which(reg.summary$rss==min(reg.summary$rss)),y =  min(reg.summary$rss), col="red", pch=20, cex=2)
> 
> plot(reg.summary$cp, xlab="Number of Variables",
+ ylab = "Mellow's Cp",
+ type="l")
> points(x = which(reg.summary$cp==min(reg.summary$cp)),y =  min(reg.summary$cp), col="red", pch=20, cex=2)
> 
> plot(reg.summary$adjr2, xlab="Number of Variables",
+ ylab = "Adjusted R-square",
+ type="l")
> points(x = which(reg.summary$adjr2==max(reg.summary$adjr2)),y =  max(reg.summary$adjr2), col="red", pch=20, cex=2)
> 
> plot(reg.summary$bic, xlab="Number of Variables",
+ ylab = "BIC",
+ type="l")
> points(x = which(reg.summary$bic==min(reg.summary$bic)),y =  min(reg.summary$bic), col="red", pch=20, cex=2)
> 
> # 최종모형식 
> coef(regfit.full,2) 
> 
> # 회귀진단
> par(mfrow=c(2,2))
> model <- lm(brozek~ weight+ free, data=fat)
> plot(model)
> 
> # VS approach
> set.seed(20141630)
> train=sample(c(TRUE,FALSE), nrow(fat), replace = TRUE)
> test=(!train)
> 
> library(leaps)
> regfit.best <- regsubsets(brozek ~ ., data = fat[train, ], nvmax = 15)
> test.mat <- model.matrix(brozek~., data = fat[test, ])
> val.error <- rep(NA, 15)
> 
> 
> for (j in 1:20) {
+ set.seed(20141630 * j)
+ train <- sample(c(TRUE, FALSE), nrow(fat), replace = TRUE)
+ test <- (!train)
+ regfit.best <- regsubsets(brozek ~ ., data = fat[train, ], nvmax = 15)
+ test.mat <- model.matrix(brozek~., data = fat[test, ])
+ val.error <- rep(NA, 15)
+ for (i in 1:15) {
+ coefi <- coef(regfit.best, id = i)
+ pred <- test.mat[, names(coefi)] %*% coefi
+ val.error[i] <- mean((fat$brozek[test] - pred) ^ 2)
+ }
+ if (j < 2) {
+ plot(
+ val.error, type = "l", ylim = c(0, max(val.error) * 1.5),
+ ylab = "MSE", xlab = "number of variables(p)"
+ )
+ } else {
+ lines(val.error, type = "l", col = j)
+ }
+ }
> 
> 
> predict.regsubsets <- function(object, newdata, id, ...) {
+ form <- as.formula(object$call [[2]])
+ mat <- model.matrix(form, newdata)
+ coefi <- coef(object, id = id)
+ xvars <- names(coefi)
+ mat[, xvars] %*% coefi
+ }
> 
> 
> # LOOCV
> k=nrow(fat)
> cv.errors <- matrix(NA, k, 15, dimnames=list(NULL, paste(1:15)))
> 
> library(leaps)
> for(j in 1:k){
+ best.fit <- regsubsets(brozek~. ,data=fat[-j,], nvmax=15)
+ for(i in 1:15){
+ pred=predict(best.fit, fat[j,],id=i)
+ cv.errors[j,i]=mean((fat$brozek[j]- pred)^2)
+ }
+ }
> 
> mean.cv.errors <- apply(cv.errors, 2, mean)
> tmp <- apply(cv.errors, 1, mean)
> 
> par(mfrow=c(1,1))
> plot(mean.cv.errors, type="l")
> points(x = which(mean.cv.errors==min(mean.cv.errors)),y =  min(mean.cv.errors), col="red", pch=20, cex=2)
> 
> 
> ## K-fold
> k <- 10
> # k=5
> library(RColorBrewer)
> library(leaps)
> set.seed(20141630)
> for (j in 1:k) {
+ folds <- sample(1:k, nrow(fat), replace = TRUE, prob = rep(1 / k, k))
+ cv.errors <- rep(NA, 15)
+ for (i in 1:15) {
+ pred <- predict(best.fit, fat[folds == k, ], id = i)
+ cv.errors[i] <- mean((fat$brozek[folds == k] - pred) ^ 2)
+ }
+ if (j < 2) {
+ plot(
+ cv.errors, type = "l", ylim = c(0, max(cv.errors) * 2),
+ ylab = "MSE", xlab = "number of variables(p)"
+ )
+ } else {
+ lines(cv.errors, type = "l", col = brewer.pal(9, "Set1")[j])
+ 1}
+ }

Data mining HW2

1. 데이터 살펴보기

2. best subset approach based on adjusted R^2, Mallows Cp, BIC

2.1선형회귀분석을 돌려보자.

2.2 몇개 변수를 추가하는 것 좋을지 살펴보자

2.3 위 결과를 이용해서 모수를 추정하여 최종모형식을 만들어보자.

2.4 회귀모형을 진단해보자.

3 validation set method

4. LOOCV

5.10-fold CV —————————–

6. Appendix_Code