고급회귀분석 실습
chapter 3, chapter 4
library(ggplot2)
dt <- data.frame(x = c(4,8,9,8,8,12,6,10,6,9),
y = c(9,20,22,15,17,30,18,25,10,20))
dt
correlation check
cor(dt$x, dt$y)
산점도 확인
plot(y~x,
data = dt,
xlab = "광고료",
ylab = "총판매액",
pch = 16,
cex = 2,
col = "darkorange")
- pch 점 모양
-
cex 점 크기
-
양의상관관계 강하네,
- 우상향이네, 단순상관선형 적용해보면 되겠다.
ggplot(dt, aes(x, y)) +
geom_point(col='steelblue', lwd=2) +
# geom_abline(intercept = co[1], slope = co[2], col='darkorange', lwd=1.2) +
xlab("광고료")+ylab("총판매액")+
# scale_x_continuous(breaks = seq(1,10))+
theme_bw() +
theme(axis.title = element_text(size = 14))
$\hat{ y} = \widehat{(E(y|X=x))} = \hat{\beta_0} + \hat{\beta_1} * x$
$H_0$ : $\beta_0$ =0 vs $H_1$ : $\beta_0 \ne 0$
$H_0$ : $\beta_1 =0$ vs $H_1$ : $\beta_1 \ne 0$
모형 적합을 한다 yhat을 찾는다.
- 회귀분석을 한다. 평균 반응을 추정한다.
lm linear model 사용
model1 <- lm(y ~ x, dt)
# lm(y ~ 0 + x, dt) beta0 없이 분석하고 싶을때
model1
설명변수 x 하나일때
- beta0hat = -2.270
- beta1hat = 2.609
summary(model1)
- 모형의 유의성 검정 자체(f검정)
- 개별 회귀계수에 대한 유의성검정(t검정)
- beta1은 유의하지 않다
- beta0는 유의하다
- f통계량은 45.24(msr/mse) p값 충분히 작아서 모형은 유의하다.
- y의 총 변동 중에 85%정도를 설명하고 있다.
- root mse(RMSE) = 2.631
6.726**2
단순선형에서만 해당
0.000149도 같음
names(model1)
model1$residuals # 보고 싶은 변수 입력해봐~
model1$fitted.values ##hat y
model1$coefficients
anova(model1) ## 회귀모형의 유의성 검정
- 설명변수의 개수가 x 자유도
- 잔차의 자유도는 n-2
a <- summary(model1)
ls(a)
summary(model1)$coef ## 회귀계수의 유의성 검정
confint(model1, level = 0.95) ##회귀계수의 신뢰구간
## beta +- t_alpha/2 (n-2) * se(beta)
qt(0.025, 8)
qt(0.975, 8)
- qt _ tquantile
model2 <- lm(y ~ 0 + x, dt)
summary(model2)
- intercept 없는 모습
- r squre가 두 번째가 높고,
- p값도 훨씬 유의하게 나옴
anova(model1)
anova(model2)
plot(y~x, data = dt,
xlab = "광고료",
ylab = "총판매액",
pch = 20,
cex = 2,
col = "darkorange")
abline(model1, col='steelblue', lwd=2)
abline(model2, col='green', lwd=2)
model들이 기울기가 살짝씩 다르다
co <- coef(model1)
ggplot(dt, aes(x, y)) +
geom_point(col='steelblue', lwd=1) +
geom_abline(intercept = co[1], slope = co[2], col='darkorange', lwd=1) +
xlab("광고료")+ylab("총판매액")+
theme_bw()+
theme(axis.title = element_text(size = 16))
# lm을 사용하지 않고 구할때
dt1 <- data.frame(
i = 1:nrow(dt),
x = dt$x,
y = dt$y,
x_barx = dt$x - mean(dt$x), # x - x평균
y_bary = dt$y - mean(dt$y)) # y - y평균
dt1$x_barx2 <- dt1$x_barx^2 # x 편차의 제곱
dt1$y_bary2 <- dt1$y_bary^2 # y편차의 제곱
dt1$xy <-dt1$x_barx * dt1$y_bary
dt1
round(colSums(dt1),3)
##hat beta0 = bar y - hat beta_1 * bar x
beta1 <- as.numeric(colSums(dt1)[8]/colSums(dt1)[6])
beta0 <- mean(dt$y) - beta1 * mean(dt$x)
cat("hat beta0 = ", beta0)
cat("hat beta1 = ", beta1)
구분할 수 있어야 한다
신뢰구간 달라진다.
## y = E(Y|x0) + epsilon 개별 y 추정
# x0 = 4.5
new_dt <- data.frame(x = 4.5)
predict(model1,
newdata = new_dt,
interval = c("confidence"), level = 0.95)
new_data=new_df 정의 안 하면 fitted value가 나온다.
confidence는 평균반응
predict(model1, newdata = new_dt,
interval = c("prediction"), level = 0.95)
prediction은 개별 y 추정
신뢰구간이 커진다. $\to$ 표준오차가 달라지기 때문
dt_pred <- data.frame(
x = 1:12,
predict(model1,
newdata=data.frame(x=1:12),
interval="confidence", level = 0.95))
dt_pred
dt_pred2 <- as.data.frame(predict(model1,
newdata=data.frame(x=1:12),
interval="prediction", level = 0.95))
dt_pred2
names(dt_pred2)[2:3] <- c('plwr', 'pupr')
plot 같이 그리게 데이터 합치기
dt_pred3 <- cbind.data.frame(dt_pred, dt_pred2[,2:3])
barx <- mean(dt$x)
bary <- mean(dt$y)
plot(y~x, data = dt,
xlab = "광고료",
ylab = "총판매액",
pch = 20,
cex = 2,
col = "grey",
ylim = c(min(dt_pred3$plwr), max(dt_pred3$pupr)))
abline(model1, lwd = 5, col = "darkorange")
lines(dt_pred3$x, dt_pred3$lwr, col = "dodgerblue", lwd = 3, lty = 2)
lines(dt_pred3$x, dt_pred3$upr, col = "dodgerblue", lwd = 3, lty = 2)
lines(dt_pred3$x, dt_pred3$plwr, col = "dodgerblue", lwd = 3, lty = 3)
lines(dt_pred3$x, dt_pred3$pupr, col = "dodgerblue", lwd = 3, lty = 3)
abline(h=bary,v=barx, lty=2, lwd=0.2, col='dark grey')
ggplot(dt_pred3, aes(x, fit)) +
geom_line(col='steelblue', lwd=2) +
xlab("")+ylab("")+
scale_x_continuous(breaks = seq(1,10))+
geom_line(aes(x, lwr), lty=2, lwd=1.5, col='darkorange') +
geom_line(aes(x, upr), lty=2, lwd=1.5, col='darkorange') +
geom_line(aes(x, plwr), lty=2, lwd=1.5, col='dodgerblue') +
geom_line(aes(x, pupr), lty=2, lwd=1.5, col='dodgerblue') +
geom_vline(xintercept = barx, lty=2, lwd=0.2, col='dark grey')+
geom_hline(yintercept = bary, lty=2, lwd=0.2, col='dark grey')+
theme_bw()
bb <- summary(model1)$sigma * ( 1 + 1/10 +(dt$x - 8)^2/46)
dt$ma95y <- model1$fitted + 2.306*bb
dt$mi95y <- model1$fitted - 2.306*bb
ggplot(dt, aes(x=x, y=y)) +
geom_point() +
geom_smooth(method=lm , color="red", fill="#69b3a2", se=TRUE) +
geom_line(aes(x, mi95y), col = 'darkgrey', lty=2) +
geom_line(aes(x, ma95y), col = 'darkgrey', lty=2) +
theme_bw() +
theme(axis.title = element_blank())
dt
dt$yhat <- model1$fitted
# fitted.values(model1) # y에 대한 추정값 구하기
dt$resid <- model1$residuals
# resid(model1)
par(mfrow=c(1,2))
plot(resid ~ x, dt, pch=16, ylab = 'Residual')
abline(h=0, lty=2, col='grey')
plot(resid ~ yhat, dt, pch=16, ylab = 'Residual')
abline(h=0, lty=2, col='grey')
par(mfrow=c(1,1))
단순선형에서는 두 plot의 차이가 없다.
- 선형성 만족
- 등분산성 나름 만족
- 정규성 아웃라이어 있는 거 같은데..
- 독립성?
library(lmtest)
dwtest(model1, alternative = "two.sided") #H0 : uncorrelated vs H1 : rho != 0
# dwtest(model1, alternative = "greater") #H0 : uncorrelated vs H1 : rho > 0
# dwtest(model1, alternative = "less") #H0 : uncorrelated vs H1 : rho < 0
p 값 커서 기각할 수 없다.
첫 번째꺼 주로 보기
qqnorm(dt$resid, pch=16)
qqline(dt$resid, col = 2)
분위수분위수 그림
- 정규분포의 실제
- 어떤 분포의 이론적 분위수와 내가 가진 sample의 분위수 비교
주로 꼬리쪽을 많이 본다.
- 이 데이터의 경우 꼬리부분이 차이가 커 보임
ggplot(dt, aes(sample = resid)) +
stat_qq() + stat_qq_line() +
theme_bw()
shapiro.test(dt$resid) ##shapiro-wilk test
#H0 : normal distributed vs H1 : not
p값 작게 나오면 정규분포라고 하기 어렵다.
- 정규성은 잔차를 넣어줬는데
- bptest는 model을 넣었다.
bptest(model1) #Breusch–Pagan test
# H0 : 등분산 vs H1 : 이분산
library(UsingR)
data(father.son)
names(father.son)
lm.fit<-lm(sheight~fheight, data=father.son)
summary(lm.fit)
아버지의 키가 아들의 키의 25%만 설명
plot(sheight~fheight,
data=father.son,
pch=16, cex=0.5,
xlab="father’s height (inches)",
ylab="son’s height (inches)")
abline(lm.fit)
amazon<-read.csv("amazon.csv")
plot(High ~Year , amazon, pch=16)
lm.fit<-lm(High~Year, data=amazon)
summary(lm.fit)
confint(lm.fit)
par(mfrow=c(1,2))
scatter.smooth(x=1:length(amazon$Year), y=residuals(lm.fit), xlab="Year")
scatter.smooth(x=predict(lm.fit), y=residuals(lm.fit), xlab=expression(hat(y)))
library(lmtest)
dwtest(lm.fit)
양의 상관관계가 있다.
- 시간순 index
- 최근 관측 데이터에 영향 많이 받는 편