투빅스 6&7기 2주차 과제 - 회귀분석 6기 장재석

rm(list=ls())

library(car)

setwd("C:\\Users\\JJS810\\Desktop")

student <-read.csv(list.files()[6],header=T)

View(student) 

str(student)

#Factor -> dummy variable

#plot & basic summary

summary(student) #missing value : 2

par(mfrow=c(2,2))

for (i in 6:9){

  plot(student[,i],student$target,main=colnames(student)[i])

}

#target과 독립변수간의 상관관계가 정확히 보이지 않는다. 

#numeric한 변수와 target간의 산점도

#correlation

colnames(student)

cov.wt(student[,c(6:9,11)],cor=T)$cor

#target과 numeric한 변수들간에 약한 상관관계가 보인다.

#0.275~0.502 정도의 양의 상관관계가 보인다.

#how to find missing value

sum(is.na(student)) #missing value : 2

colSums(is.na(student)) #Topic : 1, ParentschoolSatisfaction : 1

#remove the missing value

which(is.na(student$Topic)) #78

which(is.na(student$ParentschoolSatisfaction)) #87

student$Topic[78];student$ParentschoolSatisfaction[87]

#missing value -> omit

student <- na.omit(student) #389obs -> 387obs

#dummy -> factor

#model.matrix()

#dummy variable about 2 levels

student$gender <- ifelse(student$gender=="F",1,0) #F=1, M=0

student$semester <- ifelse(student$semester=="F",1,0) #F=1, S=0

student$responsible <- ifelse(student$responsible=="Father",1,0) #Father=1, Mum=0

student$ParentAnsweringSurvey <-ifelse(student$ParentAnsweringSurvey=="Yes",1,0) #Yes=1 ,No=0

student$StudentAbsenceDays <- ifelse(student$StudentAbsenceDays=="Above-7",1,0) #Above-7=1, Under-7=0

student$ParentschoolSatisfaction <-ifelse(student$ParentschoolSatisfaction=="Good",1,0) #Good=1, bad=0

#duumy variable about 3, 8 levels

colnames(student)

section <- matrix(0,nrow=387,ncol=2)

colnames(section) <- paste0("dumsec",letters[1:2])

for (i in 1:387 ) {

  if(student[i,2]=="A"){section[i,1]<-1}

  else if(student[i,2]=="B"){section[i,2]<-1}

  else{}

}

levels(student$Topic) #8 levels

#"Arabic"  "Biology" "English" "French"  "History" "IT"     

#"Science" "Spanish"

topic <- matrix(0,nrow=387,ncol=7)

colnames(topic)<-paste0("dum",c("arab","bio","eng","fre","his","it","Sci"))

for (j in 1:387){

  if(student[j,3]=="Arabic"){topic[j,1]<-1}

  else if(student[j,3]=="Biology"){topic[j,2]<-1}

  else if(student[j,3]=="English"){topic[j,3]<-1}

  else if(student[j,3]=="French"){topic[j,4]<-1}

  else if(student[j,3]=="History"){topic[j,5]<-1}

  else if(student[j,3]=="IT"){topic[j,6]<-1}

  else if(student[j,3]=="Science"){topic[j,7]<-1}

  else{}

}

#merge data

student <- student[,-c(2,3)]

student <- cbind(student,section, topic)

str(student)

#regression

fit <-lm(target~.,data=student) #full model

summary(fit) #0.6097

vif(fit) #sqrt(vif)>2 : multicolinearity

vif(fit)[sqrt(vif(fit))>2] #dummy라 제거하기 애매하고 4~5사이이므로 그냥 놔두자.

#일반적으로 다중공선성의 경우 VIF가 10이하이면 없다고 판단할 수 있으므로, 공선성 문제는 없다고 생각하자.

#stepwise by AIC

library(MASS)

fit.aic<-stepAIC(fit,direction='both')

fit.aic$anova

fit1 <- lm(target ~ semester + responsible + raisehand + VisITedResources + 

             AnnouncementsView + ParentAnsweringSurvey + StudentAbsenceDays + 

             dumSci, data=student)

summary(fit1) #0.6151

plot(fit1) #42,46,329 이상치이고 200후반대부터 살짝 qqline을 벗어나는 영향을 보인다.

#정규성 여부 검정 - 단일 모집단 분포의 정규성 검정

par(mfrow=c(1,1))

qqnorm(fit1$residuals) #잔차를 확인해야 한다.

hist(fit1$residuals)

#살짝 qqline이 정규성을 약간 만족하지 않는 것으로 보인다.

#정확한 회귀진단은 뒤에서 자세히 다루겠다.

#forward

fit.con <-lm(target~1,data=student)

fit.full <- lm(target~.,data=student)

fit.forward <- step(fit.con, scope=list(lower=fit.con, upper=fit.full),direction="forward")

summary(fit.forward) #0.6151

#backward

fit.backward <- step(fit.full, scope=list(lower=fit.con, upper=fit.full),direction="backward")

summary(fit.backward) #0.6151

#회귀진단

diag.student <-ls.diag(fit1)

diag.student

#정규성검정

qqPlot(fit1,labels=row.names(student),id.metho="identify",simulate=TRUE,main="Q-Q PLOT")

#leverage = 이상치 , dffits,cooks = 영향치 

myregression<-function(x){

  h=hat(model.matrix(x))  

  d=dffits(x) 

  cutoff=4/((nrow(student)-length(x$coefficients)-2))

  par(mfrow=c(1,3))

  plot(h,type="h",main="leverage plot",ylab="leverage")

  plot(d,type="h",main = "dffits",ylab="dffits")

  plot(x,which=4,cook.levels = cutoff)

  par(mfrow=c(1,1))

}

outlierTest(fit1, cutoff=0.5) #42, 46 이상치 

#샤피로테스트,더빈왓슨테스트, 등분산성 테스트 

shapiro.test(residuals(fit1)) #통계량에 의한 정규성 검정, p-value :0.04698 -> reject H0 

#H0 : 모집단은 정규분포를 따른다...샤피로 테스트에서 귀무가설을 기각하므로 모집단이 정규성을 만족하지 않는다.

#즉 정규성 부분에서 문제가 발생한다....이 문제를 어떻게 해결할까?

durbinWatsonTest(fit1) #종속변수의 독립성 p값이 의미가 없어야지 자기상관이 없는것 

crPlots(fit1) #변수들의 선형성을 확인하였고, fitting 된 모델이 어느정도 잘 맞는다.

ncvTest(fit1) #등분산성을 확인하는 테스트이고, p=0.2245947로 유의하지 않으므로 오차의 등분산성 만족 

spreadLevelPlot(fit1) #fitted vakye에 대한 스튜던트화 잔차를 시각화 

#아래의 그림에서 기존의 빨간선보다 잘 맞지 않는 것을 알 수 있다.

residplot <- function(fit, nbreaks=10) {

  z <- rstudent(fit)

  hist(z, breaks=nbreaks, freq=FALSE,xlab="Studentized Residual",

       main="Distribution of Errors")

  rug(jitter(z), col="pink")

  curve(dnorm(x, mean=mean(z), sd=sd(z)),add=TRUE, col="blue", lwd=2)

  lines(density(z)$x, density(z)$y,col="red", lwd=2, lty=2)

  legend("topright",legend = c( "Normal Curve", "Kernel Density Curve"),

         lty=1:2, col=c("blue","red"), cex=.7)

}

residplot(fit1) #정규분포 어느정도 차이가 보인다. 

#qqnorm으로 확인 했을때, 정규성이 잘 안 맞는 것을 밑의 그림으로 한번 더 확인 할 수 있다.

#선형모형 가정에 대한 전반적 검증

install.packages("gvlma")

library(gvlma)

gvmodel <- gvlma(fit1)

summary(gvmodel) #OLS 가정 만족한다.!!

#그런데 왜 정규성을 만족안할까? -> 혹시 이상치가 문제가 있는지 한번 더 확인해보자 

#high leverage points

#위에서 정규성이 안 맞는 이유에서

library(ggplot2)

library(ggiraph)

library(moonBook)

par(mfrow=c(1,1))

hat.plot <- function(fit) {

  p <- length(coefficients(fit))

  n <- length(fitted(fit))

  y=hatvalues(fit)

  name=attr(hatvalues(fit),"names")

  df=data.frame(x=1:length(y),y=as.numeric(y),name=name)

  p1<-ggplot(df,aes(x=x,y=y,tooltip=name,data_id=x))+geom_point_interactive()

  yintercept2=2*p/n

  p1<-p1+geom_hline(aes(yintercept=yintercept2),col="orange",lty="dashed")

  yintercept3=3*p/n

  p1<-p1+geom_hline(aes(yintercept=yintercept3),col="red",lty="dashed")

  ggiraph(code=print(p1))

  

}                 

hat.plot(fit1) #red line : average hat보다 2배, 3배 높은 것을 의미

#13,169,218,222,226,342 가 평균 hat보다 높게 있다.

#42,46,329 이상치

#위의 회귀진단을 바탕으로 관측치를 제거해보자.

#주의해야할 사항은 이상관측치가 흥미로운 데이터일수도 있다는 점

#따라서 원래는 이상치를 일일히 확인을 해봐야 한다.

#이상관측치 인덱스

student[c(13,42,46,169,218,222,226,329,342,354),]

plot(student$raisehand,student$target)

lines(smooth.spline(x=student$raisehand, y=student$target))

#여기서 이상치를 일일히 확인해서, 그 이유를 찾기에는 이 데이터에 대한 정보가 부족하다고 생각

#따라서 그냥 이상치를 제거하고 회귀분석을 진행하겠다.

#이상치를 제거한 모형 

student.rm <-student[-c(13,42,46,169,218,222,226,329,342,354),]

str(student.rm)

fit2 <-lm(target~.,data=student.rm)

fit.step2<- stepAIC(fit2,direction="both")

fit3 <- lm(target ~ responsible + raisehand + VisITedResources + ParentAnsweringSurvey + 

             StudentAbsenceDays + dumeng + dumit + dumSci,data=student.rm)

summary(fit3) #0.6369로 증가하고 모든 독립변수가 유의하다. 

#종속변수를 한번 표준화해보고 같은 독립변수에 회귀분석을 해보았다( 정규성문제때문에 한번해보았다.)

fit4 <- lm(scale(target) ~ responsible + raisehand + VisITedResources + ParentAnsweringSurvey + 

             StudentAbsenceDays + dumeng + dumit + dumSci,data=student.rm)

summary(fit4) #0.6369 : 위의 fit3와 같은 결과가 나온다.

par(mfrow=c(2,2))

plot(fit3)

vif(fit3) #다중공선성은 없다.

#위의 회귀진단을 다시 해보았다.

myregression(fit3) #20, 73, 329가 다시 이상치로 나온다. 

qqPlot(fit3,labels=row.names(student),id.method="identify",simulate=TRUE,main="Q-Q PLOT")

#여전히 Q-Q plot에서 약간씩 데이터가 정규성을 만족하지 않는다.

shapiro.test(residuals(fit3)) #통계량에 의한 정규성 검정, p-value :0.04334 -> reject H0 

durbinWatsonTest(fit3) #종속변수의 독립성 p값이 의미가 없어야지 자기상관이 없는것 

crPlots(fit3) #선형성 시각화

ncvTest(fit3) #등분산성 : 유의하지 않으므로 오차의 등분산성 만족 

spreadLevelPlot(fit3) #fitted vakye에 대한 스튜던트화 잔차를 시각화 

residplot(fit3)

#정규성 문제가 계속 생기기 때문에, 그 이유를 확인해보기 위해서 target 변수에 대해서 산점도를 그려보았다.

plot(student.rm$target) ; hist(student.rm$target) #완전히 정규분포를 따르지 않는다.

#왜도가 심할 경우이며, 이 경우에는 종속변수를 변환시키거나, 비모수방법인 Kruskal-Wallis test가 있다고 한다.

#interaction 교호작용까지 고려해보기 

colnames(student)

dum<-student[12:20]

#install.packages("multcomp")

#library(multcomp) #다중비교하는 패키지 

vbetahat <-vcov(fit3)

#다중비교를 위해서 전처리전의 데이터를 불러온다.

setwd("C:\\Users\\JJS810\\Desktop")

student.factor <-read.csv(list.files()[7],header=T)

str(student.factor)

student.factor<-na.omit(student.factor)

mod <-lm(target~.,data=student.factor)

mod.topic <-glht(mod,linfct=mcp(Topic="Tukey"))

summary(mod.topic) #유의한게 없다.

mod.section <- glht(mod,linfct=mcp(section="Tukey"))

summary(mod.section)

mod.gender  <- glht(mod,linfct=mcp(gender="Tukey"))

summary(mod.gender)

mod.semester <- glht(mod, linfct=mcp(semester="Tukey"))

summary(mod.semester)

mod.answer <-glht(mod, linfct=mcp(ParentAnsweringSurvey="Tukey"))

summary(mod.answer) #유의!!!!

mod.school <-glht(mod, linfct=mcp(ParentschoolSatisfaction="Tukey"))

summary(mod.school)

mod.absence <-glht(mod, linfct=mcp(StudentAbsenceDays="Tukey"))

summary(mod.absence) #유의!!!

#parentsanswer 와 studentabsence 유의 

#처음에 student에서 stepAIC로 돌린 모형으로 테스트 

summary(aov(target~.,data=student)) #교호작용은 분산분석을 통해서 알아본다. 

               StudentAbsenceDays + dumeng + dumit + dumSci+dumeng*dumit*dumSci, data=student.rm)

summary(fit.cor) #4 not defined because of singularities. 당연히 더미변수끼리는 singularity 때문에 NA로 나온다.

#probably getting that error because two or more of independent variables are perfectly collinear

colnames(student.rm)

fit.final <- lm(target ~ responsible + raisehand + VisITedResources+

                  dumeng + dumit + dumSci+ParentAnsweringSurvey*ParentAnsweringSurvey+StudentAbsenceDays*StudentAbsenceDays,data=student.rm)

summary(fit.final) #0.6356.... 결국 fit3와 같은 결과가 나온다.

#아예 처음부터 AIC로 교호작용을 고려해서 해보자.

colnames(student)

str(student)

fit.inter <- lm(target~gender+semester+responsible+raisehand*VisITedResources*AnnouncementsView*Discussion+ParentAnsweringSurvey+ParentschoolSatisfaction+StudentAbsenceDays

                +dumsecb+dumseca+dumarab+dumbio+dumeng+dumfre+dumhis+dumit+dumSci,data=student)

fit.aic.inter<-stepAIC(fit.inter,direction='both') #가변수를 제외한 numeric했던 변수끼리만 진행하였다.

#dumsecb~dumSci를 제외한 모든 변수들을 고려해보고 싶었지만, 너무 많은 경우를 고려해야되서 분석이 돌아가지 않는다.

fit.aic.inter$anova

fit.aic.step <- lm(target ~ semester + responsible + raisehand + VisITedResources + 

             AnnouncementsView + Discussion + ParentAnsweringSurvey + 

             StudentAbsenceDays + dumSci + raisehand:VisITedResources + 

             raisehand:Discussion + AnnouncementsView:Discussion, data=student) 

summary(fit.aic.step) #0.6242

Anova(fit.aic.step)

plot(fit.aic.step) #이상치 몇개를 제외하고 괜찮아 보인다. 42,46,329번

vif(fit.aic.step) #공선성 발생하기 때문에 공선성을 제거하자.

#vif문제를 해결해보자. vif가 가장 큰 변수를 순차적으로 제거한다.

fit.aic1 <-lm(target ~ semester + responsible + raisehand + VisITedResources + 

                AnnouncementsView + Discussion + ParentAnsweringSurvey + 

                StudentAbsenceDays + dumSci + raisehand:VisITedResources + 

                raisehand:Discussion ,data=student)

vif(fit.aic1)

fit.aic2 <-lm(target ~ semester + responsible + raisehand + VisITedResources + 

                AnnouncementsView + Discussion + ParentAnsweringSurvey + 

                StudentAbsenceDays + dumSci + raisehand:Discussion ,data=student)

vif(fit.aic2)

summary(fit.aic2) #0.6133으로 낮아지고, 독립변수가 유의하지 않은 것이 생겼다.

Anova(fit.aic2)

shapiro.test(residuals(fit.aic2)) #0.04715로 다시 정규성을 만족하지 않는다.

qqPlot(fit.aic2,labels=row.names(student),id.method="identify",simulate=TRUE,main="Q-Q PLOT") 

#그럼 위에서 유의하지 않은 독립변수들을 제외하고 다시 한번해보자.

#이렇게 하니까 정규성도 만족하고, R스퀘어도 괜찮고, 적당해보인다. 

fit.full <- lm(target~.,data=student.rm)

#두 모델의 F-test에서 p값이 0.9028로 유의하지 않으므로,

#fit.full에서 변수를 제거한 모형은 올바르다. 

####결론

#즉, target에 대해서, responsible이 적고, studentAbsencedays가 적을 수록 target이 올라가며,

영어, IT, 과학 부분을 학생이 배울 수록 target이 올라감을 알 수 있습니다. (크 역시 영어와 과학, IT계열이 중요하네요!!)

나머지 변수는 coefficient가 적어서, 굳히 target과 연관지어 해석하기 애매한 변수로 취급해도 될 것 같습니다.

#번외 : 아예 가능한 모든 회귀를 진행해보자. 

#all possible regression

library(leaps) #all possible regression

fit.all <-regsubsets(target~.,data=student)

summary(fit.all)

plot(fit.all,scale="adjr2") #adjusted R^2

subsets(fit.all, statistic="cp", main="Cp plot for all possible regression")

abline(1,1, lty=3,col="blue")

colnames(student)

fit.leaps <-lm(target~responsible+raisehand+VisITedResources+ParentAnsweringSurvey+StudentAbsenceDays+dumSci,data=student)

summary(fit.leaps)

plot(fit.leaps) #qqline에서 126,42,46,13번째 이상치

#bic를 기준으로 모델 생성

rs <-summary(regsubsets(target~.,data=student,nbest=ncol(student))) #nbest : number of subsets

bic.best <- which(rs$bic==min(rs$bic))

bic.vars <- names(rs$which[bic.best,])[rs$which[bic.best,]][-1] #intercept를 제외

bic.model <- paste(bic.vars, collapse="+")

bic.model <- as.formula(paste("target~",bic.model)) #식을 도출

fit.bic <- lm(bic.model, data=student.rm)

summary(fit.bic) #0.6284

#aic를 기준으로 모델 생성

n <- nrow(student)

p <- apply(rs$which,1,sum)

aic <- rs$bic-log(n)*p+2*p

aic.best <- which(aic==min(aic))

aic.vars <- names(rs$which[aic.best,])[rs$which[aic.best,]][-1]

aic.model <-paste(aic.vars,collapse = "+")

aic.model <- as.formula(paste("target~",aic.model))

fit.aic <-lm(aic.model, data=student.rm)

summary(fit.aic) #0.6345

#adjusted R^2

rs1 <- regsubsets(target~.,data=student, nbest=1, nvmax=NULL, force.in = NULL, force.out = NULL)

summary.out <- summary(rs1)

as.data.frame(summary.out$outmat)

plot(rs1, scale="adjr2", main="adjusted R^2")

which.max(summary.out$adjr2) #10

adjr.best<-names(summary.out$which[10,])[summary.out$which[10,]][-1]

adjr.model <-paste(adjr.best,collapse="+")

adjr.model <-as.formula(paste("target~",adjr.model))

fit.adjr <-lm(adjr.model, data=student.rm)

summary(fit.adjr) #0.6364

vif(fit.adjr)

shapiro.test(residuals(fit.adjr))

qqnorm(residuals(fit.adjr))

#Cp

plot(rs1, scale="Cp", main="Mellow's Cp")

which.min(summary.out$cp) #8

cp.best<-names(summary.out$which[8,])[summary.out$which[8,]][-1]

cp.model <-paste(cp.best,collapse="+")

cp.model <-as.formula(paste("target~",cp.model))

fit.cp <-lm(cp.model, data=student.rm)

summary(fit.cp) #0.6345