head(Wage)
## year age maritl race education region
## 231655 2006 18 1. Never Married 1. White 1. < HS Grad 2. Middle Atlantic
## 86582 2004 24 1. Never Married 1. White 4. College Grad 2. Middle Atlantic
## 161300 2003 45 2. Married 1. White 3. Some College 2. Middle Atlantic
## 155159 2003 43 2. Married 3. Asian 4. College Grad 2. Middle Atlantic
## 11443 2005 50 4. Divorced 1. White 2. HS Grad 2. Middle Atlantic
## 376662 2008 54 2. Married 1. White 4. College Grad 2. Middle Atlantic
## jobclass health health_ins logwage wage
## 231655 1. Industrial 1. <=Good 2. No 4.318063 75.04315
## 86582 2. Information 2. >=Very Good 2. No 4.255273 70.47602
## 161300 1. Industrial 1. <=Good 1. Yes 4.875061 130.98218
## 155159 2. Information 2. >=Very Good 1. Yes 5.041393 154.68529
## 11443 2. Information 1. <=Good 1. Yes 4.318063 75.04315
## 376662 2. Information 2. >=Very Good 1. Yes 4.845098 127.11574
str(Wage)
## 'data.frame': 3000 obs. of 11 variables:
## $ year : int 2006 2004 2003 2003 2005 2008 2009 2008 2006 2004 ...
## $ age : int 18 24 45 43 50 54 44 30 41 52 ...
## $ maritl : Factor w/ 5 levels "1. Never Married",..: 1 1 2 2 4 2 2 1 1 2 ...
## $ race : Factor w/ 4 levels "1. White","2. Black",..: 1 1 1 3 1 1 4 3 2 1 ...
## $ education : Factor w/ 5 levels "1. < HS Grad",..: 1 4 3 4 2 4 3 3 3 2 ...
## $ region : Factor w/ 9 levels "1. New England",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ jobclass : Factor w/ 2 levels "1. Industrial",..: 1 2 1 2 2 2 1 2 2 2 ...
## $ health : Factor w/ 2 levels "1. <=Good","2. >=Very Good": 1 2 1 2 1 2 2 1 2 2 ...
## $ health_ins: Factor w/ 2 levels "1. Yes","2. No": 2 2 1 1 1 1 1 1 1 1 ...
## $ logwage : num 4.32 4.26 4.88 5.04 4.32 ...
## $ wage : num 75 70.5 131 154.7 75 ...
summary(Wage)
## year age maritl race
## Min. :2003 Min. :18.00 1. Never Married: 648 1. White:2480
## 1st Qu.:2004 1st Qu.:33.75 2. Married :2074 2. Black: 293
## Median :2006 Median :42.00 3. Widowed : 19 3. Asian: 190
## Mean :2006 Mean :42.41 4. Divorced : 204 4. Other: 37
## 3rd Qu.:2008 3rd Qu.:51.00 5. Separated : 55
## Max. :2009 Max. :80.00
##
## education region jobclass
## 1. < HS Grad :268 2. Middle Atlantic :3000 1. Industrial :1544
## 2. HS Grad :971 1. New England : 0 2. Information:1456
## 3. Some College :650 3. East North Central: 0
## 4. College Grad :685 4. West North Central: 0
## 5. Advanced Degree:426 5. South Atlantic : 0
## 6. East South Central: 0
## (Other) : 0
## health health_ins logwage wage
## 1. <=Good : 858 1. Yes:2083 Min. :3.000 Min. : 20.09
## 2. >=Very Good:2142 2. No : 917 1st Qu.:4.447 1st Qu.: 85.38
## Median :4.653 Median :104.92
## Mean :4.654 Mean :111.70
## 3rd Qu.:4.857 3rd Qu.:128.68
## Max. :5.763 Max. :318.34
##
par(mfrow=c(1,3))
barplot(table(Wage$maritl), main="maritl")
barplot(table(Wage$race), main="race")
barplot(table(Wage$education), main="education")
barplot(table(Wage$region), main="region")
barplot(table(Wage$jobclass), main="jobclass")
barplot(table(Wage$health), main="health")
barplot(table(Wage$health_ins), main="health_ins")
hist(Wage$year, main="year",freq=F)
lines(density(Wage$year))
hist(Wage$age, main="age",freq=F)
lines(density(Wage$age))
hist(Wage$logwage, main="logwage",freq=F)
lines(density(Wage$logwage))
hist(Wage$wage, main="wage",freq=F)
lines(density(Wage$wage))
pairs(Wage[,c(1,2,10,11)])
corr_w<-cor(Wage[,c(1,2,10,11)])
round(corr_w,2)
## year age logwage wage
## year 1.00 0.04 0.08 0.07
## age 0.04 1.00 0.22 0.20
## logwage 0.08 0.22 1.00 0.95
## wage 0.07 0.20 0.95 1.00
corrplot(corr_w,method="number", type='lower')
5. 수치형변수상관계수(중복없이 작성) 거의 상관관계가 없음.
boxplot(wage ~ maritl, data = Wage, main = "wage & maritl")
#t.test(Wage$wage~Wage$maritl)
boxplot(wage ~ race, data = Wage, main = "wage & race")
#t.test(Wage$wage~Wage$race)
boxplot(wage ~ education, data = Wage, main = "wage & education")
#t.test(Wage$wage~Wage$education)
boxplot(wage ~ jobclass, data = Wage, main = "wage & jobclass")
t.test(Wage$wage~Wage$jobclass)
##
## Welch Two Sample t-test
##
## data: Wage$wage by Wage$jobclass
## t = -11.489, df = 2714.9, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group 1. Industrial and group 2. Information is not equal to 0
## 95 percent confidence interval:
## -20.21940 -14.32378
## sample estimates:
## mean in group 1. Industrial mean in group 2. Information
## 103.3211 120.5927
boxplot(wage ~ health, data = Wage, main = "wage & health")
t.test(Wage$wage~Wage$health)
##
## Welch Two Sample t-test
##
## data: Wage$wage by Wage$health
## t = -9.2265, df = 1934.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group 1. <=Good and group 2. >=Very Good is not equal to 0
## 95 percent confidence interval:
## -17.05452 -11.07524
## sample estimates:
## mean in group 1. <=Good mean in group 2. >=Very Good
## 101.6613 115.7262
boxplot(wage ~ health_ins, data = Wage, main = "wage & health_ins")
t.test(Wage$wage~Wage$health_ins)
##
## Welch Two Sample t-test
##
## data: Wage$wage by Wage$health_ins
## t = 18.708, df = 1989.5, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group 1. Yes and group 2. No is not equal to 0
## 95 percent confidence interval:
## 24.99464 30.84858
## sample estimates:
## mean in group 1. Yes mean in group 2. No
## 120.2383 92.3167
In this exercise, we will predict the number of applications received using the the other variables in the College data set.
Wage1<-Wage[,-c(6,10)]#region, logwage 값 제외
set.seed(1)
train_w<-sample(dim(Wage1)[1], dim(Wage1)[1]*0.7)
test_w<- (-train_w)
wage_train<-Wage1[train_w,]
wage_test<-Wage1[test_w,]
selec_for <- regsubsets(wage ~ ., data = wage_train, nvmax = ncol(Wage1)-1, method = "forward")
selec_sum <- summary(selec_for)
selec_sum
## Subset selection object
## Call: regsubsets.formula(wage ~ ., data = wage_train, nvmax = ncol(Wage1) -
## 1, method = "forward")
## 16 Variables (and intercept)
## Forced in Forced out
## year FALSE FALSE
## age FALSE FALSE
## maritl2. Married FALSE FALSE
## maritl3. Widowed FALSE FALSE
## maritl4. Divorced FALSE FALSE
## maritl5. Separated FALSE FALSE
## race2. Black FALSE FALSE
## race3. Asian FALSE FALSE
## race4. Other FALSE FALSE
## education2. HS Grad FALSE FALSE
## education3. Some College FALSE FALSE
## education4. College Grad FALSE FALSE
## education5. Advanced Degree FALSE FALSE
## jobclass2. Information FALSE FALSE
## health2. >=Very Good FALSE FALSE
## health_ins2. No FALSE FALSE
## 1 subsets of each size up to 8
## Selection Algorithm: forward
## year age maritl2. Married maritl3. Widowed maritl4. Divorced
## 1 ( 1 ) " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " "
## 3 ( 1 ) " " " " "*" " " " "
## 4 ( 1 ) " " " " "*" " " " "
## 5 ( 1 ) " " " " "*" " " " "
## 6 ( 1 ) " " "*" "*" " " " "
## 7 ( 1 ) " " "*" "*" " " " "
## 8 ( 1 ) " " "*" "*" " " " "
## maritl5. Separated race2. Black race3. Asian race4. Other
## 1 ( 1 ) " " " " " " " "
## 2 ( 1 ) " " " " " " " "
## 3 ( 1 ) " " " " " " " "
## 4 ( 1 ) " " " " " " " "
## 5 ( 1 ) " " " " " " " "
## 6 ( 1 ) " " " " " " " "
## 7 ( 1 ) " " " " " " " "
## 8 ( 1 ) " " " " " " " "
## education2. HS Grad education3. Some College education4. College Grad
## 1 ( 1 ) " " " " " "
## 2 ( 1 ) " " " " " "
## 3 ( 1 ) " " " " " "
## 4 ( 1 ) " " " " "*"
## 5 ( 1 ) " " "*" "*"
## 6 ( 1 ) " " "*" "*"
## 7 ( 1 ) " " "*" "*"
## 8 ( 1 ) "*" "*" "*"
## education5. Advanced Degree jobclass2. Information
## 1 ( 1 ) "*" " "
## 2 ( 1 ) "*" " "
## 3 ( 1 ) "*" " "
## 4 ( 1 ) "*" " "
## 5 ( 1 ) "*" " "
## 6 ( 1 ) "*" " "
## 7 ( 1 ) "*" " "
## 8 ( 1 ) "*" " "
## health2. >=Very Good health_ins2. No
## 1 ( 1 ) " " " "
## 2 ( 1 ) " " "*"
## 3 ( 1 ) " " "*"
## 4 ( 1 ) " " "*"
## 5 ( 1 ) " " "*"
## 6 ( 1 ) " " "*"
## 7 ( 1 ) "*" "*"
## 8 ( 1 ) "*" "*"
cat("Cp:",which.min(selec_sum$cp),"\n BIC:",which.min(selec_sum$bic),"\n Adj R^2:",which.min(selec_sum$adjr2))
## Cp: 8
## BIC: 7
## Adj R^2: 1
selec_plot <- function(stat, y.label, adj2 = FALSE) {
plot(stat, xlab = "Number of Variables", ylab = y.label, xaxt = "n", type = "l")
axis(side = 1, at = 1:length(stat))
if (adj2==FALSE) {
#가장 작은 값에서 한 단위의 표준오차 내의 값들의 집합
stat_1se <- min(stat) + (sd(stat) / sqrt(length(stat)))
sub <- which(stat< stat_1se)
}
else{
#가장 큰 값에서 한 단위의 표준오차 내의 값들의 집합
stat_1se <- max(stat) - (sd(stat) / sqrt(length(stat)))
sub <- which(stat > stat_1se)
}
#한 단위의 표준오차 내의 값을 가지는 것들 중 변수의 개수가 가장 작은 것 표시
abline(h = stat_1se, col = "red", lty = 2)
abline(v = sub[1], col = "green", lty = 2)
}
par(mfrow=c(1, 3))
selec_plot(selec_sum$cp, "Cp")
selec_plot(selec_sum$bic, "BIC")
selec_plot(selec_sum$adjr2, "Adjusted R2", adj2 = TRUE) # higher values are better
coef(selec_for,5)#추정계수
## (Intercept) maritl2. Married
## 89.90879 17.92762
## education3. Some College education4. College Grad
## 11.97367 23.48902
## education5. Advanced Degree health_ins2. No
## 50.09033 -21.37292
plot(selec_for,scale="Cp")
plot(selec_for,scale="bic")
plot(selec_for,scale="adjr2")
coef(selec_for,5)
## (Intercept) maritl2. Married
## 89.90879 17.92762
## education3. Some College education4. College Grad
## 11.97367 23.48902
## education5. Advanced Degree health_ins2. No
## 50.09033 -21.37292
m_w<-lm(wage~maritl+education+health_ins,data=Wage1,subset=train_w)
summary(m_w)
##
## Call:
## lm(formula = wage ~ maritl + education + health_ins, data = Wage1,
## subset = train_w)
##
## Residuals:
## Min 1Q Median 3Q Max
## -98.930 -18.902 -2.884 13.902 214.963
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 80.828 2.991 27.026 < 2e-16 ***
## maritl2. Married 20.235 1.818 11.133 < 2e-16 ***
## maritl3. Widowed 7.402 10.192 0.726 0.46781
## maritl4. Divorced 6.254 3.205 1.952 0.05110 .
## maritl5. Separated 14.173 5.539 2.559 0.01057 *
## education2. HS Grad 8.053 2.793 2.883 0.00397 **
## education3. Some College 18.539 2.958 6.267 4.45e-10 ***
## education4. College Grad 30.233 2.932 10.311 < 2e-16 ***
## education5. Advanced Degree 56.800 3.204 17.726 < 2e-16 ***
## health_ins2. No -20.652 1.635 -12.634 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 33.38 on 2090 degrees of freedom
## Multiple R-squared: 0.3299, Adjusted R-squared: 0.327
## F-statistic: 114.3 on 9 and 2090 DF, p-value: < 2.2e-16
m_w.pred<-predict(m_w,wage_test)
m_w.error<-mean((m_w.pred-wage_test$wage)^2)
m_w.error
## [1] 1339.502
coef(selec_for,1)
## (Intercept) education5. Advanced Degree
## 104.76017 45.49489