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. :73554#boxplots
par(mfrow=c(1,3))
boxplot(balance ~ default, data = Default, main = "balance & default")
boxplot(income ~ default, data = Default, main = "income & default")
#student&default 분할표 만들어보기
#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.15
- Fit a logistic regression model that uses income and balance to predict default.
dlogis<-glm(default~income+balance, data= Default, family=binomial)- 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.0455- 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_d<-sample(dim(Default)[1], dim(Default)[1]/2)
dlogis_tr<-glm(default~income+balance, data=Default, family=binomial, subset=train_d)
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 50mean(dlogis.pred!=Default$default)## [1] 0.0455train_d<-sample(dim(Default)[1], dim(Default)[1]/2)
dlogis_tr<-glm(default~income+balance, data=Default, family=binomial, subset=train_d)
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 50mean(dlogis.pred!=Default$default)## [1] 0.0455train_d<-sample(dim(Default)[1], dim(Default)[1]/2)
dlogis_tr<-glm(default~income+balance, data=Default, family=binomial, subset=train_d)
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 50mean(dlogis.pred!=Default$default)## [1] 0.0455- 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.
train_d<-sample(dim(Default)[1], dim(Default)[1]/2)
dlogis_tr<-glm(default~., data=Default, family=binomial, subset=train_d)
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 50mean(dlogis.pred!=Default$default)## [1] 0.0455error rate이 더 높아짐.
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 fiven 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
##
##
##
## #scatter plot
pairs(Weekly)
#상관계수
cor(Weekly[,-9])## Year Lag1 Lag2 Lag3 Lag4
## Year 1.00000000 -0.032289274 -0.03339001 -0.03000649 -0.031127923
## Lag1 -0.03228927 1.000000000 -0.07485305 0.05863568 -0.071273876
## Lag2 -0.03339001 -0.074853051 1.00000000 -0.07572091 0.058381535
## Lag3 -0.03000649 0.058635682 -0.07572091 1.00000000 -0.075395865
## Lag4 -0.03112792 -0.071273876 0.05838153 -0.07539587 1.000000000
## Lag5 -0.03051910 -0.008183096 -0.07249948 0.06065717 -0.075675027
## Volume 0.84194162 -0.064951313 -0.08551314 -0.06928771 -0.061074617
## Today -0.03245989 -0.075031842 0.05916672 -0.07124364 -0.007825873
## Lag5 Volume Today
## Year -0.030519101 0.84194162 -0.032459894
## Lag1 -0.008183096 -0.06495131 -0.075031842
## Lag2 -0.072499482 -0.08551314 0.059166717
## Lag3 0.060657175 -0.06928771 -0.071243639
## Lag4 -0.075675027 -0.06107462 -0.007825873
## Lag5 1.000000000 -0.05851741 0.011012698
## Volume -0.058517414 1.00000000 -0.033077783
## Today 0.011012698 -0.03307778 1.000000000- 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: 4- 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,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: 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?
wlogis.probs<-predict(wlogis, type="response")
wlogis.pred<-rep("Down",1089)
wlogis.pred[wlogis.probs>.5]<-"Up"
table(wlogis.pred,Weekly$Direction)##
## wlogis.pred Down Up
## Down 38 38
## Up 446 567#정확도
mean(wlogis.pred==Weekly$Direction)## [1] 0.5555556# Up일 때의 정확도
557/(430+557)## [1] 0.5643364#Down일 때의 정확도
54/(54+48)## [1] 0.5294118전체 정확도는 56.1%이고 Up을 예측할 때의 정확도는 56.4%로 Down을 예측할 때의 정확도인 53%보다 정확도가 더 높다.
- 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 ith observation to predict Direction using Lag1 and Lag2.
weekly_train=(Weekly$Year<2009)
wlogis_split=glm(Direction~ Lag2, data=Weekly, subset=weekly_train,family="binomial")
summary(wlogis_split)##
## Call:
## glm(formula = Direction ~ Lag2, family = "binomial", data = Weekly,
## subset = weekly_train)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.20326 0.06428 3.162 0.00157 **
## Lag2 0.05810 0.02870 2.024 0.04298 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1354.7 on 984 degrees of freedom
## Residual deviance: 1350.5 on 983 degrees of freedom
## AIC: 1354.5
##
## Number of Fisher Scoring iterations: 4- Compute the posterior probability of the market moving up for the ith observation.
wlogis_split.probs<-predict(wlogis_split,Weekly[!weekly_train,], type="response")
wlogis_split.pred<-rep("Down",dim(Weekly[!weekly_train,])[1])
wlogis_split.pred[wlogis_split.probs>.5]<-"Up"
table(wlogis_split.pred,Weekly[!weekly_train, ]$Direction)##
## wlogis_split.pred Down Up
## Down 9 5
## Up 34 56- Use the posterior probability for the ith observation in order to predict whether or not the market moves up.
#정확도
mean(wlogis_split.pred==Weekly[!weekly_train, ]$Direction)## [1] 0.625# Up일 때의 정확도
56/(56+34)## [1] 0.6222222#Down일 때의 정확도
9/(9+5)## [1] 0.6428571- Determine whether or not an error was made in predicting the direction fot the ith observation. If an error was made, then indicate this as a 1, and otherwise indicate it as a 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.
wlda<-lda(Direction~ Lag2, data=Weekly, subset=weekly_train)
wlda## Call:
## lda(Direction ~ Lag2, data = Weekly, subset = weekly_train)
##
## Prior probabilities of groups:
## Down Up
## 0.4477157 0.5522843
##
## Group means:
## Lag2
## Down -0.03568254
## Up 0.26036581
##
## Coefficients of linear discriminants:
## LD1
## Lag2 0.4414162wlda.pred<-predict(wlda,Weekly[!weekly_train,])
table(wlda.pred$class,Weekly[!weekly_train,]$Direction)##
## Down Up
## Down 9 5
## Up 34 56#정확도
mean(wlda.pred$class==Weekly[!weekly_train, ]$Direction)## [1] 0.625