statistical Data Analysis 과제 4
chaemj
2024-04-29
Default data
- In Chapter 4, we used logistic regression to predict the probability of default using income and balance on the Default data set. We will now estimate the test error of this logistic regression model using the validation set approach. Do not forget to set a random seed before beginning your analysis.
library(ISLR2)## Warning: 패키지 'ISLR2'는 R 버전 4.4.0에서 작성되었습니다library(MASS)##
## 다음의 패키지를 부착합니다: 'MASS'## The following object is masked from 'package:ISLR2':
##
## Bostonhead(Default)## default student balance income
## 1 No No 729.5265 44361.625
## 2 No Yes 817.1804 12106.135
## 3 No No 1073.5492 31767.139
## 4 No No 529.2506 35704.494
## 5 No No 785.6559 38463.496
## 6 No Yes 919.5885 7491.559##변수 설명
default: A factor with levels NO and Yes indicating whether the customer
defaulted on their debt
student: A factor with levels No and Yes indicating whether the customer
is a student
balance: The average balance that the customer has remaining on their
credit card after making their monthly payment
income: income of customer
#기초 통계량
summary(Default)## default student balance income
## No :9667 No :7056 Min. : 0.0 Min. : 772
## Yes: 333 Yes:2944 1st Qu.: 481.7 1st Qu.:21340
## Median : 823.6 Median :34553
## Mean : 835.4 Mean :33517
## 3rd Qu.:1166.3 3rd Qu.:43808
## Max. :2654.3 Max. :73554default와 student는 범주형 데이터이고 balance와 income은 연속형 데이터이다. default를 종속변수로 하고 나머지 데이터들을 독립변수로 하는 분석을 진행한다.
#boxplots
par(mfrow=c(1,3))
boxplot(balance ~ default, data = Default, main = "balance & default")
boxplot(income ~ default, data = Default, main = "income & default")
범주형 데이터인 default와 연속형 데이터인 balance와 income의 boxplot을
그려보면, default 변수에 따라 balance의 차이는 크고 income의 차이는 크지
않다는 것을 알 수 있다.
이는 채무불이행인 그룹과 아닌 그룹 간의 소득 차이는 크지 않지만,
신용카드에 남은 빚진 금액의 차이는 크다는 것을 알 수있다.
#t-test
t.test(Default$balance~Default$default)##
## Welch Two Sample t-test
##
## data: Default$balance by Default$default
## t = -48.984, df = 374.14, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group No and group Yes is not equal to 0
## 95 percent confidence interval:
## -981.7670 -905.9889
## sample estimates:
## mean in group No mean in group Yes
## 803.9438 1747.8217t.test(Default$income~Default$default)##
## Welch Two Sample t-test
##
## data: Default$income by Default$default
## t = 1.922, df = 353.62, p-value = 0.05541
## alternative hypothesis: true difference in means between group No and group Yes is not equal to 0
## 95 percent confidence interval:
## -34.38335 2988.42235
## sample estimates:
## mean in group No mean in group Yes
## 33566.17 32089.15t-test 결과:
balance와 default의 t-test의 p-value의 값은 2.2e-16으로 아주 작아서 두
그룹 간의 차이가 없다는 귀무가설을 거절할 수 있다.
income과 default의 t-test의 p-value의 값은 0.05541으로 유의수준 0.1에서
귀무가설을 기각할 수 있다.
#student&default 분할표와 카이제곱 검정
std_def=xtabs(~ student+default, data=Default)
std_def## default
## student No Yes
## No 6850 206
## Yes 2817 127chisq.test(std_def)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: std_def
## X-squared = 12.117, df = 1, p-value = 0.0004997범위 변수인 student와 default의 분할표를 만들고 카이제곱 검정을 해보면, p-value값이 0.0004997으로 귀무가설을 기각하여 채무 불이행 여부와 학생 여부가 서로 연관이 있다고 결론내릴 수 있다.
- Fit a logistic regression model that uses income and balance to predict default.
dlogis<-glm(default~income+balance, data= Default, family=binomial)
summary(dlogis)##
## Call:
## glm(formula = default ~ income + balance, family = binomial,
## data = Default)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.154e+01 4.348e-01 -26.545 < 2e-16 ***
## income 2.081e-05 4.985e-06 4.174 2.99e-05 ***
## balance 5.647e-03 2.274e-04 24.836 < 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: 2920.6 on 9999 degrees of freedom
## Residual deviance: 1579.0 on 9997 degrees of freedom
## AIC: 1585
##
## Number of Fisher Scoring iterations: 8income과 balance의 추정된 계수의 값이 모두 높게 나왔다. 추정된 값 모두 양의 값이 나왔으므로 각 독립 변수가 증가할수록 채무불이행할 확률이 높다는 것을 알 수있다.
- Using the validation set approach, estimate the test error of this model. In order to do this, you must perform the following steps:
- Split the sample set into a training set and validation set.
set.seed(1)
train_d<-sample(dim(Default)[1], dim(Default)[1]/2)
test_d<-Default[-train_d,]전체 데이터를 반으로 나누어서 반은 학습용으로 나머지는 테스트용도로 사용한다.
- Fit a multiple logistic regression model using only the training observations.
dlogis_tr<-glm(default~income+balance, data=Default, family=binomial, subset=train_d)
summary(dlogis_tr)##
## Call:
## glm(formula = default ~ income + balance, family = binomial,
## data = Default, subset = train_d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.194e+01 6.178e-01 -19.333 < 2e-16 ***
## income 3.262e-05 7.024e-06 4.644 3.41e-06 ***
## balance 5.689e-03 3.158e-04 18.014 < 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: 1523.8 on 4999 degrees of freedom
## Residual deviance: 803.3 on 4997 degrees of freedom
## AIC: 809.3
##
## Number of Fisher Scoring iterations: 8추정값들의 절댓값이 커진 것외에는 앞의 결과와 비슷한 결과가 나왔다.
- Obtain a prediction of default status for each individual in the validation set by computing the posterior provability of default for that individual, adn classifying the individual to the default category if the posterior probability is greater than 0.5.
dlogis.pred=rep("No", dim(Default)[1]/2)
dlogis.probs=predict(dlogis, test_d, type="response")
dlogis.pred[dlogis.probs>0.5]="Yes"
table(dlogis.pred, test_d$default)##
## dlogis.pred No Yes
## No 4826 107
## Yes 17 50표를 확인해보면 정답을 맞춘 경우가 그렇지 않은 경우보다 훨씬 더 많다는 것을 알 수 있다.
- Compute the validation set error, which is the fraction of the observations in the validation set that are misclassified.
#error rate
mean(dlogis.pred!=Default$default)## [1] 0.0455error rate은 4.55%로 낮게 나왔다.
- Repeat the process in (b) three times, using the different splits of the observations into a training set and a validation set. Coment on the results obtained.
train_d1<-sample(dim(Default)[1], dim(Default)[1]*0.6)
test_d1<-Default[-train_d1,]
dlogis_tr1<-glm(default~income+balance, data=Default, family=binomial, subset=train_d1)
summary(dlogis_tr1)##
## Call:
## glm(formula = default ~ income + balance, family = binomial,
## data = Default, subset = train_d1)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.192e+01 5.796e-01 -20.570 < 2e-16 ***
## income 1.989e-05 6.353e-06 3.131 0.00174 **
## balance 5.902e-03 3.046e-04 19.378 < 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: 1787.28 on 5999 degrees of freedom
## Residual deviance: 923.41 on 5997 degrees of freedom
## AIC: 929.41
##
## Number of Fisher Scoring iterations: 8dlogis1.pred=rep("No", dim(test_d1)[1])
dlogis1.probs=predict(dlogis_tr1, test_d1, type="response")
dlogis1.pred[dlogis1.probs>0.5]="Yes"
table(dlogis1.pred, test_d1$default)##
## dlogis1.pred No Yes
## No 3850 92
## Yes 22 36mean(dlogis1.pred!=test_d1$default)## [1] 0.0285위의 모델보다 income의 계수 값이 줄어들고 pvalue가 조금 커졌다. error rate은 2.85%으로 줄어들었다.
train_d2<-sample(dim(Default)[1], dim(Default)[1]*0.7)
test_d2<-Default[-train_d2,]
dlogis_tr2<-glm(default~income+balance, data=Default, family=binomial, subset=train_d2)
summary(dlogis_tr2)##
## Call:
## glm(formula = default ~ income + balance, family = binomial,
## data = Default, subset = train_d2)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.151e+01 5.179e-01 -22.218 < 2e-16 ***
## income 1.892e-05 5.898e-06 3.207 0.00134 **
## balance 5.671e-03 2.708e-04 20.940 < 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: 2070.7 on 6999 degrees of freedom
## Residual deviance: 1105.9 on 6997 degrees of freedom
## AIC: 1111.9
##
## Number of Fisher Scoring iterations: 8dlogis2.pred=rep("No", dim(test_d2)[1])
dlogis2.probs=predict(dlogis_tr2, test_d2, type="response")
dlogis2.pred[dlogis2.probs>0.5]="Yes"
table(dlogis2.pred, test_d2$default)##
## dlogis2.pred No Yes
## No 2898 64
## Yes 6 32mean(dlogis2.pred!=test_d2$default)## [1] 0.02333333error rate이 2.33%으로 줄어들었다.
train_d3<-sample(dim(Default)[1], dim(Default)[1]*0.8)
test_d3<-Default[-train_d3,]
dlogis_tr3<-glm(default~income+balance, data=Default, family=binomial, subset=train_d3)
summary(dlogis_tr3)##
## Call:
## glm(formula = default ~ income + balance, family = binomial,
## data = Default, subset = train_d3)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.193e+01 5.149e-01 -23.168 < 2e-16 ***
## income 2.152e-05 5.718e-06 3.764 0.000167 ***
## balance 5.848e-03 2.694e-04 21.705 < 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: 2266.0 on 7999 degrees of freedom
## Residual deviance: 1194.9 on 7997 degrees of freedom
## AIC: 1200.9
##
## Number of Fisher Scoring iterations: 8dlogis3.pred=rep("No", dim(test_d3)[1])
dlogis3.probs=predict(dlogis_tr3, test_d3, type="response")
dlogis3.pred[dlogis3.probs>0.5]="Yes"
table(dlogis3.pred, test_d3$default)##
## dlogis3.pred No Yes
## No 1911 53
## Yes 12 24mean(dlogis3.pred!=test_d3$default)## [1] 0.0325error rate이 3.25%로 증가했다. 위의 결과들을 종합해보면 training set 70%, test set 30%로 하는 것이 가장 적절해보인다.
- Now consider a logistic regression model that predicts the probability of default using income, balance, and a dummy variable for student. Estimate the test error for this model using the validation set approach. Comment on whether or not including a dummy variable for student leads to a reduction in the test error rate.
error rate이 가장 낮았던 모델(training set 70%, test set 30%)에 student 변수를 추가해보았다.
train_d2<-sample(dim(Default)[1], dim(Default)[1]*0.7)
test_d2<-Default[-train_d2,]
dlogis_tr2<-glm(default~income+balance+student, data=Default, family=binomial, subset=train_d2)
summary(dlogis_tr2)##
## Call:
## glm(formula = default ~ income + balance + student, family = binomial,
## data = Default, subset = train_d2)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.118e+01 6.164e-01 -18.135 <2e-16 ***
## income 3.688e-06 1.007e-05 0.366 0.7143
## balance 5.895e-03 2.882e-04 20.451 <2e-16 ***
## studentYes -6.334e-01 2.855e-01 -2.219 0.0265 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2010.0 on 6999 degrees of freedom
## Residual deviance: 1069.3 on 6996 degrees of freedom
## AIC: 1077.3
##
## Number of Fisher Scoring iterations: 8dlogis2.pred=rep("No", dim(test_d2)[1])
dlogis2.probs=predict(dlogis_tr2, test_d2, type="response")
dlogis2.pred[dlogis2.probs>0.5]="Yes"
table(dlogis2.pred, test_d2$default)##
## dlogis2.pred No Yes
## No 2880 72
## Yes 15 33mean(dlogis2.pred!=test_d2$default)## [1] 0.029student 변수를 추가했더니 income의 계수가 통계적으로 유의하지 않다고
나왔다.
income과 balance의 추정된 계수가 양수값을 가지기 때문에 두 값이 증가할
수록 채무불이행할 확률이 높고,
student의 추정계수가 음의 값이 나왔기 때문에 학생이 아닌 사람들의
채무불이행 확률이 더 높다는 것을 알 수 있다.
Weekly
- In Section 5.3.2 and 5.3.3, we saw that cv.glm() function can be used in order to compute the LOOCV test error estimate. Alternatively, one could compute those quantities using just the glm() and predict.glm() functions, and a for loop. You will now take this approach in order to compute the LOOCV error for a simple logistic regression model on the Weekly data set. Recall that in the context of classification problems, the LOOCV error is given in (5.4).
library(ISLR2)
library(MASS)
head(Weekly)## Year Lag1 Lag2 Lag3 Lag4 Lag5 Volume Today Direction
## 1 1990 0.816 1.572 -3.936 -0.229 -3.484 0.1549760 -0.270 Down
## 2 1990 -0.270 0.816 1.572 -3.936 -0.229 0.1485740 -2.576 Down
## 3 1990 -2.576 -0.270 0.816 1.572 -3.936 0.1598375 3.514 Up
## 4 1990 3.514 -2.576 -0.270 0.816 1.572 0.1616300 0.712 Up
## 5 1990 0.712 3.514 -2.576 -0.270 0.816 0.1537280 1.178 Up
## 6 1990 1.178 0.712 3.514 -2.576 -0.270 0.1544440 -1.372 Down##변수 설명
Year: The year that the observation was recorded
Lag1: Percentage return for previous week
Lag2: Percentage return for 2 weeks previous
Lag3: Percentage return for 3 weeks previous
Lag4: Percentage return for 4 weeks previous
Lag5: Percentage return for 5 weeks previous
Volume: Volume of shares traded(average number of daily shares traded in
billions)
Today: Percentage return for this week
Direction: A factor with levels Down and Up indicating whether the
market had a positive or negative return on a given week
#기초 통계량
summary(Weekly)## Year Lag1 Lag2 Lag3
## Min. :1990 Min. :-18.1950 Min. :-18.1950 Min. :-18.1950
## 1st Qu.:1995 1st Qu.: -1.1540 1st Qu.: -1.1540 1st Qu.: -1.1580
## Median :2000 Median : 0.2410 Median : 0.2410 Median : 0.2410
## Mean :2000 Mean : 0.1506 Mean : 0.1511 Mean : 0.1472
## 3rd Qu.:2005 3rd Qu.: 1.4050 3rd Qu.: 1.4090 3rd Qu.: 1.4090
## Max. :2010 Max. : 12.0260 Max. : 12.0260 Max. : 12.0260
## Lag4 Lag5 Volume Today
## Min. :-18.1950 Min. :-18.1950 Min. :0.08747 Min. :-18.1950
## 1st Qu.: -1.1580 1st Qu.: -1.1660 1st Qu.:0.33202 1st Qu.: -1.1540
## Median : 0.2380 Median : 0.2340 Median :1.00268 Median : 0.2410
## Mean : 0.1458 Mean : 0.1399 Mean :1.57462 Mean : 0.1499
## 3rd Qu.: 1.4090 3rd Qu.: 1.4050 3rd Qu.:2.05373 3rd Qu.: 1.4050
## Max. : 12.0260 Max. : 12.0260 Max. :9.32821 Max. : 12.0260
## Direction
## Down:484
## Up :605
##
##
##
## Direction은 이번주의 수익률이 양수인지 음수인지를 나타내는 범주형 데이터이고 나머지 변수들은 연속형 데이터이다.
#boxplots
par(mfrow=c(2,4))
boxplot(Year ~ Direction, data =Weekly, main = "Year & Direction")
boxplot(Lag1 ~ Direction, data =Weekly, main = "Lag1 & Direction")
boxplot(Lag2 ~ Direction, data =Weekly, main = "Lag2 & Direction")
boxplot(Lag3 ~ Direction, data =Weekly, main = "Lag3 & Direction")
boxplot(Lag4 ~ Direction, data =Weekly, main = "Lag4 & Direction")
boxplot(Lag5 ~ Direction, data =Weekly, main = "Lag5 & Direction")
boxplot(Today ~ Direction, data =Weekly, main = "Today & Direction")
Direction과 나머지 변수들 간의 boxplot을 그려보면 이번주 수익률을
나타내는 Today변수 외에는 Direction 값에 따른 큰 차이가 보이지
않는다.
# t-test
t.test(Weekly$Year~Weekly$Direction)##
## Welch Two Sample t-test
##
## data: Weekly$Year by Weekly$Direction
## t = 0.73291, df = 1039.3, p-value = 0.4638
## alternative hypothesis: true difference in means between group Down and group Up is not equal to 0
## 95 percent confidence interval:
## -0.4519156 0.9907586
## sample estimates:
## mean in group Down mean in group Up
## 2000.198 1999.929t.test(Weekly$Lag1~Weekly$Direction)##
## Welch Two Sample t-test
##
## data: Weekly$Lag1 by Weekly$Direction
## t = 1.6563, df = 1047.9, p-value = 0.09795
## alternative hypothesis: true difference in means between group Down and group Up is not equal to 0
## 95 percent confidence interval:
## -0.04378476 0.51794261
## sample estimates:
## mean in group Down mean in group Up
## 0.28229545 0.04521653t.test(Weekly$Lag2~Weekly$Direction)##
## Welch Two Sample t-test
##
## data: Weekly$Lag2 by Weekly$Direction
## t = -2.4154, df = 1053.6, p-value = 0.01589
## alternative hypothesis: true difference in means between group Down and group Up is not equal to 0
## 95 percent confidence interval:
## -0.62473558 -0.06467351
## sample estimates:
## mean in group Down mean in group Up
## -0.04042355 0.30428099t.test(Weekly$Lag3~Weekly$Direction)##
## Welch Two Sample t-test
##
## data: Weekly$Lag3 by Weekly$Direction
## t = 0.76067, df = 1058.6, p-value = 0.447
## alternative hypothesis: true difference in means between group Down and group Up is not equal to 0
## 95 percent confidence interval:
## -0.1718491 0.3894400
## sample estimates:
## mean in group Down mean in group Up
## 0.20764669 0.09885124t.test(Weekly$Lag4~Weekly$Direction)##
## Welch Two Sample t-test
##
## data: Weekly$Lag4 by Weekly$Direction
## t = 0.67288, df = 1004.5, p-value = 0.5012
## alternative hypothesis: true difference in means between group Down and group Up is not equal to 0
## 95 percent confidence interval:
## -0.1869635 0.3820924
## sample estimates:
## mean in group Down mean in group Up
## 0.2000207 0.1024562t.test(Weekly$Lag5~Weekly$Direction)##
## Welch Two Sample t-test
##
## data: Weekly$Lag5 by Weekly$Direction
## t = 0.59861, df = 1031.9, p-value = 0.5496
## alternative hypothesis: true difference in means between group Down and group Up is not equal to 0
## 95 percent confidence interval:
## -0.1965830 0.3691748
## sample estimates:
## mean in group Down mean in group Up
## 0.1878347 0.1015388t.test(Weekly$Today~Weekly$Direction)##
## Welch Two Sample t-test
##
## data: Weekly$Today by Weekly$Direction
## t = -33.685, df = 962.59, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group Down and group Up is not equal to 0
## 95 percent confidence interval:
## -3.612545 -3.214797
## sample estimates:
## mean in group Down mean in group Up
## -1.746585 1.667086위의 t-test 결과: Lag1과 Lag2, Today를 제외한 변수들의 p-value값이
0.5와 근사하거나 그 이상이므로 Direction 값에 따른 두 그룹간의 차이가
없다는 귀무가설을 기각할 수 없다.
Lag1과 Lag2, Today의 p-value값은 각각 0.09795, 0.01589,2.2e-16으로 3개
모두 0.1 유의수준에서 귀무가설을 기각할 수 있다.
따라서 로지스틱회귀분석을 할 때 Lag1과 Lag2를 사용한다. (Direction이
Today값이 양수인지 음수인지를 나타낸 것이므로 Today값을 회귀분석에
포함하는 것이 의미없음.)
- Fit a logistic regression model that predicts Direction using Lag1 and Lag2.
wlogis<-glm(Direction~ Lag1+Lag2,data=Weekly,family=binomial)
summary(wlogis)##
## Call:
## glm(formula = Direction ~ Lag1 + Lag2, family = binomial, data = Weekly)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.22122 0.06147 3.599 0.000319 ***
## Lag1 -0.03872 0.02622 -1.477 0.139672
## Lag2 0.06025 0.02655 2.270 0.023232 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1496.2 on 1088 degrees of freedom
## Residual deviance: 1488.2 on 1086 degrees of freedom
## AIC: 1494.2
##
## Number of Fisher Scoring iterations: 4Lag1의 추정된 계수는 유의하지 않다고 나왔고 Lag2의 경우는 0.01수준에서 유의성을 가진다. Lag1은 음수 값을, Lag2는 양수값을 가진 것으로 보아 Lag1 값이 증가할수록 Direction이 Up이 될 확률은 낮아지고, Lag2값이 증가할수록 그 확률이 증가한다.
- Fit a logistic regression model that predicts Direction using Lag1 and Lag2 using all but the first observation.
wlogis<-glm(Direction~ Lag1+Lag2,data=Weekly[-1, ],family=binomial)
summary(wlogis)##
## Call:
## glm(formula = Direction ~ Lag1 + Lag2, family = binomial, data = Weekly[-1,
## ])
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.22324 0.06150 3.630 0.000283 ***
## Lag1 -0.03843 0.02622 -1.466 0.142683
## Lag2 0.06085 0.02656 2.291 0.021971 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1494.6 on 1087 degrees of freedom
## Residual deviance: 1486.5 on 1085 degrees of freedom
## AIC: 1492.5
##
## Number of Fisher Scoring iterations: 4- Use the model from (b) to predict the direction of the first observation. You can do this by predicting that the first observation will go up if P(Direction=“Up”|Lag1, Lag2)>0.5. Was this observation correctly classified?
ifelse(predict(wlogis,Weekly[1,],type="response")>.5, "Up","Down")## 1
## "Up"print(Weekly$Direction[1])## [1] Down
## Levels: Down Up첫번째 관측치를 제외한 데이터로 로지스틱 회귀모델을 학습하고 첫번째 관측치를 예측해보았다. 예측 결과는 “Up”으로 나왔지만, 첫 번째 값의 실제값은 “Down”이다. 즉, 예측이 틀렸음을 알 수 있다.
- Write a for loop from \(i\)=1 to \(i\)=n, where n is the number of observations in the data set, that performs each of the following steps:
- Fit a logistic regression model using all but the \(i\)th observation to predict Direction
using Lag1 and Lag2.
- Compute the posterior probability of the market moving up for the
\(i\)th observation.
- Use the posterior probability for the \(i\)th observation in order to predict
whether or not the market moves up.
- Determine whether or not an error was made in predicting the direction for the \(i\)th observation. If an error was made, then indicate this as a 1, and otherwise indicate it as a 0.
pred_errors <-c()
for (i in 1:nrow(Weekly)) {
wlogis_i <- glm(Direction ~ Lag1 + Lag2, data = Weekly[-i, ], family = binomial)
wlogis_i.pred <- ifelse(predict(wlogis_i, Weekly[i, ],type="response")>.5,"Up","Down")
pred_errors[i]<-as.numeric(wlogis_i.pred!=Weekly[i,"Direction"])#예측값이 실제값과 다르면 1, 같으면 0
}
pred_errors## [1] 1 1 0 1 0 1 0 0 0 1 1 0 1 0 1 0 1 0 1 0 0 0 1 1 1 1 1 1 0 1 1 1 1 0 1 0 0
## [38] 0 1 0 1 0 0 1 0 1 1 1 0 1 0 0 0 1 0 0 1 1 0 0 0 0 1 0 1 1 0 0 1 0 1 1 0 0
## [75] 0 1 0 1 1 0 0 1 1 0 1 1 0 0 1 0 0 1 1 1 0 0 0 0 0 1 0 1 1 0 0 1 0 1 0 0 1
## [112] 1 0 0 1 0 0 1 0 0 1 1 1 1 0 0 0 1 0 1 0 1 1 0 0 0 1 1 1 0 0 0 1 0 0 0 0 0
## [149] 0 1 1 1 0 1 0 0 1 1 0 1 0 0 1 1 0 0 1 0 0 1 0 0 1 1 1 0 1 0 1 0 0 0 0 0 0
## [186] 0 0 1 1 0 1 0 1 0 1 0 1 0 0 1 0 0 1 0 0 1 0 1 0 1 1 1 0 0 1 1 0 1 0 0 1 1
## [223] 0 0 0 1 1 1 0 1 0 1 0 1 0 0 0 1 1 0 1 0 1 0 1 0 1 0 1 1 0 1 0 0 1 0 0 1 0
## [260] 0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 1 0 0 1 0 1 1 0 0 0 0 0 1 0 1 0
## [297] 0 1 0 0 0 1 0 0 1 1 0 0 1 0 0 0 0 1 0 1 1 0 0 1 0 1 0 1 1 0 0 0 1 0 1 0 0
## [334] 1 1 1 1 0 1 0 0 1 0 0 0 1 0 1 0 1 0 0 0 0 0 1 1 0 0 1 0 0 1 0 0 0 1 1 0 1
## [371] 1 1 1 1 0 0 0 1 0 0 0 0 0 0 1 0 1 1 0 0 1 1 0 0 0 0 0 1 0 0 1 1 1 0 1 0 1
## [408] 0 1 1 1 0 1 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 0 1 0 1 0 0 0 0 0 1 1 1 1
## [445] 0 1 1 0 1 0 1 1 0 1 0 0 1 0 0 1 1 0 0 0 0 1 1 0 0 1 0 1 0 0 0 1 0 0 1 0 0
## [482] 0 1 1 1 0 1 0 0 0 1 0 1 1 1 0 0 0 0 1 1 1 0 1 1 0 1 0 0 0 1 0 1 0 0 0 1 0
## [519] 1 1 0 0 1 1 0 0 0 1 1 0 1 0 1 1 1 1 1 0 0 0 1 0 0 0 1 1 0 1 0 0 0 1 1 1 1
## [556] 1 1 0 1 0 1 0 0 1 0 0 1 1 1 0 0 0 1 1 1 1 1 1 1 1 1 1 0 1 0 0 1 0 0 1 0 1
## [593] 0 0 1 0 0 1 0 1 1 0 1 1 1 0 1 0 1 0 1 0 0 0 1 0 1 0 1 0 1 1 0 1 1 0 1 0 1
## [630] 0 1 1 1 1 0 1 1 0 0 0 1 1 1 1 0 1 1 1 0 1 0 0 0 1 1 1 1 1 1 0 1 0 0 1 0 0
## [667] 0 1 1 0 1 0 1 1 1 1 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 1 0 1 1
## [704] 0 0 0 0 1 0 1 0 1 0 1 0 0 1 1 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 1 1 1 0 1 1 0
## [741] 1 1 1 1 0 0 0 1 1 1 1 1 1 0 1 0 0 0 0 0 0 1 0 1 1 1 0 0 0 1 0 0 1 0 0 0 1
## [778] 1 1 0 0 0 1 0 0 1 1 1 0 0 1 0 1 0 1 0 0 1 0 0 1 0 0 0 0 0 1 0 1 1 0 0 1 1
## [815] 0 1 1 1 0 0 0 0 0 1 1 0 0 1 0 0 1 0 1 0 0 0 1 1 0 1 1 0 1 0 1 0 1 1 0 0 1
## [852] 1 1 0 1 1 0 0 0 1 0 1 0 1 0 1 0 0 0 0 0 1 0 0 1 1 0 0 1 0 1 0 1 1 0 1 0 1
## [889] 1 0 0 0 1 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 1 1 1 1 1 0 1 1 0 0 0 0 0 1 0 0 1
## [926] 0 0 0 0 1 0 1 1 1 0 0 1 1 0 1 1 1 1 0 1 0 1 0 1 0 1 0 1 0 0 1 1 1 1 1 0 1
## [963] 0 0 0 1 1 1 0 1 1 1 1 0 0 0 0 1 1 0 0 0 0 1 0 0 1 1 1 0 0 1 1 1 0 1 0 0 0
## [1000] 0 1 0 0 1 0 1 0 0 1 1 1 1 0 1 0 0 1 0 0 1 0 0 1 1 0 1 1 1 0 1 1 0 0 0 1 0
## [1037] 1 0 1 1 1 1 0 0 1 0 0 0 0 0 0 1 0 1 1 0 0 0 0 0 1 1 1 0 0 0 1 0 1 1 1 0 0
## [1074] 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0LOOCV방법으로 학습 후 예측을 해본 결과 위의 결과가 나왔다. 예측값과 실제값이 다르면 1, 같으면 0이 출력된다.
- Take the average of the n numbers obtained in (d) 4) in order to obtain the LOOCV estimate for the test error. Comment on the results.
mean(pred_errors)## [1] 0.4499541error rate이 44.995%로 나왔다.