[함수의 구조] lm(formula, data, na.action)
[기능] 선형 모형을 적합한다.
[option][함수의 구조] glm(formula, family, data, na.action)
[기능] glm은 일반화선형모형을 적합하기 위해 사용
[option][함수의 구조] step(object, direction = c(“both”, “backward”, “forward”))
[기능] step() 함수는 AIC를 이용한 모형 선택을 위해 사용한다.
[option]####(4) 예측함수
[structure] predict(object, newdata, type)
[function] 다양한 모형 적합 결과로부터 예측값을 계산할 때 사용한다.
[option]
# importing data
prod <- read.csv("./data/productivityREG.csv", header = T)
head(prod)
범주형 변수 변환
# factoring predictor variables
prod$quarter <- factor(prod$quarter)
prod$department <- factor(prod$department)
prod$day <- factor(prod$day)
prod$team <- factor(prod$team)
회귀모형 적합
# fitting a linear regression model
fit.all <- lm(productivity ~ ., data = prod)“direction =” : 변수 선택 방법, ’both’는 단계적 변수 선택, ’backward’는 후진소거, ’forward’는 전진 선택법
# selecting model
fit.step <- step(fit.all, direction = "both")## Start: AIC=-3977.14
## productivity ~ quarter + department + day + team + target + smv +
## wip + over_time + incentive + numworkers
##
## Df Sum of Sq RSS AIC
## - day 5 0.10729 12.695 -3979.2
## - wip 1 0.02050 12.608 -3977.6
## - over_time 1 0.02508 12.613 -3977.3
## <none> 12.588 -3977.1
## - quarter 4 0.23909 12.827 -3967.5
## - numworkers 1 0.17880 12.766 -3965.9
## - smv 1 0.27299 12.861 -3959.1
## - department 1 0.36581 12.953 -3952.3
## - team 11 0.76326 13.351 -3944.0
## - target 1 0.73169 13.319 -3926.3
## - incentive 1 1.36454 13.952 -3882.8
##
## Step: AIC=-3979.2
## productivity ~ quarter + department + team + target + smv + wip +
## over_time + incentive + numworkers
##
## Df Sum of Sq RSS AIC
## - wip 1 0.02537 12.720 -3979.3
## <none> 12.695 -3979.2
## - over_time 1 0.03878 12.734 -3978.3
## + day 5 0.10729 12.588 -3977.1
## - quarter 4 0.21936 12.914 -3971.2
## - numworkers 1 0.18567 12.880 -3967.6
## - smv 1 0.25487 12.950 -3962.6
## - department 1 0.38288 13.078 -3953.4
## - team 11 0.75550 13.450 -3947.1
## - target 1 0.72922 13.424 -3928.9
## - incentive 1 1.35667 14.052 -3886.2
##
## Step: AIC=-3979.33
## productivity ~ quarter + department + team + target + smv + over_time +
## incentive + numworkers
##
## Df Sum of Sq RSS AIC
## <none> 12.720 -3979.3
## + wip 1 0.02537 12.695 -3979.2
## - over_time 1 0.03445 12.755 -3978.8
## + day 5 0.11216 12.608 -3977.6
## - quarter 4 0.21427 12.934 -3971.7
## - numworkers 1 0.19492 12.915 -3967.1
## - smv 1 0.28616 13.006 -3960.5
## - department 1 0.37045 13.091 -3954.5
## - team 11 0.76557 13.486 -3946.6
## - target 1 0.70411 13.424 -3930.9
## - incentive 1 1.64287 14.363 -3867.6
위의 변수 선택에 의해서 어떤 변수가 제외되거나 선택되었는지 확인
# names(fit.step)
fit.step$anova입력변수 10개 중에서 ‘day’, ’wip’을 제외한 8개를 가진 모형의 AIC가 가장 작아서 최종 선택 모형이 되었다.
summary(fit.step)##
## Call:
## lm(formula = productivity ~ quarter + department + team + target +
## smv + over_time + incentive + numworkers, data = prod)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.42738 -0.04496 0.00437 0.06740 0.33788
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.682e-01 6.287e-02 5.856 6.62e-09 ***
## quarterQuarter2 -1.400e-05 1.005e-02 -0.001 0.998889
## quarterQuarter3 -1.080e-02 1.177e-02 -0.918 0.359096
## quarterQuarter4 -2.310e-02 1.167e-02 -1.978 0.048176 *
## quarterQuarter5 5.591e-02 1.993e-02 2.805 0.005140 **
## departmentsweing -1.637e-01 3.173e-02 -5.159 3.04e-07 ***
## team2 -5.849e-03 1.777e-02 -0.329 0.742076
## team3 1.015e-02 1.893e-02 0.536 0.592040
## team4 -4.137e-03 1.839e-02 -0.225 0.822105
## team5 4.814e-03 1.960e-02 0.246 0.805988
## team6 -4.602e-02 2.134e-02 -2.157 0.031271 *
## team7 -5.645e-02 1.954e-02 -2.888 0.003968 **
## team8 -6.962e-02 1.845e-02 -3.775 0.000171 ***
## team9 -5.156e-02 1.770e-02 -2.913 0.003666 **
## team10 -3.348e-02 1.790e-02 -1.871 0.061701 .
## team11 -8.259e-02 2.026e-02 -4.076 4.97e-05 ***
## team12 -1.197e-02 2.007e-02 -0.596 0.551090
## target 5.745e-01 8.077e-02 7.113 2.29e-12 ***
## smv -4.527e-03 9.983e-04 -4.535 6.54e-06 ***
## over_time -2.973e-06 1.889e-06 -1.573 0.115997
## incentive 2.596e-03 2.389e-04 10.865 < 2e-16 ***
## numworkers 2.920e-03 7.802e-04 3.742 0.000194 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.118 on 914 degrees of freedom
## Multiple R-squared: 0.3333, Adjusted R-squared: 0.3179
## F-statistic: 21.76 on 21 and 914 DF, p-value: < 2.2e-16예를 들어, 입력변수 incentive의 회귀계수 추정값은 양수이므로 incentive가 증가함에 따라 목표변수 productivity가 증가한다. 범주형 변수인 quarter는 5개의 범주로 구성되어 있기 때문에 4개의 가변수가 생성되었다. 각 입력변수의 P-value는 0.05보다 작은 경우에 유의하다고 할 수 있는데, 이의 가변수 중 하나 이상이 유의하여 다른 가변수가 유의하므로 제거되지 않고 모형에 포함된다. R-squared는 33.33% (Adjusted R-squared는 31.79%)로 적합한 선형회귀모형으로 데이터를 설명할 수 있는 부분이 약 33%이고, F-검정의 p-value가 아주 작아서 모형이 적합하다고 할 수 있다.
2-3에서는 최종모형으로 목표변수의 예측값을 계산하여 예측값을 산출한다. 이로부터 생성된 실제값과 예측값 차이의 제곱합 평균 MSE 0.01358993과 절대치 평균 MAE 0.08336836을 예측력의 평가측도로 사용할 수 있다. 이 값은 작을 수록 더 예측력이 높다고 할 수 있다.
# making predictions
pred.reg <- predict(fit.step, newdata = prod, type = "response")
# pred.reg# evaluation
mean((prod$productivity - pred.reg)^2) # MSE## [1] 0.01358993mean(abs(prod$productivity - pred.reg)) # MAE## [1] 0.08336836먼저 와인 품질데이터를 호출한다. 이 데이터에서 이항형 변수 quality는 두 개의 반응값 1과 0을 가지고 나머지 11개는 입력 변수로 수치형이다. 이항형 목표변수의 경우에 로지스틱회귀모형을 적합하므로 glm 함수에서 모형 공식 ’quality ~ .’를 사용하여 모든 입력인자를 모형에 적합시킨다. 이때 분포는 이항으로 지정 (family = binomial)한다. AIC가 가장 작은 최종 모형은 step 함수를 이용하여 수행한 단계적 선택법 (direction = “both”)에 의해 유의한 입력 변수들이 선택된다. 단계적으로 density, residsugar, fixed 입력변수가 제외되고 나머지 8개만 남게 된다.
# importing data
wine <- read.csv("./data/winequalityCLASS.csv", header = T)로지스틱 회귀 모형 적합
# fitting a logistic regression model
fit.all <- glm(quality ~ ., family = binomial, data = wine)“family = binomial” : 분포는 이항분포로 지정
# AIC --> selecting model
fit.step <- step(fit.all, direction = "both") # stepwise variable selection## Start: AIC=1249.78
## quality ~ fixed + volatile + citric + residsugar + chlorides +
## freeSD + totalSD + density + pH + sulphates + alcohol
##
## Df Deviance AIC
## - density 1 1225.8 1247.8
## - residsugar 1 1226.0 1248.0
## - fixed 1 1226.2 1248.2
## <none> 1225.8 1249.8
## - chlorides 1 1228.2 1250.2
## - pH 1 1228.6 1250.6
## - freeSD 1 1230.5 1252.5
## - citric 1 1232.4 1254.4
## - totalSD 1 1237.1 1259.1
## - volatile 1 1247.3 1269.3
## - alcohol 1 1275.1 1297.1
## - sulphates 1 1283.4 1305.4
##
## Step: AIC=1247.84
## quality ~ fixed + volatile + citric + residsugar + chlorides +
## freeSD + totalSD + pH + sulphates + alcohol
##
## Df Deviance AIC
## - residsugar 1 1226.0 1246.0
## - fixed 1 1227.4 1247.4
## <none> 1225.8 1247.8
## - chlorides 1 1228.2 1248.2
## - pH 1 1229.4 1249.4
## + density 1 1225.8 1249.8
## - freeSD 1 1230.5 1250.5
## - citric 1 1232.4 1252.4
## - totalSD 1 1237.2 1257.2
## - volatile 1 1247.5 1267.5
## - sulphates 1 1288.6 1308.6
## - alcohol 1 1347.2 1367.2
##
## Step: AIC=1246.05
## quality ~ fixed + volatile + citric + chlorides + freeSD + totalSD +
## pH + sulphates + alcohol
##
## Df Deviance AIC
## - fixed 1 1227.5 1245.5
## <none> 1226.0 1246.0
## - chlorides 1 1228.8 1246.8
## + residsugar 1 1225.8 1247.8
## - pH 1 1229.9 1247.9
## + density 1 1226.0 1248.0
## - freeSD 1 1230.8 1248.8
## - citric 1 1232.6 1250.6
## - totalSD 1 1238.7 1256.7
## - volatile 1 1248.1 1266.1
## - sulphates 1 1289.2 1307.2
## - alcohol 1 1351.7 1369.7
##
## Step: AIC=1245.45
## quality ~ volatile + citric + chlorides + freeSD + totalSD +
## pH + sulphates + alcohol
##
## Df Deviance AIC
## <none> 1227.5 1245.5
## - chlorides 1 1229.8 1245.8
## + fixed 1 1226.0 1246.0
## + density 1 1226.8 1246.8
## + residsugar 1 1227.4 1247.4
## - freeSD 1 1232.5 1248.5
## - citric 1 1232.7 1248.7
## - pH 1 1238.3 1254.3
## - totalSD 1 1241.8 1257.8
## - volatile 1 1248.2 1264.2
## - sulphates 1 1295.2 1311.2
## - alcohol 1 1351.9 1367.9fit.step$anovadirection = ‘both’, 단계적 선택법에 의해 유의미한 입력변수 선택/제외
아래의 2-5에서 보듯이 최종모형은 8개의 입력변수로 구성되어 있고 회귀계수의 부호와 크기로 해석한다. 예를 들어 휘발산(volatile)의 회귀계수가 음수이므로 이것의 수치가 높을 수록 와인의 품질(quality)이 낮을 가능성이 높다. volatile 수치 한 단위가 증가할 수록 오즈(odds)가 exp(-2.625798)=0.027238179배로 줄어든다. Null deviance는 상수항만 가진 모형의 이탈도, Residual deviance는 최종모형의 이탈도를 의미하고, 최종모형의 AIC는 1245.5이다.
summary(fit.step)##
## Call:
## glm(formula = quality ~ volatile + citric + chlorides + freeSD +
## totalSD + pH + sulphates + alcohol, family = binomial, data = wine)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8269 -0.8366 0.2804 0.8323 2.2659
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.598450 2.246265 -1.602 0.109162
## volatile -2.625798 0.587804 -4.467 7.93e-06 ***
## citric -1.304024 0.570926 -2.284 0.022369 *
## chlorides -8.134165 5.330701 -1.526 0.127033
## freeSD 0.024636 0.010981 2.244 0.024862 *
## totalSD -0.014548 0.003919 -3.712 0.000205 ***
## pH -2.029157 0.620285 -3.271 0.001070 **
## sulphates 5.609315 0.722816 7.760 8.47e-15 ***
## alcohol 0.933382 0.091576 10.192 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1647.5 on 1193 degrees of freedom
## Residual deviance: 1227.5 on 1185 degrees of freedom
## AIC: 1245.5
##
## Number of Fisher Scoring iterations: 4아래 2-6에서 보듯 최종모형이 저장된 fit.step을 predict 함수에 적용하여 목표변수의 예측값을 계산하고 임계치(cutoff) 0.5보다 크면 품질(quality)이 우수(=1)로, 아니면 품질이 보통(=0)으로 예측한다.
적합한 모형에서 얻은 예측값과 관측값으로 만든 정오분류표(confusion matrix)는 2-6과 같이 만들 수 있다. 실제 0인데 0으로 예측한 경우가 402개, 1인데 1로 분류한 경우가 495개다. 반면 0에서 1로, 1에서 0으로 오분류한 경우가 각각 147개, 150개이다. 예측 정확도는 약 75.1%, 즉 오분류률은 24.9%이다. 민감도와 특이도는 각각 76.7%, 73.2%이다. 현재까지는 1로 예측할 확률이 임계치 0.5보다 클 경우에 1로, 0.5 이하일 경우에 0으로 예측하였다. 다른 임계치를 선택하기 위해서는 cutoff를 0.5가 아닌 다른 값으로 지정하면 된다.
# making predictions
p <- predict(fit.step, newdata = wine, type = "response")
cutoff <- 0.5 # cutoff
yhat <- ifelse(p > cutoff, 1, 0)# evaluation
tab <- table(wine$quality, yhat, dnn = c("Observed", "Predicted"))
print(tab)## Predicted
## Observed 0 1
## 0 402 147
## 1 150 495sum(diag(tab)) / sum(tab) # accuracy## [1] 0.7512563tab[2, 2] / sum(tab[2, ]) # sensitivity## [1] 0.7674419tab[1, 1] / sum(tab[1, ]) # specificity## [1] 0.7322404