한림대학교 동탄성심병원 통계워크샵-회귀분석

한림대학교 동탄성심병원 2차 통계워크샵 - 회귀분석

L_dontan_reg.R

#################################################################

###### 단순선형회귀 ########

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)