options(digits = 4) HSGPA <- c(3.0, 3.2, 2.8, 2.5, 3.2, 3.8, 3.9, 3.8, 3.5, 3.1) FGPA <- c(2.8, 3.0, 2.8, 2.2, 3.3, 3.3, 3.5, 3.7, 3.4, 2.9) SATV <- c(500, 550, 450, 400, 600, 650, 700, 550, 650, 550) scholar <- data.frame(FGPA, HSGPA, SATV) # collect into a data frame # install.packages("psych") # library(psych) describe(scholar) # provides descrptive information about each variable corrs <- cor(scholar) # find the correlations and set them into an object called 'corrs' corrs # print corrs pairs(scholar) # pairwise scatterplots attach(scholar) m1 <- lm(FGPA ~ SATV, data = scholar) summary(m1)
> options(digits = 4) > HSGPA <- c(3.0, 3.2, 2.8, 2.5, 3.2, 3.8, 3.9, 3.8, 3.5, 3.1) > FGPA <- c(2.8, 3.0, 2.8, 2.2, 3.3, 3.3, 3.5, 3.7, 3.4, 2.9) > SATV <- c(500, 550, 450, 400, 600, 650, 700, 550, 650, 550) > > scholar <- data.frame(FGPA, HSGPA, SATV) # collect into a data frame > describe(scholar) # provides descrptive information about each variable Error in describe(scholar) : could not find function "describe" > > corrs <- cor(scholar) # find the correlations and set them into an object called 'corrs' > corrs # print corrs FGPA HSGPA SATV FGPA 1.0000 0.9226 0.8144 HSGPA 0.9226 1.0000 0.8745 SATV 0.8144 0.8745 1.0000 > > pairs(scholar) # pairwise scatterplots > attach(scholar) The following objects are masked _by_ .GlobalEnv: FGPA, HSGPA, SATV > m1 <- lm(FGPA ~ SATV, data = scholar) > summary(m1) Call: lm(formula = FGPA ~ SATV, data = scholar) Residuals: Min 1Q Median 3Q Max -0.2804 -0.1305 -0.0566 0.0350 0.6481 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.95633 0.54419 1.76 0.1169 SATV 0.00381 0.00096 3.97 0.0041 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.27 on 8 degrees of freedom Multiple R-squared: 0.663, Adjusted R-squared: 0.621 F-statistic: 15.8 on 1 and 8 DF, p-value: 0.00412
위 linear modeling 에서 아래와 같은 선형모델을 얻는다.
y hat (FGPA) = 0.95633 + 0.00381 * SATV
p <- data.frame(c(400,500,600,700,800)) colnames(p) <- "SATV" pred <- predict(m1, p) pred[2]-pred[1] pred[5]-pred[4]
> p <- data.frame(c(400,500,600,700,800)) > colnames(p) <- "SATV" > pred <- predict(m1, p) > pred[2]-pred[1] 2 0.381 > pred[5]-pred[4] 5 0.381
아래는 원 데이터를 d 로 바꾼후 독립변인인 SATV를 평균점수를 0으로 하는 SATV_centered 로 바꾸어서 다시 regression을 한 것이다. m1_centered의 결과를 보면
d <- scholar d$SATV_centered <- d$SATV - mean(d$SATV) mean(d$SATV) m1_centered <- lm(FGPA ~ SATV_centered, data=d) summary(m1_centered)
> d <- scholar > d$SATV_centered <- d$SATV - mean(d$SATV) > mean(d$SATV) [1] 560 > m1_centered <- lm(FGPA ~ SATV_centered, data=d) > summary(m1_centered) Call: lm(formula = FGPA ~ SATV_centered, data = d) Residuals: Min 1Q Median 3Q Max -0.2804 -0.1305 -0.0566 0.0350 0.6481 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.09000 0.08531 36.22 3.7e-10 *** SATV_centered 0.00381 0.00096 3.97 0.0041 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.27 on 8 degrees of freedom Multiple R-squared: 0.663, Adjusted R-squared: 0.621 F-statistic: 15.8 on 1 and 8 DF, p-value: 0.00412
SATV의 기울기 계수인 0.00381은 모델 직선의 (데이터를 관통하는) 위치를 알려주는 것이고 이 선을 중심으로 데이터들이 포진하게 된다. 그리고 선에서 데이터 까지의 거리가 에러 값이다. 또한 이 거리는 표준편차 값을 갖게 되고 이를 이용하여 표준오차 (standard error) 값을 구해볼 수 있다. 즉 이 표준오차 값은 에러값 들이 선을 중심으로 얼마나 잘 포진되어 있는 지를 보여주는 지표가 된다.
이 논리도 연구자는 계수 값을 (b에 해당하는 계수) 표준오차 값으로 (standard error)값으로 나누고 이를 t 계산 값으로 (t-calculated value) 삼아서 significance 테스트를 할 수 있다. 표준오차가 작은 경우는 선을 중심으로 실제 데이터 값들이 좁게 모여 있음을 의미하므로 높은 t 값을 얻을 수 있겠다. 따라서 계수의 역할에 대한 통계학적인 의미가 있다고 판단한다.
> m1 <- lm(FGPA ~ SATV, data = scholar) > summary(m1) Call: lm(formula = FGPA ~ SATV, data = scholar) Residuals: Min 1Q Median 3Q Max -0.2804 -0.1305 -0.0566 0.0350 0.6481 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.95633 0.54419 1.76 0.1169 SATV 0.00381 0.00096 3.97 0.0041 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.27 on 8 degrees of freedom Multiple R-squared: 0.663, Adjusted R-squared: 0.621 F-statistic: 15.8 on 1 and 8 DF, p-value: 0.00412
위에서 SATV의
pt(3.97, df = 8, lower.tail = F)
를 이용하여 probability를 구한다. 0.0041
혹은 0.00412
이다. b1은 기울기이다 (regression coefficient). xi 지점에서 y값을 살펴보면 실제 y값에서 (y value) 기울기 선에 위치한 y값을 (y fitted value 혹은 predicted value of y) 뺀 값은 우리가 이야기한 error에 혹은 residual에 해당하는 값이다. 이 에러들이 사선을 중심으로 어떻게 포진해 있는가를 보기 위해서 b1 기울기 주변에 모인 residual 값들의 standard error 값을 알아보려면: residual의 분산값을 x의 Sum of Square 값에 제곱근을 씌운 값으로 나누어 준다. 즉,
\begin{eqnarray*} \displaystyle s_{b_{1}} & = & \sqrt {\frac {MSE}{SS_{X}}} \\ & = & \displaystyle \sqrt { \frac{1}{n-2} * \frac{SSE}{SS_{X}}} \\ & = & \displaystyle \sqrt { \frac{1}{n-2} * \frac{\Sigma{(Y-\hat{Y})^2}}{\Sigma{(X_{i} - \bar{X})^2}}} \\ \end{eqnarray*}
> sum(resid(m1)^2) [1] 0.5822 > sse <- sum(resid(m1)^2) > ssx <- sum((SATV - mean(SATV))^2) > sqrt((1/8)*(sse/ssx)) [1] 0.0009598 >
혼란스러운 것을 정리하기 위해서:
library(visreg) visreg(m1)
# load necessary libraries library(ggpubr) # create plot with regression line, regression equation and R^2 ggplot(scholar, aes(x = SATV, y = FGPA)) + geom_smooth(method = "lm") + xlim(380, 720) + geom_point() + theme_minimal()
ggplot(scholar) + aes(x = SATV, y = FGPA, size = HSGPA) + geom_point() + scale_color_gradient() + labs( y = "First Year GPA", x = "SAT score", size = "Hi School GPA" ) + theme_minimal()
m2 <- lm(FGPA~SATV+HSGPA, data=scholar) summary(m2)
> m2 <- lm(FGPA~SATV+HSGPA, data=scholar) > summary(m2) Call: lm(formula = FGPA ~ SATV + HSGPA, data = scholar) Residuals: Min 1Q Median 3Q Max -0.2431 -0.1125 -0.0286 0.1269 0.2716 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.233102 0.456379 0.51 0.625 SATV 0.000151 0.001405 0.11 0.917 HSGPA 0.845192 0.283816 2.98 0.021 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.192 on 7 degrees of freedom Multiple R-squared: 0.851, Adjusted R-squared: 0.809 F-statistic: 20.1 on 2 and 7 DF, p-value: 0.00126 > >
아래 회귀식에서 (regression equation) b 값들이 통계학적으로 유의미한지 살펴보는 방법
Y hat (FGPA) = 0.233102 + 0.000151 SATV + 0.845192 HSGPA
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.233102 0.456379 0.51 0.625 SATV 0.000151 0.001405 0.11 0.917 HSGPA 0.845192 0.283816 2.98 0.021 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
위 분석을 받아들이기 망서리게 되는 이유는 simple regression에서 SATV의 FGPA에 대한 설명력이 유의미하다는 것을 밝혔음에도 불구하고 두개의 독립변인을 동시에 고려할 때에는 중요하지 않은 변인이 되기 때문이다. 이에 대한 설명을 하기 위해서 partial and semipartial correlation 문서를 살펴본다.