###### 단순선형회귀 ########
women
plot(weight~height, data=women)
fit <- lm(weight~height, data=women)
abline(fit, col="blue")
summary(fit) #-87.51667은 y절편이라고함.
cor.test(women$height, women$weight) # 상관계수 구하는 함수
cor.test(women$weight, women$height)
0.9954948^2

plot(fit)
par(mfrow=c(2,2)) 
plot(fit)
par(mfrow=c(1,1))

### 다항회귀(polynomial regression)
fit2 <- lm(weight~ height + I(height^2) , data=women)
summary(fit2)
lines(women$height, fitted(fit2), col="red")
plot(fit2)

fit3 <- lm(weight~height + I(height^2) + I(height^3), data = women)
plot(fit3)

fit4 <- lm(weight~ height + I(height^2), data=women[-c(13,15),])
plot(fit4)

########  다중회귀분석 ##########
state.x77
View(state.x77)
states <- as.data.frame(state.x77[ ,  c("Murder","Population","Illiteracy","Income","Frost")])
states
View(states)
fit = lm(Murder~Population +Illiteracy + Income + Frost , data=states)
plot(fit)
summary(fit)

#install.packages("car", dependencies = TRUE)
library(car)
vif(fit)
sqrt(vif(fit))

###### 이상관측치 #######
influencePlot(fit, id.method = "identify")

states["Nevada",   ]
fitted(fit)["Nevada"]
residuals(fit)["Nevada"]

##### 회귀모형의 교정 ######
states
summary(powerTransform(states$Murder))
ncvTest(fit)
spreadLevelPlot(fit)
boxTidwell(Murder~ Population + Illiteracy , data= states)

##### 예측 변수 선택 #####
fit1 <- lm(Murder~ ., data=states)
summary(fit1)
fit2 <- lm(Murder ~ Population + Illiteracy , data=states)
summary(fit2)

#AIC (Akaike's An information Criterion)
AIC(fit1, fit2)

###### stepwise regression (Backward stepwise regression, Forward stepwise regression)
#Backward stepwise regression
full.model = lm(Murder~. , data = states)
reduced.model = step(full.model, direction = "backward")
summary(reduced.model)

#Forward stepwise regression
min.model = lm(Murder~1, data = states)
fwd.model  <- step(min.model, direction="forward", scope = (Murder~Population+ Illiteracy + Income
                                                            + Frost), trace=0)
summary(fwd.model)

## all subset regression
library(leaps)
leaps <- regsubsets(Murder~ Population + Illiteracy + Income + Frost, data= states, nbest=4)
plot(leaps, scale="adjr2")


########## Logistic Regression ########
require(survival)
str(colon)
colon1 <- na.omit(colon)
View(colon)
View(colon1)
result <- glm(status ~ rx+sex+age+obstruct+perfor + adhere + nodes + differ + extent + surg, family = binomial, data=colon1)
summary(result)

reduced.model = step(result)
summary(reduced.model)

require(moonBook)
extractOR(reduced.model)
fit = glm(formula = status ~ rx + obstruct + adhere + nodes + extent + surg, family = binomial, data = colon1)
fit.od = glm(formula = status ~ rx + obstruct + adhere + nodes + extent + surg, family = quasibinomial, data = colon1) 
pchisq(summary(fit.od)$dispersion*fit$df.residual, fit$df.residual, lower = F)
#0.2803691이 값이 0.05보다 크다면 과산포는 없다고 확신할 수 있습니다.

?ORplot()
plot()
ORplot(fit, main = "Plot for Odds Ratios")
ORplot(fit, type=2, show.OR=FALSE, show.CI=TRUE, pch=15, lwd=2, col=c("darkblue", "red"), main="Plot of OR" )
ORplot(fit, type=3, show.OR=FALSE, show.CI=TRUE, pch=15, lwd=2, col=c("darkblue", "red"), main="Plot of OR")


##### Poisson Regression #######

install.packages("robust")
library(robust)
data(breslow.dat, package = "robust")
summary(breslow.dat)

install.packages("qcc")
library(qcc)
qcc.overdispersion.test(breslow.dat$sumY, type="poisson")

fit = glm(sumY ~ Base + Age + Trt, family = quasipoisson, data= breslow.dat)
summary(fit)

install.packages("moonBook")
library(moonBook)

extractOR(fit)
extractOR(fit, digits = 3)
ORplot(fit, type = 2, show.CI=TRUE, main="Plot for Quasipoisson")


#### Survival Analysis #####

require(survival)
data("colon")
View(colon)
colon <- na.omit(colon)
rm(colon1)
str(colon)

colon$TS <- Surv(colon$time,colon$status==1)
fit = survfit(TS ~ rx, data = colon)
plot(fit)
plot(fit,col=1:3, lty=1:3)
legend("topright", legend=levels(colon$rx), col=1:3, lty = 1:3)


###### Cumulative hazard Start #####
plot(fit, col=1:3, lty=1:3, fun ="cumhaz", mark.time=FALSE, ylab="Cumulative hazard")
legend("topleft", legend = levels(colon$rx), col=1:3, lty =1:3)
###### Cumulative hazard End #######


#### Log-rank test ####
survdiff(Surv(time, status ==1)~rx, data=colon)

#### Cox Regression #####
out = coxph(Surv(time, status ==1)~rx, data=colon)
summary(out)

######  Hazard ratios of all individual variables #####
colon$TS <- Surv(colon$time,colon$status==1)
out = coxph(colon$TS~rx, data=colon )
#install.packages("moonBook")
require(moonBook)
attach(colon)
out = mycph(TS~.-id-study-time-status-etype, data=colon)
out
out2 =coxph(TS~. -id-study-time-status-etype, data=colon)
final =step(out2, direction = "backward")
HRplot(out, type =2, show.CI = TRUE)