다중 회귀 분석과 독립 변수의 상관

모집단 : 

 

3차원 그래프(회귀평면  )

 

 

표본 크기 50의 표본에 의한 회귀평면 추정(독립변수 간의 상관관계가 없을 때)

표본의  와 추정된 회귀평면

 

 

표본의  과  에 큰 상관이 없다면 추정된 회귀평면은 상당히 정확합니다.

 

표본 크기 50의 표본에 의한 회귀평면 추정(독립변수 간에 높은 상관이 있을 때)

표본의  와 추정된 회귀평면

 

 과  의 상관이 커짐에 따라 추정된 회귀평면의 정확성이 매우 떨어짐을 알 수 있습니다.

 

그리고 이렇게 표본 회귀계수의 오차가 커지는 것은 1차 회귀식의 추정에서 독립변수의 범위가 제한되는 것과 비슷한 이유임을 그래프에서 알 수가 있습니다.  과  가 큰 상관을 보일 경우에 특정한 방향에서 자료를 보면 모든 자료가 한정된 범위 안에 분포되어 있습니다.

 

아래의 R 코드는

모집단에서 그리고 크기 50의 표본에서  과  의 상관이

추정되는 회귀평면(특히 신뢰구간)에 어떤 영향을 주는지 보여 줍니다.

다중 회귀분석에서 Collinearity 혹은 Multi-Collinearity를 확인하기 위해 사용되는

VIF도 보여줍니다.

 

3차원 그래프의 경우,

마우스의 왼쪽 버튼을 누르고 마우스를 움직이면 시점이 변화하고,

마우스의 오른쪽 버튼을 누르고 마우스를 움직이면 확대/축소가 됩니다.

 

우선 “car”와 “rgl” 패키지를 인스톨해야 하며,

“#=====” 단위로 끊어서 실행시키면 됩니다.

 

Multiple Regression and Corr b2 independent var.R

#====== Function Def and reading library =====

predictFunc1 <- function(x1, x2) {

 

predict(mylm,newdata=data.frame(x1, x2),interval = c("confidence"),

level = 0.95,type="response")[,2] }

 

predictFunc2 <- function(x1, x2) {

predict(mylm,newdata=data.frame(x1, x2),interval = c("confidence"),

level = 0.95,type="response")[,3] }

 

library(car)

library(rgl)

#=====================================

 

#==================================

open3d(); n <- 1000;

x1<-rnorm(n, mean=55, sd=10)

x2<-rnorm(n, mean=55, sd=10)

 

x11(title="plot of x1 and x2"); plot(x1,x2)

 

y <- 1*x1 + -1*x2 + rnorm(n, sd=10); mylm <- lm(y~x1+x2)

 

plot3d(x1,x2,y, type="p", col="red", xlab="X1", ylab="X2", zlab="Y", site=5, lwd=15)

 

ranges <- rgl:::.getRanges()

x <- seq(ranges$xlim[1], ranges$xlim[2], length=9)

y <- seq(ranges$ylim[1], ranges$ylim[2], length=9)

 

planes3d(a=coef(mylm)[2], b=coef(mylm)[3], c=-1, d=coef(mylm)[1], color="red")

 

#dfbetaPlots(mylm)

print("Correlation and VIFs for population")

cor(x1, x2)

vif(mylm)

#=================

 

#====================================

n <- 50;

x1<-rnorm(n, mean=55, sd=10)

x2<-rnorm(n, mean=55, sd=10)

 

plot(x1,x2)

 

y <- 1*x1 + -1*x2 + rnorm(n, sd=10); mylm <- lm(y~x1+x2)

 

plot3d(x1,x2,y, type="p", col="red", xlab="X1", ylab="X2", zlab="Y", site=5, lwd=15)

 

ranges <- rgl:::.getRanges()

x <- seq(ranges$xlim[1], ranges$xlim[2], length=9)

y <- seq(ranges$ylim[1], ranges$ylim[2], length=9)

 

planes3d(a=coef(mylm)[2], b=coef(mylm)[3], c=-1, d=coef(mylm)[1], color="red", title="Sample 1")

 

z <- outer(x,y,predictFunc1); surface3d(x, y, z, alpha=.2)

z <- outer(x,y,predictFunc2); surface3d(x, y, z, alpha=.2)

 

#dfbetaPlots(mylm)

#dfbetaPlots(mylm)

print("Correlation and VIFs for sample 1")

cor(x1, x2)

vif(mylm)

#=================

 

#====================================

n <- 50;

x1<-rnorm(n, mean=55, sd=10)

x2<-0.86*x1+rnorm(n, 0, 5)

 

plot(x1,x2)

 

y <- 1*x1 + -1*x2 + rnorm(n, sd=10); mylm <- lm(y~x1+x2)

 

plot3d(x1,x2,y, type="p", col="red", xlab="X1", ylab="X2", zlab="Y", site=5, lwd=15)

 

ranges <- rgl:::.getRanges()

x <- seq(ranges$xlim[1], ranges$xlim[2], length=9)

y <- seq(ranges$ylim[1], ranges$ylim[2], length=9)

 

planes3d(a=coef(mylm)[2], b=coef(mylm)[3], c=-1, d=coef(mylm)[1], color="red")

 

z <- outer(x,y,predictFunc1); surface3d(x, y, z, alpha=.2)

z <- outer(x,y,predictFunc2); surface3d(x, y, z, alpha=.2)

 

#dfbetaPlots(mylm)

#dfbetaPlots(mylm)

print("Correlation and VIFs for sample 1")

cor(x1, x2)

vif(mylm)

#=================

 

#====================================

n <- 50;

x1<-rnorm(n, mean=55, sd=10)

x2<-0.91*x1+rnorm(n, 0, 3)

 

plot(x1,x2)

 

y <- 1*x1 + -1*x2 + rnorm(n, sd=10); mylm <- lm(y~x1+x2)

 

plot3d(x1,x2,y, type="p", col="red", xlab="X1", ylab="X2", zlab="Y", site=5, lwd=15)

 

ranges <- rgl:::.getRanges()

x <- seq(ranges$xlim[1], ranges$xlim[2], length=9)

y <- seq(ranges$ylim[1], ranges$ylim[2], length=9)

 

planes3d(a=coef(mylm)[2], b=coef(mylm)[3], c=-1, d=coef(mylm)[1], color="red")

 

z <- outer(x,y,predictFunc1); surface3d(x, y, z, alpha=.2)

z <- outer(x,y,predictFunc2); surface3d(x, y, z, alpha=.2)

 

#dfbetaPlots(mylm)

#dfbetaPlots(mylm)

print("Correlation and VIFs for sample 1")

cor(x1, x2)

vif(mylm)

#=================