AI Analytics, by YD Lee

회귀분석 Regression, Part I

Topics:

Part I:
- iris 자료
- 단순 선형회귀
Part II
- 사영과 회귀분석
- 다중 선형회귀
Part III
- 가변수와 선형모형
- 모형의 다양성
Part IV
- 모형 선택
- 회귀 진단

1. iris 자료

Subspecies of iris

setosa versicolor virginica

Parts of flower : Sepal & Petal

iris data :
- Wikipedia 설명 [iris data, wiki]
- R Help [help(iris)]

3종류의 붓꽃, setosa, versicolor, virginica 각 50개씩 전체 150개의 꽃에 대하여 꽃받침의 길이와 폭, 꽃잎이 길이와 폭을 측정한 자료

dim(iris)

## [1] 150   5

head(iris)

##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

편리한 사용을 위하여, 별도로 irix 자료를 만들어 사용한다

irix 자료는 iris 자료에서 다음과 같이 변수이름을 간략하게 만든 자료다
- sl : Sepal.Length (꽃받침의 길이)
- sw : Sepal.Width (꽃받침의 폭)
- pl : Petal.Length (꽃잎의 길이)
- pw : Petal.Width (꽃잎의 폭)
- sp : Species (붓꽃의 종류)

irix <- iris
names(irix) <- c('sl','sw','pl','pw','sp')
irix[c(1,50,51,100,101,150),]

##      sl  sw  pl  pw         sp
## 1   5.1 3.5 1.4 0.2     setosa
## 50  5.0 3.3 1.4 0.2     setosa
## 51  7.0 3.2 4.7 1.4 versicolor
## 100 5.7 2.8 4.1 1.3 versicolor
## 101 6.3 3.3 6.0 2.5  virginica
## 150 5.9 3.0 5.1 1.8  virginica

library(GGally,warn.conflicts = FALSE, quietly=TRUE)

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

ggpairs(irix,aes(col=sp),upper='blank', progress=FALSE)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

setosa 자료

iris 자료 중, setosa 50개에 대한 자료만을 뽑아서 setosa 데이터를 구성하자

setosa 에서, 꽃받침의 폭 sw 와 꽃받침의 길이 sl 사이의 관계를 살펴보자.

아래 그림에서 50 개의 점들은 sw 가 커질수록 sl 값도 커지는 관계가 있음을 보여준다.

파란선은 50개의 점들을 대표하는 선으로, sw 와 sl 사의 관계를 함축한다.

library(ggplot2)
setosa <- irix[1:50,-5]  # 50 setosa only, drop 5th variable sp
p <- ggplot(setosa,aes(x=sw, y=sl)) + geom_point()
p <- p + geom_smooth(method='lm', se=FALSE)
p

## `geom_smooth()` using formula 'y ~ x'

2. 단순 선형회귀

회귀모형의 가정

관측점들과 가정:
- 관측하여 알 수 있는 것: 50개의 $(sw,sl)$ 순서쌍으로 표현되는 관측점들
- 가정되는 대상: 직접 관측할 수 없고 추상적으로 가정되는 대상
  - ‘sw 와 sl 사이의 함수적 관계’
  - ’sw 와 sl 사이의 함수적 관계’의 형태
  - ’관측점들’과 ’sw 와 sl 사이의 함수적 관계’에 대한 관계

단순 선형회귀모형에서의 가정:

일반화 하여, 변수 sl 를 $y$, 변수 sw 를 $x$ 라고 부르고, 관측값의 개수를 $n$ 이라 하자
- 모형:
  - $y_j= \alpha + \beta x_j + \epsilon_j$, $\quad j=1,\ldots,n$ 이고,
  - $\epsilon_j ~\stackrel{iid}{\sim} ~N(0, \sigma^2 ), \quad j=1,\ldots,n,~~$
  - 관측가능한 값 중 $x_j$ 는 비확률적인 값이다
- 해설:
  - 두 변수 $x$ 와 $y$ 사이에는 ’직선적 관계’가 있다. 즉, $y=\alpha +\beta \, x$ 다
    
    이 직선의 정확한 식은 알 수 없다. 즉, 절편 $\alpha$와 기울기 $\beta$ 는 모른다
  - $n$ 개의 관측값 $(x_j , y_j)$ 은 ‘직선적 관계’와 ’오차’ $\epsilon_j$ 로부터 얻어진다
  - 관측가능한 값 중 $y_j$ 가 확률적 값이 되는 이유는 오차 $\epsilon_j$ 가 확률적이기 때문이다
  - 보통, 확률변수는 대문자로 표시하고, 그 확률변수가 관측된 값을 나타내기 위하여 소문자를 사용한다.
    
    이에 따르면, $Y_j$ 는 확률변수를 의미하고 $y_j$ 는 그 관측된 값을 의미하게 된다.
    
    여기서는, 큰 혼동이 없는 범위 내에서, 이를 구별하지 않고 소문자를 사용하여 $y_j$ 라고 하기로 한다.

단순 선형회귀모형의 표현:

$x$ 가 비확률적인 값으로 가정되므로, 다음 세 가지 표현은 동일한 의미다

$\qquad \qquad \qquad y_j= \alpha + \beta x_j + \epsilon_j$, $~~ \epsilon_j ~\stackrel{iid}{\sim} ~N(0, \sigma^2 ), \quad j=1,\ldots,n$

$\qquad \qquad \qquad y_j ~\stackrel{independent}{\sim} ~N(~\alpha+\beta x_j , ~\sigma^2 ~), \quad j=1,\ldots,n$

$\qquad \qquad \qquad y_j ~\stackrel{independent}{\sim} ~N(~\mu_j , ~\sigma^2 ~), ~~ \mu_j = \alpha+\beta x_j , \quad j=1,\ldots,n$

’정규분포 모집단 모형’과의 비교:
- 정규분포 모집단 모형: $y_j ~\stackrel{iid}{\sim} ~N(~\mu, ~\sigma^2 ~), ~~ j=1,\ldots,n$
- 회귀모형에서는 모평균 $\mu$ 가 $j$ 에 따라 변화되고, 즉, $\mu_j$ 이고,
  
  그 변화의 원인이 되는 값은 $x_j$ 이고, 관측에 의하여 알 수 있다

아래 그림에서,
- 빨간선은 가정되는 $x$와 $y$ 사이의 직선적 관계를 의미하는 선으로
  
  직접 관측할 수 없고 그 정확한 값이 얼마인지 알 수 없다
  
  아래 그림의 빨간선은 실재하는 선이 아니고 가상적인 예로 그린 것이다
- 50개의 검정색 관측점들은 관측에 의하여 얻어지는 값들이다
- 파란선은 관측된 점들로부터 ’추정되는 직선적 관계’를 의미하고
  
  관측점들에 의하여 구해질 수 있는 직선이다
- 파란선과 빨간선은 그 값과 성격이 다른 직선이고 그 차이가 얼마인지 알 수 없다

res <- lm(sl~sw,setosa)
p+geom_abline(intercept=coef(res)[1]*0.9 ,slope=coef(res)[2]*1.12,col='red', linetype='dashed')

## `geom_smooth()` using formula 'y ~ x'

회귀분석 용어

변수와 모수:
- 독립변수: 설명변수, $x$ 변수, 종속변수의 변동원인을 제공하는 변수
- 종속변수: $y$ 변수, 결과로 나타나는 변수
- 오차: $\epsilon$, 확률적인 값, 직접 관측되는 값은 아님
- 모수: 회귀모형에서 추정하여야 할 모수는
  - $\alpha$: 절편, 직선의 y-절편
  - $\beta$: 기울기, 직선의 기울기
  - $\sigma^2$: 분산, 오차의 분산
- 회귀식: $E(y|x)$ 를 의미함. 가정으로부터, 직선 $y=\alpha+\beta x$ 를 의미함

추정값과 예측값:
- 모수 추정값 :
  - $\hat{\alpha}$: 절편 추정값
  - $\hat{\beta}$: 기울기 추정값
  - $\hat{\sigma^2}$: 분산 추정값
- 추정값: $y$의 관측값들 $y_j$에 대한 추정값, 즉, $\hat{y}_j = \hat{\alpha} + \hat{\beta} x_j$
- 추정오차(잔차) : $e_j = y_j - \hat{y}_j$, $~~\hat{\epsilon}_j$ 라고도 함
- 추정회귀식: $\hat{y}(x) = \widehat{E(y|x)} = \hat{\alpha}+\hat{\beta} x$
- 예측값 Predictor: $\hat{y}(x)$
  - $x$ 가 주어진 경우의, $y$에 대한 예측값, 즉 $y(x)$ 를 예측하는 값
  - 예측값 $\hat{y}(x)$ 은, 기호로는 추정회귀식과 동일하게 표현되지만
    
    예측값으로서의 $\hat{y}(x)$은, 확률적인 값 $E(y|x)+\epsilon$ 에 대한 predictor 를 의미

통계 용어 정리

estimator, predictor, forecastor:
- estimator : 추정값, 추정량, 비확률적인 값에 대하여 사용
- predictor : 예측값, 예측량, 확률적인 값에 대하여 사용
- forecastor : 선행예측값, 시간에 따라 발생하는 학률적인 값을 발생 이전 예측한 값

외생변수, 내생변수: 구조방정식모형에서 사용되는 용어
- 구조방정식모형: 독립변수 $x$ 가 확률적인 경우를 고려하는, 회귀모형이 발전된 모형
- 외생변수 exogenous variable : 독립변수, 변동원인을 외부로 내보내는 변수
- 내생변수 endogenous variable : 종속변수, 변동원인을 내부로 받아들이는 변수

회귀모형에서의 추정

관측점들과 직선 사이의 차이: $~~\epsilon_j = y_j - (\alpha + \beta x_j )$
- $\epsilon_j ~\stackrel{iid}{\sim}~ N(0, \sigma^2 ), ~~j=1,\ldots,n$
- $\sum_j^n \epsilon_j ~\sim~ N(0, n \sigma^2 )$
- $\bar{\epsilon}_n = (1/n) \sum_j^n \epsilon_j ~\sim~ N \left( 0, (1/n) \sigma^2 \right)$
- $E(\bar{\epsilon}_n )= 0$, $~~Var( \bar{\epsilon}_n) = \sigma^2/n$
- $\sum_j^n \epsilon_j^2 ~\sim~ \sigma^2 \cdot \chi^2 (n)$
- $E\left(\sum_j^n \epsilon_j^2 \right )= n \sigma^2$, $~~Var\left( \sum_j^n \epsilon_j^2 \right )= 2 n \sigma^4$

관측점들을 가장 잘 요약(설명)하는 직선이란?
- 관측점들과 직선 사이의 거리가 가장 작게 만드는 직선
- 관측점들과 직선 사이의 거리:
  - 절대값거리: $L(\alpha, \beta)=\sum_j \left| y_j -(\alpha + \beta x_j) \right|$
  - 제곱거리: $S(\alpha, \beta)=\sum_j \left( y_j -(\alpha + \beta x_j) \right)^2$
- 최소자승법: $S(\alpha, \beta)$ 를 최소로 만드는 직선 (즉, $\alpha$ 와 $\beta$)을 선택
- 최소자승 추정량: $S(\alpha, \beta)$ 를 최소로 만드는 $\alpha$ 와 $\beta$ 값
  
  이를 기호로 각각 $\hat{\alpha}$, $\hat{\beta}$ 라고 하자. 즉, $S(\hat{\alpha}, \hat{\beta}) \le S(\alpha, \beta), ~~\forall \alpha, \beta$
- 최소자승법의 장점: 일반적으로 사용하는 방법 (좋은 방법)
  - 미분이론을 적용하여, 간단하게 식을 유도할 수 있다
  - 정규분포와 관련되어 이론을 발전시킬 수 있다
- 절대값거리 최소화 방법의 장점: 특별한 경우에 사용
  - 계산을 간단히 처리할 수 있다 (LP 알고리즘 이용)
  - 이중지수분포와 잘 어울리는 방법
  - 이상치 outlier 에 덜 민감하다

최소자승(추정)법:
- $S(\alpha, \beta)$ 를 최소화 하는 값 $\hat{\alpha}$, $\hat{\beta}$ 구하기
- $S(\alpha, \beta)$ 를 최소화 하는 $\alpha$, $\beta$ 로 미분하면
  - $(\partial/\partial \alpha) S(\alpha, \beta)= -2 \, \sum_j \left( y_j -(\alpha + \beta x_j) \right)$
  - $(\partial/\partial \beta) S(\alpha, \beta) = -2 \sum_j \, x_j \left( y_j -(\alpha + \beta x_j) \right)$
- $\hat{\alpha}$ 과 $\hat{\beta}$ 은 위의 두 미분 식을 0으로 만드는 연립방정식의 근
  - $\sum_j \left( y_j -(\hat{\alpha} + \hat{\beta} x_j) \right) =0,~~$ (식1)
  - $\sum_j \, x_j \left( y_j -(\hat{\alpha} + \hat{\beta} x_j) \right)=0,~~$ (식2)
- (식1)을 정리하여 다시 표현하면, $\hat{\alpha} + \bar{x}_n \hat{\beta} = \bar{y}_n,~~$ (식3)
- (식1)에 $\bar{x}_n$ 를 곱하여 (식2)에서 빼고, (식3)을 대입하여, 다시 정리하면,
  - $\sum_j (x_j - \bar{x}_n)^2 \hat{\beta} = \sum_j (x_j -\bar{x}_n )(y_j -\bar{y}_n)$
- 즉, 이로부터 $\hat{\alpha}$ 과 $\hat{\beta}$ 은 다음과 같이 구해진다
  - $\hat{\beta} = \left[\sum_j (x_j -\bar{x}_n )(y_j -\bar{y}_n)\right]/\left[\sum_j (x_j - \bar{x}_n)^2 \right]$
  - $\hat{\alpha} = \bar{y}_n - \bar{x}_n \hat{\beta}$
최소자승값의 성질:
- 관측점들과 추정직선 사이의 차이: $~e_j = y_j - (\hat{\alpha} + \hat{\beta} x_j )$
- $S = S(\hat{\alpha}, \hat{\beta}) = \sum_j e_j^2$
  
  $S$ : 잔차제곱합, SSE, Sum of Squared Error, 적합결여도
- $S(\hat{\alpha}, \hat{\beta}) \le S(\alpha, \beta)$
- $S(\hat{\alpha}, \hat{\beta}) ~\sim~ \sigma^2 \cdot \chi^2 (n-2),~~$ (why?)
- $E(S) = E\left(\sum_j^n e_j^2 \right )= (n-2) \sigma^2$,
- $Var(S) = Var\left( \sum_j^n e_j^2 \right )= 2 (n-2) \sigma^4$

오차분산의 추정량:
- $\hat{\sigma}^2 = (\sum_j e_j^2 )/ (n-2) = S/(n-2)$
- 불편추정량 : $E(\hat{\sigma}^2 ) = \sigma^2$

단순 선형회귀 모형의 특수한 경우:
- 정규분포 모집단 모형: $~~y_j = \alpha +\epsilon_j ,~~ \epsilon_j ~\stackrel{iid}{\sim}~ N(0,\sigma^2 )$
- 원점을 지나는 단순회귀 모형: $~~y_j = \beta x_j +\epsilon_j,~~ \epsilon_j ~\stackrel{iid}{\sim}~ N(0,\sigma^2 )$

setosa 데이터에서의 단순 선형회귀 적합 결과

회귀계수 추정값

$\hat{\alpha}=2.639$
$\hat{\beta}=0.6905$

( res <- lm(sl~sw,setosa) )

## 
## Call:
## lm(formula = sl ~ sw, data = setosa)
## 
## Coefficients:
## (Intercept)           sw  
##      2.6390       0.6905

ANOVA Table

두 모형의 비교:
- Reduced Model: $y_j = \mu +\epsilon_j,~~ \mu= \alpha,~~$ (정규분포 모집단 모형)
- Full Model: $y_j = \mu_j +\epsilon_j, ~~\mu_j = \alpha+\beta x_j , ~~$(단순 선형회귀 모형)
- Reduced Model 은, Full Model 에서 $\beta$ 값이 0 으로 고정된 특수한 경우

각 모형에서의 (최소)잔차제곱합:
- 최소잔차제곱합: 모형과 관측자료와의 (최소)거리, 적합결여도
- Reduced Model : $~~S_1 = \sum_j (y_j -\hat{\mu})^2$
- Full Model: $~~S_2 = \sum_j (y_j -\hat{\mu_j})^2$
- $S_2 \le S_1,~~$ (why?)
- $SST = S_1 ~\sim~ \sigma^2 \cdot \chi^2 (n-1)$
- $SSE= S_2 ~\sim~ \sigma^2 \cdot \chi^2 (n-2)$
- $SSR= S_1 - S_2 ~\sim~ \sigma^2 \cdot \chi^2 (1)$
- $SST = SSR + SSE$
- $SSR= S_1 - S_2= \sum_j (\hat{\mu}_j- \hat{\mu})^2, ~~$ (why?)
- $SSR$ 은 두 모형에서의 적합결여도의 차이
- $SSR$ 이 작다면, 즉 두 모형의 적합효과의 차이가 크지 않다면,
  
  단순한 모형인 Reduced Model 이 더 좋은 모형이라는 의미가 됨.
  
  반대로, SSR 이 값이 커지면, Full Model 이 더 좋다는 의미

ANOVA Table : Analysis of Variance Table
- SST, SSR, SSE 값을 계산하고 자유도와 함께 제시한 표
- MSR = SSR/1, MSE=SSE/(n-2)
- F = MSR/MSE 는, $H_0$ 가 사실일 때, 자유도 (1, n-2) 인 F 분포를 따름
- 정규분포 모집단 모형($\beta=0$)과 단순 선형회귀모형($\beta \neq 0$)에 대한 비교
- 가설 : $H_0 : \beta=0 ~~ \mbox{vs.}~~ H_1 : \beta \neq 0$
- p-value 가 작을수록 $~H_1~$ 을 지지함
  
  p-value 가 alpha-value (보통은 0.05) 보다 작으면, $~H_1~$ 으로 결론을 내림

Source	Df	Sunm Sq	Mean Sq	F value	Pr(F>0)
$x$	1	SSR	MSR	F	p-value
Residuals	n-2	SSE	MSE
Total	n-1	SST

ANOVA Table for setosa

setosa 자료에 대하여, 독립변수 sw, 종속변수 sl 을 이용한 단순 선형회귀모형의 ANOVA Table
SSR = 3.3569, 자유도 1, MSR = 3.3569
SSE = 2.7313, 자유도 48, MSE = 0.0569
$\hat{\sigma}^2 = 0.0569 = MSE$ (MSE 는 오차분산 $\sigma^2$에 대한 추정량)
p-value : $6.71 \times 10^{-10}$

anova(res)

## Analysis of Variance Table
## 
## Response: sl
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## sw         1 3.3569  3.3569  58.994 6.71e-10 ***
## Residuals 48 2.7313  0.0569                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Summary Table

5 number summary of $~~e_j,~j=1,\ldots,n$

Table for Regression Coefficients :
- Estimate: $\hat{\alpha}$, $\hat{\beta}$, 모수 추정량
- Variance of the estimators:
  - $q_n = \sum_j x_j^2$, $~~~s_n = \sum_j (x_j - \bar{x}_n )^2$
  - $Var(\hat{\alpha})= (\sigma^2/n) \cdot ( q_n /s_n )$
  - $Var(\hat{\beta})= (\sigma^2/n) \cdot ( n / s_n )$
  - $\widehat{Var}(\hat{\alpha})= (\hat{\sigma}^2/n) \cdot ( q_n /s_n )$
  - $\widehat{Var}(\hat{\beta})= (\hat{\sigma}^2/n) \cdot ( n / s_n )$
  - 위와 같이 계산되는 이유는 ’다중 선형회귀모형’에서 설명함
- Std. Error (of Estimated Coefficients) : 추정량 표본오차
  - $se(\hat{\alpha}) = \sqrt{\widehat{Var}(\hat{\alpha})}$
  - $se(\hat{\beta}) = \sqrt{\widehat{Var}(\hat{\beta})}$
- t-value: 통계량(t-value)은 자유도 $n-2$ 인 t-분포를 따름
  - 통계량: $t = \hat{\alpha}/se(\hat{\alpha})$
    
    가설: $H_0 : \alpha =0, ~~\mbox{vs.}~~ H_1: \alpha \neq 0$
  - 통계량: $t = \hat{\beta}/se(\hat{\beta})$
    
    가설: $H_0 : \beta =0, ~~\mbox{vs.}~~ H_1: \beta \neq 0$
  - $t$ 의 절대값, $|t|$ 가 클수록, $H_1$ 을 지지함
  - p-value: Pr(>|t|) 로 표시됨. 작을수록 (보통은 0.05 이하), $H_1$ 을 지지함

Residual standard error : $\hat{\sigma} =\sqrt{MSE}$

Multiple R-squared: $R^2 = SSR/SST = 1- SSE/SST$
Adjusted R-squared: $Adj. R^2 = 1- MSE/MST$, $~~ MST=SST/(n-1)$

Summary Table for setosa

Estimate: $\hat{\alpha}$=2.639, $\hat{\beta}$=0.6905
$se(\hat{\alpha})$ = 0.31, $se(\hat{\beta})$ = 0.0899
$T = \hat{\alpha}/se(\hat{\alpha})$ = 0.8513

$T = \hat{\beta}/se(\hat{\beta})$ =0.7681

$\hat{\sigma}$ = 0.2385,
$R^2$ = 0.5514, $Adj. R^2$= 0.542

summary(res)

## 
## Call:
## lm(formula = sl ~ sw, data = setosa)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.52476 -0.16286  0.02166  0.13833  0.44428 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.6390     0.3100   8.513 3.74e-11 ***
## sw            0.6905     0.0899   7.681 6.71e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2385 on 48 degrees of freedom
## Multiple R-squared:  0.5514, Adjusted R-squared:  0.542 
## F-statistic: 58.99 on 1 and 48 DF,  p-value: 6.71e-10

R 에서 lm 결과의 출력물

추정잔차의 자유도: df.residual, res$df.residual
회귀계수: res$coefficients, coeff(res)
추정잔차: res$residuals, resid(res)
추정값 : res$fitted.values, predict(res)

attributes(res)

## $names
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"        
## 
## $class
## [1] "lm"

res$df.residual

## [1] 48

res$coefficients

## (Intercept)          sw 
##   2.6390012   0.6904897

coef(res)

## (Intercept)          sw 
##   2.6390012   0.6904897

head(res$residuals)

##           1           2           3           4           5           6 
##  0.04428474  0.18952960 -0.14856834 -0.17951937 -0.12476423  0.06808885

head(resid(res))

##           1           2           3           4           5           6 
##  0.04428474  0.18952960 -0.14856834 -0.17951937 -0.12476423  0.06808885

head(res$fitted.values)

##        1        2        3        4        5        6 
## 5.055715 4.710470 4.848568 4.779519 5.124764 5.331911

head(predict(res))

##        1        2        3        4        5        6 
## 5.055715 4.710470 4.848568 4.779519 5.124764 5.331911

trees 자료에 대한 단순 선형회귀

trees 자료: 31 그루의 체리 나무에 대한 직경(Girth), 높이(Height), 부피(Volume) 측정 자료

다음은, trees 자료에서, Volume 을 종속변수로 하고, Girth 를 독립변수로 하는 단순 선형회귀 분석이다

dim(trees)

## [1] 31  3

head(trees)

##   Girth Height Volume
## 1   8.3     70   10.3
## 2   8.6     65   10.3
## 3   8.8     63   10.2
## 4  10.5     72   16.4
## 5  10.7     81   18.8
## 6  10.8     83   19.7

res <- lm(Volume~Girth,trees)
summary(res)

## 
## Call:
## lm(formula = Volume ~ Girth, data = trees)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.065 -3.107  0.152  3.495  9.587 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -36.9435     3.3651  -10.98 7.62e-12 ***
## Girth         5.0659     0.2474   20.48  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.252 on 29 degrees of freedom
## Multiple R-squared:  0.9353, Adjusted R-squared:  0.9331 
## F-statistic: 419.4 on 1 and 29 DF,  p-value: < 2.2e-16

anova(res)

## Analysis of Variance Table
## 
## Response: Volume
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Girth      1 7581.8  7581.8  419.36 < 2.2e-16 ***
## Residuals 29  524.3    18.1                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

H <- cut(trees$Height,3)
p <- ggplot(trees)+geom_point(aes(x=Girth, y=Volume, color= H))
p+geom_abline(intercept=coef(res)[1],slope=coef(res)[2],col='red')

원점을 지나는 회귀모형

원점을 지나는 회귀모형:
- y-절편에 해당하는 $\alpha$ 값이 없는 (0으로 설정된) 모형
- $y_j = \beta x_j +\epsilon_j,~~\epsilon_j ~\stackrel{iid}{\sim}~ N(0,\sigma^2 ), ~~j=1,2,\ldots,n$

원점을 지나는 회귀모형에서의 ANOVA Table:
- $y$-절편이 있는 모형에 대한 ANOVA Table 에서의 비교 모형:
  - Reduced Model (Model A): $~~y_j = \alpha + \epsilon_j,$
  - Full Model (Model AB): $~~y_j = \alpha + \beta x_j +\epsilon_j,$
- 원점을 지나는 모형에 대한 ANOVA Table 에서의 비교 모형:
  - Reduced Model (Model 0): $~~y_j = \epsilon_j,$
  - Full Model (Model B): $~~y_j = \beta x_j +\epsilon_j,$
- 원점을 지나는 모형에서의 제곱합과 결정계수:
  - SST = SSE of Reduced Model = $\sum_j y_j^2$
  - SSE = SSE of Full Model = $\sum_j (y_j- \hat{\beta} x_j )^2$
  - $R^2 = SSR/SST$, $~~SSR=SST-SSE$
  - $R^2$ : $y$-절편이 있는 경우와 비교하여 사용하지 말자

’원점을 지나는 회귀모형’에 대한 해석에서 주의할 점:
- $R^2$ 이 크다고 하여 좋은 모형이라고 말할 수 없다
- 원점이 ’지렛점’과 같은 역할을 하여, 왜곡된 추정이 될 수 있다
- ‘지렛점’ : 하나의 관측점이 회귀분석 결과에 큰 영향을 미치는 경우

’$y$-절편이 있는 모형’과 ’원점을 지나는 모형’에 대한 비교 방법?
- Model B 는, Model AB 에서 $\alpha=0$ 인 조건이 부과된 특수한 경우이므로
  
  Reduced Model 과 Full Model 관계가 성립한다
- 때문에 ANOVA 에 근거하여, 위 두 모형을 가설검정 방법으로 비교할 수 있다.
  
  즉, $H_0 : \alpha=0$, vs. $H_1 : \alpha \neq 0~$ 에 대한 가설검정이 가능하다

Model B 와 Model AB 를 비교하는 다른 방법?
- 가설검정 방법은, 단순한 모형을 기본모형($H_0$)으로 잡고 기본모형을 선호하여 결론을 내린다
- 가설검정 방법을 사용하면, 기본모형인 Reduced Model 에 편향되는 결론을 얻게 되고,
  
  복잡한 모형인 Full Model 의 장점이 명확한 경우만 Full Model 을 선택하게 된다
- Model B 와 Model AB 의 비교에서, $y$-절편이 없는 Model B 는 기본모형($H_0$)으로
  
  설정하기에는 보편성이 부족한 모형이라고 보인다.
- 때문에, Model B 와 Model AB 를 비교하기 위해서는 가설검정 이외의 방법도 필요하다
- 예를 들면, AIC 혹은 BIC 등등 일반적인 모형선택 방법이 사용될 수 있다

R 의 lm 함수에서 모형설정 방법, 사용례:
- lm(y~1) : $~y_j = \alpha +\epsilon_j,~$ 설명변수 없이, y-절편만 있는 모형
- lm(y~ x) : $~y_j = \alpha+\beta x_j +\epsilon_j,~$ y-절편 있는 모형
- lm(y~ 1+x), lm(y~ x+1) : 위의 y-절편 있는 모형과 동일
- lm(y~ 0+x), lm(y~ x+0) : $~y_j = \beta x_j +\epsilon_j,~$ 원점을 지나는 모형
- lm(y~ -1+x), lm(y~x-1) : 위의 원점을 지나는 모형과 동일

원점을 지나는 회귀모형을 적용하는 경우 문제가 될 만한 경우에 해당하는 자료를
```
시뮬레이션 방법으로 생성하여 이를 이용하여 비교하여 살펴보자 
```

library(ggplot2)
set.seed(777)
x<- sample(seq(2,3,l=10),20,repl=TRUE )
y<- 20*(sqrt(x)+0.1*rnorm(20))

b<- coef(lm(y~x))
s<- coef(lm(y~x-1))

u<-  seq(0,3.2,l=200)
w1<- sqrt(u)*coef(lm(y~I(sqrt(x))-1))
w2<- b[1] + b[2]*u
w3<- s*u

ggplot()+geom_point(aes(x=x,y=y)) +xlim(0,3.5) +ylim(0,42)+
         geom_line(aes(x=u,y=w1), col='orange') +
         geom_line(aes(x=u,y=w2), col='red') + 
         geom_line(aes(x=u,y=w3), col='blue')

summary( res0 <- lm(y~x-1) )

## 
## Call:
## lm(formula = y ~ x - 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0805 -1.9582 -0.3559  1.9082  7.1360 
## 
## Coefficients:
##   Estimate Std. Error t value Pr(>|t|)    
## x  11.9987     0.2559   46.89   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.014 on 19 degrees of freedom
## Multiple R-squared:  0.9914, Adjusted R-squared:  0.991 
## F-statistic:  2199 on 1 and 19 DF,  p-value: < 2.2e-16

summary( res1<- lm(y~x)    )

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5182 -1.2646 -0.4388  1.3524  5.2445 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   15.771      4.781   3.299  0.00399 **
## x              6.051      1.815   3.333  0.00370 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.445 on 18 degrees of freedom
## Multiple R-squared:  0.3817, Adjusted R-squared:  0.3473 
## F-statistic: 11.11 on 1 and 18 DF,  p-value: 0.003699

위 시뮬레이션 자료에 대하여 ANOVA 를 이용한 가설검정 방법을 적용해 보자
- anova(res0) : Model 0 과 Model B 의 비교 (결론: Model B 선호)
- anova(res1) : Model A 와 Model AB 의 비교 (결론: Model AB 선호)
- anova(res0,res1) : Model B 와 Model AB 의 비교 (결론: Model AB 선호)

AIC, BIC 에 의한 비교:
- Model B 와 Model AB 의 비교 (결론: Model AB 선호)

anova(res0)       # Comparison: Model 0, Model B

## Analysis of Variance Table
## 
## Response: y
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## x          1 19976.2 19976.2  2198.6 < 2.2e-16 ***
## Residuals 19   172.6     9.1                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(res1)       # Comparison: Model A, Model AB

## Analysis of Variance Table
## 
## Response: y
##           Df  Sum Sq Mean Sq F value   Pr(>F)   
## x          1  66.417  66.417  11.111 0.003699 **
## Residuals 18 107.595   5.977                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(res0,res1)  # Comparison: Model B, Model AB

## Analysis of Variance Table
## 
## Model 1: y ~ x - 1
## Model 2: y ~ x
##   Res.Df    RSS Df Sum of Sq     F   Pr(>F)   
## 1     19 172.63                               
## 2     18 107.59  1    65.038 10.88 0.003994 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

c(AIC(res0), BIC(res0))

## [1] 103.8662 105.8577

c(AIC(res1), BIC(res1))

## [1] 96.41029 99.39749

모형선택 기준 AIC 와 BIC 에 대하여는 나중에 학습하게 됨

AIC, BIC 는 낮은 값을 갖는 모형이 더 좋은 모형임

’y-절편 이 있는 모형’과 ’원점을 지나는 모형’의 비교 선택:
- $R^2$ 로 두 모형을 비교하면 안됨.
- ANOVA 에 의하여 두 모형을 직접 비교할 수 있음.
- 차후에 학습하게 될, AIC, BIC 등의 모형선택 기준을 적용할 수 있음
- ’원점을 지나는 모형’을 선택하려면,
  - 논리적으로 원점을 지나는 것이 타당해야 하고,
  - 관측점들이 원점에서 멀리 떨어져 있는 경우가 아니어야 하고,
  - SSE 는 항상 커지게 되지만, 그 변화량이 가능한 작아야 하고,
  - MSE 의 변화가 크지 않거나 혹은 작아지는 경우가 바람직

Lawn roller data

Lawn roller data
- DAAG 패키지에 포함된 자료
- 잔디밭 롤러의 무게와 파인 깊이에 관한 10회의 기록
- 변수:
  - weight: 사용한 롤러의 무게
  - depression: 잔디밭이 파인 깊이

비교 모형:
- rout1: y-절편이 있는 모형
- rout0: y-절편이 없는 모형(원점을 지나는 모형)

Lawn roller data 에서 원점을 지나는 모형:
- y-절편이 없는 모형이, y-절편이 있는 모형에 비하여
- SSE 는 더 커지지만, residual 자유도가 1 증가한 것에 비하여,
  
  SSE 의 증가가 적어서 MSE 는 줄어들고 있다
- 관측점들이 원점에서 약간 떨어져 있으나, 아주 멀리 떨어져 있는 것은 아니다
- ANOVA 에 의하여 y-절편이 없는 모형이, y-절편이 있는 모형을 비교하면
  
  y-절편이 없는 모형을 선호하는 결론에 이른다
- 두 모형에 대한 AIC, BIC 값을 살펴보면, y-절편이 없는 모형이 더 낮은 값이다
  
  이런 점들에 근거하여, y-절편이 없는 모형이 더 좋다고 볼 수 있다

# install.packages("DAAG")
library(ggplot2)
ggplot(DAAG::roller)+geom_point(aes(x=weight,y=depression))+xlim(0,13)+ylim(0,32)

rout0<- lm(depression ~ -1+weight,data=DAAG::roller)
rout1<- lm(depression ~ weight,   data=DAAG::roller)
anova(rout0)

## Analysis of Variance Table
## 
## Response: depression
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## weight     1 2637.32  2637.3  63.862 2.233e-05 ***
## Residuals  9  371.68    41.3                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(rout1)

## Analysis of Variance Table
## 
## Response: depression
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## weight     1 657.97  657.97  14.503 0.005175 **
## Residuals  8 362.93   45.37                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(rout0,rout1)

## Analysis of Variance Table
## 
## Model 1: depression ~ -1 + weight
## Model 2: depression ~ weight
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1      9 371.68                           
## 2      8 362.93  1    8.7433 0.1927 0.6723

c( AIC(rout0), BIC(rout0) )

## [1] 68.53317 69.13834

c( AIC(rout1), BIC(rout1) )

## [1] 70.29512 71.20288

ss0<-anova(rout0)$'Sum Sq'
ss1<-anova(rout1)$'Sum Sq'
ss0[1]/sum(ss0)   # R^2 no-intercept model

## [1] 0.8764782

ss1[1]/sum(ss1)   # R^2 with-intercept model

## [1] 0.6444963

실습 과제

trees 자료에서,
- Volume 을 종속변수로 하고, Height 를 독립변수로 하는 단순 선형회귀 분석 수행
- Summary Table 의 값을 보고 그 결과의 의미를 설명하시오.
- ANOVA Table 의 값을 보고 그 결과의 의미를 설명하시오.
- 다음과 같이, 회귀분석 결과를 res 에 저장하고, 변수, y, yhat, dfr 을 설정하고,
  - res <- lm(Volume~ Height, trees)
  - y <- trees$Volume
  - yhat <- predict(res)
  - dfr <- res$df.residual
- ANOVA Table 에 있는 모든 값들을 y, yhat, dfr 을 이용하여 구하시오.

trees 자료에서,
- ggplot2 패키지를 이용하여, (Height, Volume) 산점도를 그리고 회귀직선을 그리시오.
- ggplot2 패키지를 이용하여, (Girth, Volume) 산점도를 그리고 회귀직선을 그리시오.
- 다음 두 모형을 비교하여 평가하고,근거를 제시하시오.
  - lm(Volume~ Height, trees)
  - lm(Volume~ Height-1, trees)
- 다음 두 모형을 비교하여 평가하고,근거를 제시하시오.
  - lm(Volume~ Girth, trees)
  - lm(Volume~ Girth-1, trees)