회귀분석 Regression, Part III


Topics:


  • Part I:

    • iris 자료

    • 단순 선형회귀

  • Part II

    • 사영과 회귀분석

    • 다중 선형회귀

  • Part III

    • 가변수와 선형모형

    • 모형의 다양성

  • Part IV

    • 모형 선택

    • 회귀 진단



1. 가변수와 선형모형


  • 집단별 모형화의 필요성:

    • 아버지의 키와 자식의 키:

      • 딸인 경우: \(~Y_j = \alpha_0 + \beta \, x_j + \epsilon_j\),\(~~j \in\) 딸

      • 아들인 경우: \(~Y_j = \alpha_1 + \beta \, x_j + \epsilon_j\), \(~~j \in\) 아들

    • 붓꽃의 종류에 따른 sl 과 sw 의 관계:

      • setosa 의 경우 : \(~sl_j = \alpha_0 + \beta \, sw_j + \epsilon_j\), \(~~j = 1, \ldots,50\)

      • versicolor 의 경우 : \(~sl_j = \alpha_1 + \beta \, sw_j + \epsilon_j\), \(~~j = 1, \ldots,50\)

      • virginica 의 경우 : \(~sl_j = \alpha_2 + \beta \, sw_j + \epsilon_j\), \(~~j = 1, \ldots,50\)


  • 가변수를 이용한 모형의 통합

    • 가변수 dummy variable :

      • s 변수: 딸은 0, 아들은 1

      • (\(s_1\), \(s_2\)) 변수: setosa 는 (0,0), versicolor 는 (1,0), virginca 는 (0,1)

    • 통합 표현된 모형:

      • 아버지와 자식의 키: \(~Y_l = \gamma_0 + \gamma_1 \, s_l + \beta \, x_l + \epsilon_l\), \(~~l=1,\ldots, n\)

        • \(\alpha_0 = \gamma_0\)

        • \(\alpha_1 = \gamma_0 + \gamma_1\)

      • 붓꽃의 sl, sw 관계: \(~sl_l = \gamma_0 + \gamma_1 \, s_{1l} + \gamma_2 \, s_{2l} + \beta \, sw_l + \epsilon_l\), \(~~l=1,\ldots, 150\)

        • \(\alpha_0 = \gamma_0\)

        • \(\alpha_1 = \gamma_0 + \gamma_1\)

        • \(\alpha_2 = \gamma_0 + \gamma_2\)

    • 변수의 수식표현:

      • \(~s_l = I(l \in \mbox{아들})~\)

      • \(~s_{1l} = I(l \in \mbox{versicolor})~\), \(~s_{2l} = I(l \in \mbox{virginica})\)


  • 가변수의 활용:

    • 가변수의 개수: \(~g+1~\) 개의 집단이 있다면, \(~g~\) 개의 가변수 필요

    • 가변수를 도입하여, 다중 회귀모형과 ANOVA 모형을 통합할 수 있음

    • 선형모형: ANOVA 모형과 회귀모형을 통합한 모형



ANOVA


  • ANOVA 모형/분석:

    • \(g+1\) 개의 집단이 있고, 각 집단 \(~i=0,1,\ldots,g~\) 에 대하여,

      \(Y_{ij}= \alpha + \beta_{i}+ \epsilon_{ij}\), \(~~\epsilon_{ij}~\stackrel{iid}{\sim}~ N(0,\sigma^2)\), \(~~j=1,\ldots,k_i\)

    • \(g+1\) 개 집단 각각의 모평균에 차이가 있는 지를 알아보기 위한 통계적인 모형

    • 두 집단의 모평균에 차이가 있는 지를 알아보기 위한 t-검정 방법을 확장한 방법

    • 가변수로만 구성된 다중 선형회귀모형,

      • \(l=l(i,j)\), \(~~l=1,2,\ldots, n,~~\) 여기서 \(~n=\sum_{i=0}^g k_i\)

      • 모형: \(~~Y_l= \alpha + \beta_1 \, s_{1l} + \beta_2 \, s_{2l}+ \cdots+ \beta_g \, s_{gl}+ \epsilon_l\)

      • 가변수 \(s_{il}~\): \(~l\) 이 \(i\) 집단에 해당하면 1, 아니면 0, \(~i=1,\ldots,g\)

      • 각 집단의 모평균 \(~\mu_i~\),

        • \(\mu_0=\alpha\)

        • \(\mu_i = \alpha+\beta_i\), \(~~i=1,\ldots,g\)

    • Balanced ANOVA: 각 집단에서의 표본 크기 \(k_i\) 가 동일한 경우

    • One/two way ANOVA:

      • one-way: \(~Y_{ij}=\alpha +\beta_i + \epsilon_{ij}\)

      • two-way: \(~Y_{ijk}=\alpha +\beta_i + \gamma_j+ \epsilon_{ijk}\)

        • w/o replication: \(~k=1,~~~\) w replication: \(~k>1\)



sleep 데이터에 대한 ANOVA


  • sleep 데이터:

    • 10명의 환자에 대하여 수면 시간의 증가를 기록한 자료

      • extra: 수면시간 증가

      • group: 약 복용여부

      • ID: 환자 번호

    • William Gosset 은 1908년 Biometrika 에 Student 라는 가명으로 논문을 발표하여,

      t-분포의 개념과 그 사용을 제안함. [Student 논문], [참고 wiki]

    • sleep 데이터는 그 논문에 등장하는 예제 자료임. [help(sleep)]

    • sleep 데이터에 대한 해석:

      • 얼핏보면, 약 복용 여부에 따른 2 그룹에서의 수면시간 증가의 차이를

        살펴보기 위한 자료, 1-way ANOVA 를 하여야 하는 자료로 보임.

      • 잘 살펴보면, 환자를 의미하는 ID 변수가 있어서, 동일한 환자들에 대한

        random block 을 설정할 필요가 있음. 2-way ANOVA 를 하여야 하는 경우.

      • t-test 로 설명하면, 쌍체 (paired) 자료의 경우이므로, 동일한 환자에 대한

        그룹 1 과 그룹 2 사이의 차이를 구하여 처리해야할 자료.

        차이값 1 그룹에 대한 모평균이 0 인지를 검정하여야 하는 경우.


anova(lm(extra~ factor(group), sleep))
## Analysis of Variance Table
## 
## Response: extra
##               Df Sum Sq Mean Sq F value  Pr(>F)  
## factor(group)  1 12.482 12.4820  3.4626 0.07919 .
## Residuals     18 64.886  3.6048                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(extra~ factor(group)+factor(ID), sleep))
## Analysis of Variance Table
## 
## Response: extra
##               Df Sum Sq Mean Sq F value   Pr(>F)   
## factor(group)  1 12.482 12.4820 16.5009 0.002833 **
## factor(ID)     9 58.078  6.4531  8.5308 0.001901 **
## Residuals      9  6.808  0.7564                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
t.test(extra ~ group, paired=T, sleep)
## 
##  Paired t-test
## 
## data:  extra by group
## t = -4.0621, df = 9, p-value = 0.002833
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -2.4598858 -0.7001142
## sample estimates:
## mean difference 
##           -1.58
t.test(extra ~ group, sleep)
## 
##  Welch Two Sample t-test
## 
## data:  extra by group
## t = -1.8608, df = 17.776, p-value = 0.07939
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
##  -3.3654832  0.2054832
## sample estimates:
## mean in group 1 mean in group 2 
##            0.75            2.33
t.test(extra ~ group, var.equal=TRUE, sleep)
## 
##  Two Sample t-test
## 
## data:  extra by group
## t = -1.8608, df = 18, p-value = 0.07919
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
##  -3.363874  0.203874
## sample estimates:
## mean in group 1 mean in group 2 
##            0.75            2.33



가변수와 대비


  • 대비 contrast 의 필요성:

    • 예를 들어, 세 집단 \(~i=1,2,3~\) 인 경우를 살펴보자

    • 가변수 \(s_1\), \(s_2\), \(s_3\) : 각 집단에 속하는 경우에만 1, 다른 경우는 0

    • 평균 \(~\mu_i =\alpha + \beta_1 s_1 + \beta_2 s_2 +\beta_3 s_3\) 이면

      \(\mu_i = \alpha + \beta_i\), \(~i=1,2,3\)

    • 평균은 3개: \(~\mu_i\), \(~i=1,2,3\)

      추정해야 할 모수는 4 개: \(~\alpha\), \(\beta_i\), \(i=1,2,3\)

      4개의 모수를 동시에 추정할 수 없음.

      모수들 사이에 제약식 하나를 추가로 부가 해야 함.

      이 제약식이 contrast 임.

    • 집단의 개수 \(g+1\) 보다 가변수의 개수 \(g\) 가 하나 더 적다는 말은

      대비가 부가되어 있다는 말과 같음

    • 다음 4 개의 벡터 \(\mathbf{1}\), \(s_1\), \(s_2\), \(s_4\) 가 독립이 아니므로

      \(\alpha\), \(\beta_1\),\(\beta_2\), \(\beta_3\) 사이에 제약식을 주어, 하나를 제거하겠다는 의미임

    \(\qquad \mathbf{1} = \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix}, \quad s_1 = \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}, \quad s_2= \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix},\quad s_3 = \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix}\)


  • contrast 예 (세 집단의 경우):

    • \(y\)-절편 없는 모형:

      • \(\alpha=0\): 벡터 \(\mathbf{1}\) 이 빠지게 되므로, 가변수 (\(s_1\), \(s_2\),\(s_3\)) 모두 사용 \(\\[.1in]\)

    • \(y\)-절편 있는 모형:

      • treatment cotrast: \(\beta_1 =0,~\) 가변수 (\(s_2\),\(s_3\)) 사용

      • SAS contrast: \(\beta_3 =0,~\) 가변수 (\(s_1\),\(s_2\)) 사용

      • sum contrast: 가변수 (\(s_1 -s_3\), \(s_2 - s_3\)) 사용

      • helmert contrast: 가변수 (\(s_2 -s_1\), \(2 s_3 -s_1 - s_2\)) 사용

      • contrast matrix representation: 위 4 가지 순서대로 …


\(\qquad \qquad \qquad \begin{pmatrix} \beta_1 \\ \beta_2 \\ \beta_3 \end{pmatrix}= \begin{pmatrix} 0 & 0 \\ 1 & 0 \\ 0 & 1 \end{pmatrix}\begin{pmatrix} \delta_1 \\ \delta_2 \end{pmatrix},\qquad \quad \begin{pmatrix} \beta_1 \\ \beta_2 \\ \beta_3 \end{pmatrix}= \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ 0 & 0 \end{pmatrix}\begin{pmatrix} \delta_1 \\ \delta_2 \end{pmatrix}\)


\(\qquad \qquad \qquad \begin{pmatrix} \beta_1 \\ \beta_2 \\ \beta_3 \end{pmatrix}= \left( \begin{array}{rr} 1 & 0 \\ 0 & 1 \\ -1 & -1 \end{array} \right) \begin{pmatrix} \delta_1 \\ \delta_2 \end{pmatrix},\qquad \begin{pmatrix} \beta_1 \\ \beta_2 \\ \beta_3 \end{pmatrix}= \left(\begin{array}{rr} -1 & -1 \\ 1 & -1 \\ 0 & 2 \end{array}\right) \begin{pmatrix} \delta_1 \\ \delta_2 \end{pmatrix}\)



  • 4개의 집단에 대한 contrast:


contr.treatment(4)
##   2 3 4
## 1 0 0 0
## 2 1 0 0
## 3 0 1 0
## 4 0 0 1
contr.sum(4)
##   [,1] [,2] [,3]
## 1    1    0    0
## 2    0    1    0
## 3    0    0    1
## 4   -1   -1   -1
contr.helmert(4)
##   [,1] [,2] [,3]
## 1   -1   -1   -1
## 2    1   -1   -1
## 3    0    2   -1
## 4    0    0    3
contr.SAS(4)
##   1 2 3
## 1 1 0 0
## 2 0 1 0
## 3 0 0 1
## 4 0 0 0



  • R에서의 가변수 사용:

    • 명목변수: R에서는 factor 변수라고 함

    • 순서변수: R에서는 ordered 변수라고 함

    • 회귀분석에서, 설명변수로 사용되는 경우,

      • factor 변수: 자동으로 가변수를 대입

      • ordered 변수: polynomial contrast 로 변환된 변수 대입

    • R 함수:

      • options()$contrasts

      • contrast()

      • levels()

      • relevel()


  • factor 변수와 ordered 변수:

    • ordered 변수는 factor 변수이면서, factor levels 선후 관계가 있는 경우


mtcars$cyl
##  [1] 6 6 4 6 8 6 8 4 4 6 6 8 8 8 8 8 8 4 4 4 4 8 8 8 8 4 4 4 8 6 8 4
factor(mtcars$cyl)
##  [1] 6 6 4 6 8 6 8 4 4 6 6 8 8 8 8 8 8 4 4 4 4 8 8 8 8 4 4 4 8 6 8 4
## Levels: 4 6 8
ordered(mtcars$cyl)
##  [1] 6 6 4 6 8 6 8 4 4 6 6 8 8 8 8 8 8 4 4 4 4 8 8 8 8 4 4 4 8 6 8 4
## Levels: 4 < 6 < 8
is.factor(mtcars$cyl)
## [1] FALSE
is.ordered(mtcars$cyl)
## [1] FALSE
is.factor(factor(mtcars$cyl))
## [1] TRUE
is.ordered(ordered(mtcars$cyl))
## [1] TRUE
is.ordered(factor(mtcars$cyl))
## [1] FALSE
is.factor(ordered(mtcars$cyl))
## [1] TRUE



  • levels() 와 relevel()

    • levels() 를 이용하여, factor 변수 level 의 이름을 바꿀 수 있다

    • level 의 나열 순서는 디폴트는 알파벳 순서로 결정되고,

    • relevel() 을 이용하여 reference level 을 바꿀 수 있다

    • reference level 은 맨 앞에 오는 level 이다

    • 전체 level 순서를 바꿀 때는 factor(x,levels= …) 과 같은 방법 사용

    • 불필요한 level 값을 제거할 때는, x[,drop=T] 와 같은 방법 사용


iris<- iris; rm(iris)

names(iris)
## [1] "Sepal.Length" "Sepal.Width"  "Petal.Length" "Petal.Width"  "Species"
names(iris) <- c('sl','sw', 'pl','pw', 'sp')
names(iris)
## [1] "sl" "sw" "pl" "pw" "sp"
table(iris$sp)
## 
##     setosa versicolor  virginica 
##         50         50         50
(smp <-sample(iris$sp,30) )
##  [1] setosa     virginica  virginica  virginica  versicolor versicolor
##  [7] versicolor setosa     virginica  versicolor versicolor versicolor
## [13] setosa     setosa     versicolor versicolor setosa     setosa    
## [19] virginica  setosa     virginica  virginica  versicolor versicolor
## [25] setosa     setosa     setosa     versicolor virginica  setosa    
## Levels: setosa versicolor virginica
levels(smp)
## [1] "setosa"     "versicolor" "virginica"
levels(smp)<- c('st','vc','vg')

levels(smp)
## [1] "st" "vc" "vg"
smp
##  [1] st vg vg vg vc vc vc st vg vc vc vc st st vc vc st st vg st vg vg vc vc st
## [26] st st vc vg st
## Levels: st vc vg
sort(smp)
##  [1] st st st st st st st st st st st vc vc vc vc vc vc vc vc vc vc vc vg vg vg
## [26] vg vg vg vg vg
## Levels: st vc vg


smp <- relevel(smp, ref='vg')

levels(smp)
## [1] "vg" "st" "vc"
smp
##  [1] st vg vg vg vc vc vc st vg vc vc vc st st vc vc st st vg st vg vg vc vc st
## [26] st st vc vg st
## Levels: vg st vc
sort(smp)
##  [1] vg vg vg vg vg vg vg vg st st st st st st st st st st st vc vc vc vc vc vc
## [26] vc vc vc vc vc
## Levels: vg st vc
smp <- factor(smp, levels=levels(smp)[c(2,3,1)])

levels(smp)
## [1] "st" "vc" "vg"
sort(smp)
##  [1] st st st st st st st st st st st vc vc vc vc vc vc vc vc vc vc vc vg vg vg
## [26] vg vg vg vg vg
## Levels: st vc vg
tmp<- sort(smp)[1:10]
tmp
##  [1] st st st st st st st st st st
## Levels: st vc vg
tmp[,drop=T]
##  [1] st st st st st st st st st st
## Levels: st



  • 상수항(y-절편항) 있는 경우와 없는 경우 lm() 의 결과

    • 상수항이 있는 경우는, 각 level 집단의 평균은 상수항 만큼을 더한 값이다

    • 상수항이 있는 경우, refernce level 이 상수항으로 대체 된다


  • factor 변수가 있는 경우,

    • 각 집단의 효과(평균)를 알아볼 때는, 절편항이 없는 모형이 편리

    • ANOVA table 을 볼 때는, 절편항이 있는 모형이 편리


( res <- lm(sl~sp-1,iris) )
## 
## Call:
## lm(formula = sl ~ sp - 1, data = iris)
## 
## Coefficients:
##     spsetosa  spversicolor   spvirginica  
##        5.006         5.936         6.588
( rez <- lm(sl~sp,iris) )
## 
## Call:
## lm(formula = sl ~ sp, data = iris)
## 
## Coefficients:
##  (Intercept)  spversicolor   spvirginica  
##        5.006         0.930         1.582
coef(rez)[1]
## (Intercept) 
##       5.006
coef(rez)[1] + coef(rez)[2]
## (Intercept) 
##       5.936
coef(rez)[1] + coef(rez)[3]
## (Intercept) 
##       6.588
anova(res)
## Analysis of Variance Table
## 
## Response: sl
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## sp          3 5184.9 1728.30  6521.7 < 2.2e-16 ***
## Residuals 147   39.0    0.27                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(rez)
## Analysis of Variance Table
## 
## Response: sl
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## sp          2 63.212  31.606  119.26 < 2.2e-16 ***
## Residuals 147 38.956   0.265                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



  • options()$contrasts 를 통하여, contrast 설정 시스템 디폴트를 확인

    • 시스템 디폴트는

      • factor 변수는 contr.treatment 이고

      • ordered 변수는 contr.ploy 이다


  • lm() 의 결과에서 상수항(y-절편)으로 대체되는 level 은,

    • contr.treatment 를 사용하는 경우, reference level 이 되고

    • contr.SAS 를 사용하는 경우, 순서가 맨 뒤인 level 이 된다


  • y-절편과 Contrast 에 따른 회귀모형식의 변화:

    • w/o-intercept model:

    [old lecture note]

    • contr.treatment model:

    [old lecture note]

    • contr.sum model:

    [old lecture note]


options()$contrasts
##         unordered           ordered 
## "contr.treatment"      "contr.poly"
levels(iris$sp)
## [1] "setosa"     "versicolor" "virginica"
lm(sl~sp-1,iris)
## 
## Call:
## lm(formula = sl ~ sp - 1, data = iris)
## 
## Coefficients:
##     spsetosa  spversicolor   spvirginica  
##        5.006         5.936         6.588
contrasts(iris$sp) <- 'contr.SAS'
contrasts(iris$sp)
##            setosa versicolor
## setosa          1          0
## versicolor      0          1
## virginica       0          0
( res <- lm(sl~sp,iris)  )
## 
## Call:
## lm(formula = sl ~ sp, data = iris)
## 
## Coefficients:
##  (Intercept)      spsetosa  spversicolor  
##        6.588        -1.582        -0.652
coef(res)[1]+ coef(res)[2]              # setosa
## (Intercept) 
##       5.006
coef(res)[1]+ coef(res)[3]              # versicolor
## (Intercept) 
##       5.936
coef(res)[1]                            # virginica
## (Intercept) 
##       6.588
contrasts(iris$sp)<- 'contr.sum'
contrasts(iris$sp)
##            [,1] [,2]
## setosa        1    0
## versicolor    0    1
## virginica    -1   -1
(res <-lm(sl~sp,iris) )
## 
## Call:
## lm(formula = sl ~ sp, data = iris)
## 
## Coefficients:
## (Intercept)          sp1          sp2  
##     5.84333     -0.83733      0.09267
coef(res)[1]+ coef(res)[2]              # setosa
## (Intercept) 
##       5.006
coef(res)[1]+ coef(res)[3]              # versicolor
## (Intercept) 
##       5.936
coef(res)[1]-coef(res)[2]-coef(res)[3]  # virginica
## (Intercept) 
##       6.588


lm(sl~sp-1,iris)
## 
## Call:
## lm(formula = sl ~ sp - 1, data = iris)
## 
## Coefficients:
##     spsetosa  spversicolor   spvirginica  
##        5.006         5.936         6.588
contrasts(iris$sp)<- 'contr.helmert'
contrasts(iris$sp)
##            [,1] [,2]
## setosa       -1   -1
## versicolor    1   -1
## virginica     0    2
( res <- lm(sl~sp,iris) )
## 
## Call:
## lm(formula = sl ~ sp, data = iris)
## 
## Coefficients:
## (Intercept)          sp1          sp2  
##      5.8433       0.4650       0.3723
coef(res)[1] - coef(res)[2] - coef(res)[3]   # setosa
## (Intercept) 
##       5.006
coef(res)[1] + coef(res)[2] - coef(res)[3]   # versicolor
## (Intercept) 
##       5.936
coef(res)[1] + 2*coef(res)[3]                # virginica
## (Intercept) 
##       6.588
contrasts(iris$sp)<- 'contr.treatment'
contrasts(iris$sp)
##            versicolor virginica
## setosa              0         0
## versicolor          1         0
## virginica           0         1
( res <-lm(sl~sp,iris)  )
## 
## Call:
## lm(formula = sl ~ sp, data = iris)
## 
## Coefficients:
##  (Intercept)  spversicolor   spvirginica  
##        5.006         0.930         1.582
coef(res)[1]                   # setosa
## (Intercept) 
##       5.006
coef(res)[1] + coef(res)[2]    # versicolor
## (Intercept) 
##       5.936
coef(res)[1] + coef(res)[3]    # virginica
## (Intercept) 
##       6.588



  • factor 변수가 있는 경우의 ANOVA Table:

    • \(y\)-절편이 없는 모형:

      • factor 변수는 집단 수 \(g\) 만큼의 가변수로 바뀌어 대입되므로,

      • factor 변수의 자유도는 집단 수 \(g\) 와 동일하다

    • \(y\)-절편이 있는 모형:

      • factor 변수는 집단 수 \(g-1\) 만큼의 가변수로 바뀌어 대입되므로,

      • factor 변수의 자유도는 집단 수 \(g-1\) 과 동일하다


iris<-iris; rm(iris)
names(iris) <- c('sl','sw', 'pl','pw', 'sp')
levels(iris$sp)<- c('st','vc','vg')

s1 <- s2 <- s3 <-rep(0,150)
s1[iris$sp=='st'] <- 1      # 1 when setosa 
s2[iris$sp=='vc'] <- 1      # 1 when versicolor 
s3[iris$sp=='vg'] <- 1      # 1 when virginica

irix<- data.frame(iris,s1,s2,s3)

lm(sl~sw+sp-1,iris) 
## 
## Call:
## lm(formula = sl ~ sw + sp - 1, data = iris)
## 
## Coefficients:
##     sw    spst    spvc    spvg  
## 0.8036  2.2514  3.7101  4.1982
lm(sl~sw+s1+s2+s3-1,irix)
## 
## Call:
## lm(formula = sl ~ sw + s1 + s2 + s3 - 1, data = irix)
## 
## Coefficients:
##     sw      s1      s2      s3  
## 0.8036  2.2514  3.7101  4.1982
lm(sl~sw+sp,iris)
## 
## Call:
## lm(formula = sl ~ sw + sp, data = iris)
## 
## Coefficients:
## (Intercept)           sw         spvc         spvg  
##      2.2514       0.8036       1.4587       1.9468
lm(sl~sw+s2+s3,irix)
## 
## Call:
## lm(formula = sl ~ sw + s2 + s3, data = irix)
## 
## Coefficients:
## (Intercept)           sw           s2           s3  
##      2.2514       0.8036       1.4587       1.9468


anova(lm(sl~sw+sp-1,iris))  # w/o-intercept model
## Analysis of Variance Table
## 
## Response: sl
##            Df Sum Sq Mean Sq  F value    Pr(>F)    
## sw          1 4996.7  4996.7 26050.62 < 2.2e-16 ***
## sp          3  199.2    66.4   346.15 < 2.2e-16 ***
## Residuals 146   28.0     0.2                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sl~sw+sp,iris))    # w-intercept model
## Analysis of Variance Table
## 
## Response: sl
##            Df Sum Sq Mean Sq  F value  Pr(>F)    
## sw          1  1.412   1.412   7.3628 0.00746 ** 
## sp          2 72.752  36.376 189.6512 < 2e-16 ***
## Residuals 146 28.004   0.192                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



2. 모형의 다양성


  • 선형모형:

    • 모형이 추정할 (평균) 모수들에 대하여 선형함수인 경우

    • 선형함수: (절편 없는) 일차함수의 개념을 다차원 공간으로 확장한 것

    • 관측변수에 대한 선형함수를 말하는 것은 아니나, 간혹 혼용되는 경우가 있음

    • 예) 선형모형:

      • \(y_j = \beta_0 + \beta_1 x_j + \beta_2 x_j^2 + \beta_3 \, \log(x_j ) + \epsilon_j\)

      • \(\log(y_j) = \beta_0 + \beta_1 \sqrt{x_j } + \epsilon_j\)

      • 위 두 경우 모두, 독립변수 \(x\)에 대한 관점에서는 비선형임

    • 예) 비선형모형:

      • \(y_j = \beta_0 \cdot (1- \exp \{ | \beta_1 - \beta_2 x_j |^\gamma \} ) + \epsilon_j\)

      • \(y_j = \beta_1 x_j /(\beta_2 + x_j) + \epsilon_j\)


  • R에서 선형모형의 설정 방법:

    • 0, 1 : 상수항 (y-절편)의 없고 있음을 나타냄

    • 마침표(.) : 종속변수로 사용되지 않은, data.frame 내에 있는 모든 변수

    • +, - : 독립변수의 추가와 제거

    • I() : 함수적 변환을 통한 새로운 변수

    • 콜론(:) : 교호작용

    • *, /, ^ : 변수들 사이의 교호작용을 포함하는 복합변수 생성

    • offset() : 계수 추정을 하지 않고 추가하는 변수



lm 모형의 설정: trees 데이터의 예


names(trees) <- c('g', 'h', 'v')


lm(v~ ., trees)
## 
## Call:
## lm(formula = v ~ ., data = trees)
## 
## Coefficients:
## (Intercept)            g            h  
##    -57.9877       4.7082       0.3393
# by just seeing the coefficients, we can tell the model

coef(lm(v~ ., trees))
## (Intercept)           g           h 
## -57.9876589   4.7081605   0.3392512
coef(lm(v~ .-1, trees))
##          g          h 
##  5.0440083 -0.4773192
coef(lm(v~ .+0, trees))
##          g          h 
##  5.0440083 -0.4773192



  • 시험할 모형:

    • \(v_j = \beta_0 + \beta_1 \,g_j + \beta_2 \, h_j + \epsilon_j\)

    • \(v_j = \beta_0 + \beta_1 \,g_j + \beta_2 \, h_j + \beta_3 \, g_j \, h_j + \epsilon_j\)

    • \(v_j = \beta_0 + \beta_1 \,g_j + \beta_2 \, h_j + \beta_3 \, g_j \, h_j + \beta_4 \, g_j^2 \, h_j + \epsilon_j\)

    • \(v_j = \beta_0 + \beta_1 \, g_j^2 \, h_j + \epsilon_j\)

    • \(v_j = \beta_0 \, g_j^2 \, h_j + \epsilon_j\)


rez <- lm(v~ g+h , trees)
res <- lm(v~ g+h +I(g*h) , trees)
rew <- lm(v~ g+h + I(g*h) + I(g^2*h), trees)
rev <- lm(v~ I(g^2*h), trees) 
req <- lm(v~ -1+ I(g^2*h) , trees)  

summary(rez)
## 
## Call:
## lm(formula = v ~ g + h, data = trees)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.4065 -2.6493 -0.2876  2.2003  8.4847 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -57.9877     8.6382  -6.713 2.75e-07 ***
## g             4.7082     0.2643  17.816  < 2e-16 ***
## h             0.3393     0.1302   2.607   0.0145 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.882 on 28 degrees of freedom
## Multiple R-squared:  0.948,  Adjusted R-squared:  0.9442 
## F-statistic:   255 on 2 and 28 DF,  p-value: < 2.2e-16
summary(res)
## 
## Call:
## lm(formula = v ~ g + h + I(g * h), data = trees)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.5821 -1.0673  0.3026  1.5641  4.6649 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 69.39632   23.83575   2.911  0.00713 ** 
## g           -5.85585    1.92134  -3.048  0.00511 ** 
## h           -1.29708    0.30984  -4.186  0.00027 ***
## I(g * h)     0.13465    0.02438   5.524 7.48e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.709 on 27 degrees of freedom
## Multiple R-squared:  0.9756, Adjusted R-squared:  0.9728 
## F-statistic: 359.3 on 3 and 27 DF,  p-value: < 2.2e-16


summary(rew)
## 
## Call:
## lm(formula = v ~ g + h + I(g * h) + I(g^2 * h), data = trees)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.0401 -1.0784 -0.1449  1.7335  4.2313 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  3.626023  45.875470   0.079    0.938
## g           -0.523697   3.713835  -0.141    0.889
## h           -0.064359   0.801283  -0.080    0.937
## I(g * h)     0.009384   0.079104   0.119    0.906
## I(g^2 * h)   0.002010   0.001211   1.659    0.109
## 
## Residual standard error: 2.625 on 26 degrees of freedom
## Multiple R-squared:  0.9779, Adjusted R-squared:  0.9745 
## F-statistic: 287.7 on 4 and 26 DF,  p-value: < 2.2e-16
summary(rev)
## 
## Call:
## lm(formula = v ~ I(g^2 * h), data = trees)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.6195 -1.1002 -0.1656  1.7451  4.1976 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.977e-01  9.636e-01  -0.309     0.76    
## I(g^2 * h)   2.124e-03  5.949e-05  35.711   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.493 on 29 degrees of freedom
## Multiple R-squared:  0.9778, Adjusted R-squared:  0.977 
## F-statistic:  1275 on 1 and 29 DF,  p-value: < 2.2e-16
summary(req)
## 
## Call:
## lm(formula = v ~ -1 + I(g^2 * h), data = trees)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.6696 -1.0832 -0.3341  1.6045  4.2944 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## I(g^2 * h) 2.108e-03  2.722e-05   77.44   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.455 on 30 degrees of freedom
## Multiple R-squared:  0.995,  Adjusted R-squared:  0.9949 
## F-statistic:  5996 on 1 and 30 DF,  p-value: < 2.2e-16


anova(rew)
## Analysis of Variance Table
## 
## Response: v
##            Df Sum Sq Mean Sq   F value    Pr(>F)    
## g           1 7581.8  7581.8 1100.5828 < 2.2e-16 ***
## h           1  102.4   102.4   14.8618 0.0006814 ***
## I(g * h)    1  223.8   223.8   32.4934 5.352e-06 ***
## I(g^2 * h)  1   19.0    19.0    2.7534 0.1090651    
## Residuals  26  179.1     6.9                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(rev)
## Analysis of Variance Table
## 
## Response: v
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## I(g^2 * h)  1 7925.8  7925.8  1275.3 < 2.2e-16 ***
## Residuals  29  180.2     6.2                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(req)
## Analysis of Variance Table
## 
## Response: v
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## I(g^2 * h)  1  36144   36144  5996.4 < 2.2e-16 ***
## Residuals  30    181       6                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(req,rev)
## Analysis of Variance Table
## 
## Model 1: v ~ -1 + I(g^2 * h)
## Model 2: v ~ I(g^2 * h)
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     30 180.83                           
## 2     29 180.24  1   0.59318 0.0954 0.7596


  • req 와 rev 에 대한 anova 함수 적용 결과

    p-value 가 매우 커서, 단순한 모형인 y-절편이 없는 모형을 선호


c(AIC(rew), BIC(rew))
## [1] 154.3487 162.9527
c(AIC(rev), BIC(rev))
## [1] 148.5429 152.8448
c(AIC(req), BIC(req))
## [1] 146.6447 149.5127


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

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



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


bz<- coef(rez)
bs<- coef(res)
bq<- coef(req)

vz <- function(w){ bz[1]+bz[2]*w[1]+bz[3]*w[2] }
vs <- function(w){ bs[1]+bs[2]*w[1]+bs[3]*w[2] +bs[4]*w[1]*w[2] }
vq <- function(w){ bq[1]*w[1]*w[1]*w[2] }

gsn <- seq( min(trees$g),max(trees$g), l=31)
hsn <- seq( min(trees$h),max(trees$h), l=51) 

ww <- expand.grid(g=gsn, h=hsn)

zmatz <- t(matrix(apply( ww, 1, vz), ncol=51 ))
zmats <- t(matrix(apply( ww, 1, vs), ncol=51 ))
zmatq <- t(matrix(apply( ww, 1, vq ), ncol=51 ))

pp <- plot_ly(z=~zmats,x=gsn, y=hsn) %>% add_surface(showscale=FALSE, opacity=0.8) %>%
           
         add_surface(showscale=FALSE, z=~zmatz,x=gsn, y=hsn,color=I('red')) %>%
  
         add_markers(x=trees$g, y=trees$h, z=trees$v, color=I('orange'))  

pp %>% layout(title="trees data", 
               scene=list(xaxis=list(title="g"), yaxis=list(title="h"), 
                          zaxis=list(title="v")
              ))


pp <- plot_ly(z=~zmatq,x=gsn, y=hsn) %>% add_surface(showscale=FALSE, opacity=0.8) %>%
           
         add_surface(showscale=FALSE, z=~zmatz,x=gsn, y=hsn, color=I('red')) %>%
  
         add_markers(x=trees$g, y=trees$h, z=trees$v, color=I('orange'))  

pp %>% layout(title="trees data", 
               scene=list(xaxis=list(title="g"), yaxis=list(title="h"), 
                          zaxis=list(title="v")
              ))



mtcars 데이터에서의 교호작용


  • 교호작용interaction 이 없는 모형과 교호작용이 있는 모형의 예:

    • 모형:

      • 교호작용이 없는 모형: 추정된 두 회귀직선의 기울기가 같음

      • 교호작용이 있는 모형: 추정된 두 회귀직선의 기울기가 다름

    • 효과:

      • 회귀직선의 기울기: disp 가 mpg 에 미치는 효과를 의미함,

      • 교호작용: 회귀직선의 기울기 (효과) 가 am 의 값에 따라 달라짐


head(mtcars)
##                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1


is.factor(mtcars$am)
## [1] FALSE
mtcars$amx <- factor(mtcars$am)
is.factor(mtcars$amx)
## [1] TRUE


library(ggplot2)

pp <- ggplot(mtcars)+geom_point(aes(x=disp,y=mpg,col=amx))

pp+ stat_smooth(aes(x=disp,y=mpg,col=amx), method='lm', se=FALSE) # +facet_wrap(~amx)
## `geom_smooth()` using formula 'y ~ x'


rez <- lm(mpg~ amx+disp, mtcars)           # w/o-interaction model
res <- lm(mpg~ amx+disp +am:disp, mtcars)  # w-interaction model

summary(rez)
## 
## Call:
## lm(formula = mpg ~ amx + disp, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.6382 -2.4751 -0.5631  2.2333  6.8386 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 27.848081   1.834071  15.184 2.45e-15 ***
## amx1         1.833458   1.436100   1.277    0.212    
## disp        -0.036851   0.005782  -6.373 5.75e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.218 on 29 degrees of freedom
## Multiple R-squared:  0.7333, Adjusted R-squared:  0.7149 
## F-statistic: 39.87 on 2 and 29 DF,  p-value: 4.749e-09
summary(res)
## 
## Call:
## lm(formula = mpg ~ amx + disp + am:disp, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.6056 -2.1022 -0.8681  2.2894  5.2315 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 25.157064   1.925053  13.068 1.94e-13 ***
## amx1         7.709073   2.502677   3.080  0.00460 ** 
## disp        -0.027584   0.006219  -4.435  0.00013 ***
## disp:am     -0.031455   0.011457  -2.745  0.01044 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.907 on 28 degrees of freedom
## Multiple R-squared:  0.7899, Adjusted R-squared:  0.7674 
## F-statistic: 35.09 on 3 and 28 DF,  p-value: 1.27e-09


a0 <- coef(rez)[1]
a1 <- coef(rez)[1]+coef(rez)[2]
b <-  coef(rez)[3]
  
pp + geom_abline(intercept=a0, slope=b, color='darkturquoise', size=1) +
      
    geom_abline(intercept=a1, slope=b, color='tomato', size=1)



  • 교호작용 콜론(:)의 계산:

    • 1:x 는 1 과 같은 의미,

    • x:x, x:x:x 는 x 와 같은 의미

    • x:z:x 는 x:z 와 같은 의미

    • x:z 는 z:x 와 같고, I(x*z) 를 대입한 것과 같은 결과



  • 양적변수 사이의 교호작용의 의미


coef( lm(v~ g+h +I(g*h) , trees) )
## (Intercept)           g           h    I(g * h) 
##  69.3963156  -5.8558479  -1.2970834   0.1346544
coef( lm(v~ g+h + g:h , trees) )
## (Intercept)           g           h         g:h 
##  69.3963156  -5.8558479  -1.2970834   0.1346544



  • 교호작용 콜론(:)의 계산 예1:

    • mtcars 데이터 에서 2 범주변수의 예

      • am : 0,1 값을 갖는 양적변수 (가변수 역할)

      • amx : 0, 1 값을 갖는 명목변수


coef( lm(mpg~ am +disp, mtcars) )
## (Intercept)          am        disp 
## 27.84808111  1.83345825 -0.03685086
coef( lm(mpg~ amx +disp, mtcars) )
## (Intercept)        amx1        disp 
## 27.84808111  1.83345825 -0.03685086
coef( lm(mpg~ am +disp + I(am*disp), mtcars) )
##  (Intercept)           am         disp I(am * disp) 
##  25.15706407   7.70907298  -0.02758360  -0.03145482
coef( lm(mpg~ am +disp + am:disp, mtcars) )
## (Intercept)          am        disp     am:disp 
## 25.15706407  7.70907298 -0.02758360 -0.03145482
coef( lm(mpg~ amx +disp + amx:disp, mtcars) )
## (Intercept)        amx1        disp   amx1:disp 
## 25.15706407  7.70907298 -0.02758360 -0.03145482



lm 모형의 다양성: iris 데이터의 예


  • 교호작용 콜론(:)의 계산 예2:

    • iris 데이터, 다범주 명목변수의 예

      • sp : ‘st’,‘v’c’,‘vg’ 값을 갖는 명목변수

      • s1, s2, s3 : 0, 1 값을 갖는 양적변수 형태의 가변수


names(iris) <- c('sl','sw', 'pl','pw', 'sp')
levels(iris$sp)<- c('st','vc','vg')

s1 <- s2 <- s3 <-rep(0,150)
s1[iris$sp=='st'] <- 1      # 1 when setosa 
s2[iris$sp=='vc'] <- 1      # 1 when versicolor 
s3[iris$sp=='vg'] <- 1      # 1 when virginica

irix<- data.frame(iris,s1,s2,s3)

contrasts(iris$sp)
##    vc vg
## st  0  0
## vc  1  0
## vg  0  1
coef(  lm(sl~  1 + sp + sw + sp:sw , iris)   )
## (Intercept)        spvc        spvg          sw     spvc:sw     spvg:sw 
##   2.6390012   0.9007335   1.2678352   0.6904897   0.1745880   0.2110448
coef(  lm(sl~  1 + s2 + s3 + sw + I(sw*s2) + I(sw*s3), irix)  )
## (Intercept)          s2          s3          sw  I(sw * s2)  I(sw * s3) 
##   2.6390012   0.9007335   1.2678352   0.6904897   0.1745880   0.2110448
coef(  lm(sl~  0 + sp + sp:sw , iris)   )
##      spst      spvc      spvg   spst:sw   spvc:sw   spvg:sw 
## 2.6390012 3.5397347 3.9068365 0.6904897 0.8650777 0.9015345
coef(  lm(sl~  0 + s1 + s2 + s3 + I(sw*s1) + I(sw*s2) + I(sw*s3), irix)  )
##         s1         s2         s3 I(sw * s1) I(sw * s2) I(sw * s3) 
##  2.6390012  3.5397347  3.9068365  0.6904897  0.8650777  0.9015345



  • 복합항의 대입 ( * 의 경우) :

    • sw*pl 은 sw+pl+sw:pl 과 같은 의미

    • pl*sw 는 pl+sw+pl:sw 와 같은 의미

    • 위 두 경우는 결과에 나타나는 변수 순서만 다를 뿐 동일한 모형

    • sp*sw 는 sp+sw+sp:sw 와 같고, 1, spvc, spvg, sw, spvc:sw, spvg:sw 가 포함됨

    • sp*sw-1 는 -1+sp+sw+sp:sw 와 같고, spst, spvc, spvg, sw, spvc:sw, spvg:sw 가 포함됨

    • sp*sw-1-sw 는 -1+sp+sp:sw 와 같고, spst, spvc, spvg, spst:sw, spvc:sw, spvg:sw 가 포함됨


coef(lm(sl~ sw*pl,iris))
## (Intercept)          sw          pl       sw:pl 
##  1.40438275  0.84995691  0.71845958 -0.07701327
coef(lm(sl~ pl*sw,iris))
## (Intercept)          pl          sw       pl:sw 
##  1.40438275  0.71845958  0.84995691 -0.07701327
coef(lm(sl~ sw*sp,iris))
## (Intercept)          sw        spvc        spvg     sw:spvc     sw:spvg 
##   2.6390012   0.6904897   0.9007335   1.2678352   0.1745880   0.2110448
coef(lm(sl~ sp*sw,iris))
## (Intercept)        spvc        spvg          sw     spvc:sw     spvg:sw 
##   2.6390012   0.9007335   1.2678352   0.6904897   0.1745880   0.2110448
coef(lm(sl~ sp*sw-1,iris))
##      spst      spvc      spvg        sw   spvc:sw   spvg:sw 
## 2.6390012 3.5397347 3.9068365 0.6904897 0.1745880 0.2110448
coef(lm(sl~ sp*sw-1-sw,iris))
##      spst      spvc      spvg   spst:sw   spvc:sw   spvg:sw 
## 2.6390012 3.5397347 3.9068365 0.6904897 0.8650777 0.9015345



  • 복합항의 대입 ( / 의 경우) :

    • sw/pl 은 sw+sw:pl 과 같은 의미

    • pl/sw 는 pl+pl:sw 와 같은 의미

    • 1/sw 은 1+1:sw 와 같은 의미이므로, 결국 1 이 됨

    • sw/1 는 sw+sw:1 과 같은 의미이므로, 1, sw 가 됨

    • sw/sp 는 sw + sw:sp 와 같은 의미로, 1, sw, sw:spvc, sw:spvg 항이 포함됨

    • sp/sw 는 sp + sp:sw 와 같은 의미로, 1, spvc, spvg, spst:sw, spvc:sw, spvg:sw 항이 포함됨

    • sp/sw-1 는 -1+sp + sp:sw 와 같은 의미로, spst, spvc, spvg, spst:sw, spvc:sw, spvg:sw 항이 포함됨


coef( lm(sl~ sw/pl,iris) )
## (Intercept)          sw       sw:pl 
##   3.9267848   0.1009329   0.1440569
coef( lm(sl~ pl/sw,iris) )
## (Intercept)          pl       pl:sw 
##  4.24809988  0.03403512  0.13145758
coef( lm(sl~ 1/sw,iris) )
## (Intercept) 
##    5.843333
coef( lm(sl~ sw/1,iris) )
## (Intercept)          sw 
##   6.5262226  -0.2233611
coef( lm(sl~ sw/sp,iris) )
## (Intercept)          sw     sw:spvc     sw:spvg 
##   3.3578882   0.4832626   0.4466483   0.6007515
coef( lm(sl~ sp/sw,iris) )
## (Intercept)        spvc        spvg     spst:sw     spvc:sw     spvg:sw 
##   2.6390012   0.9007335   1.2678352   0.6904897   0.8650777   0.9015345
coef( lm(sl~ sp/sw-1,iris) )
##      spst      spvc      spvg   spst:sw   spvc:sw   spvg:sw 
## 2.6390012 3.5397347 3.9068365 0.6904897 0.8650777 0.9015345



  • 복합항의 대입 ( ^ 의 경우) :

    • (sw+pl)^3 은 (sw+pl):(sw+pl):(sw+pl) 과 같은 의미로,

      전개하면 sw:sw:sw + sw:sw:pl + sw:pl:pl + pl:pl:pl 과 같다.

      따라서, 1, sw, pl, sw:pl 항이 포함된다

    • (sw+sp)^3 은 (sw+sp):(sw+sp):(sw+sp) 와 같은 의미로,

      1, sw, spvc, spvg, sw:spvc, sw:spvg 항이 포함된다

    • .^2 에서, 점(.) 은 종속변수로 사용된 변수를 제외한 모든 변수를 의미하므로,

      (sw+pl+pw+sp)^2 과 같다. 이를 정리하면,

      sw+pl+pw+sp+sw:pl+sw:pw+sw:sp+pl:pw+pl:sp+pw:sp 와 같고,

      1, sw, pl, pw, spvc, spvg, sw:pl, sw:pw, sw:spvc, sw:spvg,

      pl:pw, pl:spvc, pl:spvg, pw:spvc, pw:spvg 항들이 포함된다.


coef( lm(sl~ (sw+pl)^3,iris) )
## (Intercept)          sw          pl       sw:pl 
##  1.40438275  0.84995691  0.71845958 -0.07701327
coef(  lm(sl~ (sw+pl):(sw+pl):(sw+pl),iris) )
## (Intercept)          sw          pl       sw:pl 
##  1.40438275  0.84995691  0.71845958 -0.07701327
coef( lm(sl~ (sw+sp)^3,iris) )
## (Intercept)          sw        spvc        spvg     sw:spvc     sw:spvg 
##   2.6390012   0.6904897   0.9007335   1.2678352   0.1745880   0.2110448
coef( lm(sl~ .^2,iris) )
## (Intercept)          sw          pl          pw        spvc        spvg 
##   1.6997802   0.8300898   0.3178192   2.5827178  -2.7193358  -6.1704174 
##       sw:pl       sw:pw     sw:spvc     sw:spvg       pl:pw     pl:spvc 
##  -0.0148670  -0.6111773   0.4300016   0.8287899  -0.1194774   0.7397631 
##     pl:spvg     pw:spvc     pw:spvg 
##   0.9033624  -1.0069681  -0.2863915



coef( lm(sl~ (sw+pl+pw+sp)^3, iris) )
## (Intercept)          sw          pl          pw        spvc        spvg 
##  -0.5564029   1.2610318   2.0532798   2.2297899  11.3195632  10.3873300 
##       sw:pl       sw:pw     sw:spvc     sw:spvg       pl:pw     pl:spvc 
##  -0.3710377   0.4043050  -4.2044566  -5.0382366  -0.3865094  -4.3819250 
##     pl:spvg     pw:spvc     pw:spvg    sw:pl:pw  sw:pl:spvc  sw:pl:spvg 
##  -3.4780311  -5.1536246  -7.3548682  -0.5121047   1.5629123   1.3591013 
##  sw:pw:spvc  sw:pw:spvg  pl:pw:spvc  pl:pw:spvg 
##   0.4907695   1.7890774   1.7661169   1.6249096
coef( lm(sl~ (sw+pl+pw+sp+log(sw)+log(pl)+log(pw))^2, iris) )
##     (Intercept)              sw              pl              pw            spvc 
##     -97.5971570     123.1649962      -6.6737522       1.0705062     -17.9819837 
##            spvg         log(sw)         log(pl)         log(pw)           sw:pl 
##     -39.6294327    -129.5973203      34.8626954      -5.8154524       4.3218692 
##           sw:pw         sw:spvc         sw:spvg      sw:log(sw)      sw:log(pl) 
##     -12.4926314     -21.7963924     -29.6480766     -37.5424713      10.2568382 
##      sw:log(pw)           pl:pw         pl:spvc         pl:spvg      pl:log(sw) 
##       4.8585458      -2.1902224      -3.8122168      -9.1347913      -8.8123216 
##      pl:log(pl)      pl:log(pw)         pw:spvc         pw:spvg      pw:log(sw) 
##       3.6338651       6.8773112       2.9340706      11.8840411      36.8639325 
##      pw:log(pl)      pw:log(pw)    spvc:log(sw)    spvg:log(sw)    spvc:log(pl) 
##      -4.4048917       0.4765823      81.5878652     104.2385012       0.9411779 
##    spvg:log(pl)    spvc:log(pw)    spvg:log(pw) log(sw):log(pl) log(sw):log(pw) 
##      24.1251859      -3.1181759     -13.9260811     -53.7006555     -14.7325107 
## log(pl):log(pw) 
##      -6.4626084



  • 다양한 선형모형의 구성:

    • iris 데이터에 포함된 설명변수는 sw, pl, pw, sp 네 개다

    • 선형모형에는 이 설명변수들로부터 얻어지는 다양한 함수형태가 포함될 수 있으므로

      매우 많은 수의 설명변수를 모형에 포함시킬 수 있고, 다양한 종류의 선형모형을 구성할 수 있다

      선형모형에서 설명변수를 추가하면서 매우 다양한 모형을 구성할 수 있지만,

      선형모형 잔차의 자유도는 \(n-r\) 이므로 ( 여기서 \(~r=rank(X)~\)),

      추정가능한 모수의 개수가 관측 자료의 개수 \(n\) 이상이 될 수 없다는 절대적 한계가 있다

    • 절대적 한계와 별도로, 가능하다면 작은 개수의 설명변수를 사용하는 것이 바람직 하다

      • Ocam’s razor: 현상을 설명하는 모형 중 가장 단순한 모형이 가장 설득력이 높다

      • AIC : 사용한 독립변수 개수 이상으로 log-likelihood 값을 증가시켜야 한다

    • 바람직한 모형선택 방법에 대하여는 이후 별도의 주제로 학습하게 된다



offset 함수를 이용한 가설검정


  • 독립변수 항에서, offset() 함수가 적용된 변수는 그 회귀 계수를 추정하지 않고 모형에 대입한다

    • 예:

      • R model: ’y ~ x1 + offset(1.5*x2)’

      • 의미: \(~Y_j = \beta_0 + \beta_1 \, x_{1j} + 1.5 \, x_{2j} + \epsilon_j\)


  • 선형모형의 다양한 표현법을 이용한 일반적인 선형가설의 검정:

    • 가설과 F-검정:

      • Reduced Model (Model R): \(~~H_0 : C \boldsymbol{\beta}= \mathbf{d}\)

      • Full Model (Model F): \(~~H_1 : C \boldsymbol{\beta} \neq \mathbf{d}\)

      • \(Q(\mathbf{d}) = SSE(Model \,\,R) - SSE(Model \,\, F)~\) 인 관계가 성립



  • 가설검정 예1:

    • 모형: \(~y_j = \beta_0 + \beta_1 x_{1j}+\beta_2 x_{2j}+\beta_3 x_{3j}+\epsilon_j\)

    • 가설의 표현 (귀무가설):

      • \(H_0 : ~\beta_1 + \beta_3 = 0.1\), 이고 \(~\beta_2 = 0.75\)

      • \(H_0 : C \boldsymbol{\beta}= \mathbf{d}~~\), 여기서 \(~ C=\begin{pmatrix} 0 & 1 & 0 & 1 \\ 0 & 0 & 1 & 0 \end{pmatrix}, ~~~ \mathbf{d}=\begin{pmatrix} 0.1 \\ 0.75 \end{pmatrix}\)

    • 모형에의 대입:

      • \(~\beta_1 + \beta_3 = 0.1~\) 은, \(~\beta_3 = 0.1-\beta_1~\) 이므로, I(sw-pw)+offset(0.1*pw)

      • \(~\beta_2 = 0.75\) 은, offset(0.75*pl)


coef( res1<-lm(sl~ sw+pl+pw, iris) )
## (Intercept)          sw          pl          pw 
##   1.8559975   0.6508372   0.7091320  -0.5564827
coef( res0<- lm(sl~ I(sw-pw)+offset(0.1*pw)+offset(0.75*pl), iris) )
## (Intercept)  I(sw - pw) 
##   1.5667025   0.7202354
anova(res0,res1)
## Analysis of Variance Table
## 
## Model 1: sl ~ I(sw - pw) + offset(0.1 * pw) + offset(0.75 * pl)
## Model 2: sl ~ sw + pl + pw
##   Res.Df    RSS Df Sum of Sq     F Pr(>F)
## 1    148 14.585                          
## 2    146 14.445  2   0.13952 0.705 0.4958


  • 위 에서 두 모형의 적합 결과 res0 와 res1 을 비교하는 anova table 은

    res1 에서의 잔차제곱합은 0.13952 만큼 줄어들고, p-value 가 0.4958 임을 보여주고 있다.

    이 결과에 의하면, 가설 \(~H_0~\) 가 지지된다.



  • 가설검정 예2:

    • 모형: \(~y_j = \beta_0 + \beta_1 x_{1j}+\beta_2 x_{2j}+\epsilon_j\)

    • 가설의 표현 (귀무가설):

      • \(H_0 : ~\beta_1 - \beta_2 =0\), 이고 \(~\beta_0 = 2.3\)

      • \(H_0 : C \boldsymbol{\beta}= \mathbf{d}~~\), 여기서 \(~ C=\begin{pmatrix} 0 & 1 & -1 \\ 1 & 0 & 0 \end{pmatrix}, ~~~ \mathbf{d}=\begin{pmatrix}0 \\ 2.3 \end{pmatrix}\)

    • 모형에의 대입:

      \(y\)-절편 \(\beta_0\) 를 상수 2.3 으로 설정하기 위한 두 가지 방법:

      • 종속변수 항 \(y\) 에서 2.3 을 빼고, 원점을 지나는 회귀모형을 적용한다

      • 1 벡터에 해당하는 변수를 만들어, 예를 들어 one 이라는 이름을 붙이고,

        offset 함수를 이용하여 이 변수에 2.3 을 곱하여 모형에 추가한다


coef( res1<-lm(sl~ sw+pl, iris) )
## (Intercept)          sw          pl 
##   2.2491402   0.5955247   0.4719200
coef( lm(sl - 2.3 ~ 0+I(sw+pl), iris) )
## I(sw + pl) 
##  0.5169859
one<- rep(1,150)

coef( res0<- lm(sl~ 0+I(sw+pl)+offset(2.3*one), iris) )
## I(sw + pl) 
##  0.5169859
anova(res0,res1)
## Analysis of Variance Table
## 
## Model 1: sl ~ 0 + I(sw + pl) + offset(2.3 * one)
## Model 2: sl ~ sw + pl
##   Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
## 1    149 17.854                                
## 2    147 16.329  2    1.5247 6.8632 0.001414 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


  • 위 에서 두 모형의 적합 결과 res0 와 res1 을 비교하는 anova table 은

    res1 에서의 잔차제곱합이 1.5247 만큼 줄어들고 있음을 보여준다.

    p-value 가 0.0014 로 매우 작은 값이므로, 가설 \(~H_0~\) 를 기각하는 것을 지지한다



비선형 모형과 자료변환 방법


  • 비선형 회귀모형 및 적합방법:

    • 비선형 회귀모형:

      • 꼭 필요한 경우가 아니라면 사용하지 않는 것이 바람직

      • 비선형 회귀모형으로 해결하지 말고,

        • 가능하면 선형모형으로 해결하고,

        • 선형모형을 사용할 수 없는 경우라면,

          • 일반화선형모형을 적용하여 해결 (이후에 학습)

          • 자료변환 방법을 통하여 해결

    • 비선형 회귀모형 적합 방법:


  • 자료변환 방법:

    • 관측된 자료를 변환 하는 방법

      • \(x\)-변수 변환 방법: 다양한 함수적 변환 적용 가능

      • \(y\)-변수 변환 방법능: log 변환, 제곱근 변환, Box-Cox 변환

      • 양측 변환 방법: 독립변수와 종속변수 값의 물리적 동일성 유지가 필요한 경우

    • Box-Cox 변환법: [참고 wiki]

\(\qquad \qquad \qquad \quad T(y) = \left\{ \begin{array}{ll} (y^\lambda -1)/\lambda, & \lambda \neq 0 , \\ \log(y), & \lambda=0 \end{array} \right.\)


  • Cobb Douglas 모형:

    • \(Y= A L^{\beta_1} K^{\beta_2} \cdot U\)

    • \(\log(Y_j )= \alpha+ {\beta_1} \, \log(L_j) + \beta_2 \, \log(K_j)+\log(U_j), ~\) (여기서 \(~\alpha=\log(A)~\))

    • \(\log(U_j ) ~\stackrel{iid}{\sim}~ Normal(0 , \sigma^2 ),~\) (\(~U_j~\): 로그노말 분포를 따름 [참고 wiki] )



Forbes 데이터: 자료변환 방법의 예



MASS::forbes
##       bp  pres
## 1  194.5 20.79
## 2  194.3 20.79
## 3  197.9 22.40
## 4  198.4 22.67
## 5  199.4 23.15
## 6  199.9 23.35
## 7  200.9 23.89
## 8  201.1 23.99
## 9  201.4 24.02
## 10 201.3 24.01
## 11 203.6 25.14
## 12 204.6 26.57
## 13 209.5 28.49
## 14 208.6 27.76
## 15 210.7 29.04
## 16 211.9 29.88
## 17 212.2 30.06


( res <- lm(log(pres)~ bp, MASS::forbes))
## 
## Call:
## lm(formula = log(pres) ~ bp, data = MASS::forbes)
## 
## Coefficients:
## (Intercept)           bp  
##    -0.97087      0.02062


library(ggplot2)

ggplot(MASS::forbes)+geom_point(aes(x=bp,y=pres))+stat_smooth(method='lm', aes(x=bp,y=pres),col='blue')
## `geom_smooth()` using formula 'y ~ x'



trees 데이터: 자료변환 방법의 예


  • trees 데이터의 예: 다음 모형들을 비교하여 살펴보자

    • res0: \(~v ~\sim~ -1 +I(g^2*h)\)

    • res1: \(~log(v) ~\sim~ log(g) +log(h)\)

    • res2: \(~log(v) ~\sim~ I(2*log(g) +log(h))\)


res0 <- lm(v~ -1 + I(g^2 *h), trees )
res1 <- lm(log(v)~ 1 + log(g) + log(h), trees )
res2 <- lm(log(v)~ 1 + I(2*log(g) + log(h)), trees )

summary(res0)
## 
## Call:
## lm(formula = v ~ -1 + I(g^2 * h), data = trees)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.6696 -1.0832 -0.3341  1.6045  4.2944 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## I(g^2 * h) 2.108e-03  2.722e-05   77.44   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.455 on 30 degrees of freedom
## Multiple R-squared:  0.995,  Adjusted R-squared:  0.9949 
## F-statistic:  5996 on 1 and 30 DF,  p-value: < 2.2e-16
summary(res1)
## 
## Call:
## lm(formula = log(v) ~ 1 + log(g) + log(h), data = trees)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.168561 -0.048488  0.002431  0.063637  0.129223 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6.63162    0.79979  -8.292 5.06e-09 ***
## log(g)       1.98265    0.07501  26.432  < 2e-16 ***
## log(h)       1.11712    0.20444   5.464 7.81e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08139 on 28 degrees of freedom
## Multiple R-squared:  0.9777, Adjusted R-squared:  0.9761 
## F-statistic: 613.2 on 2 and 28 DF,  p-value: < 2.2e-16
summary(res2)
## 
## Call:
## lm(formula = log(v) ~ 1 + I(2 * log(g) + log(h)), data = trees)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.167714 -0.048284 -0.007407  0.064576  0.137522 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            -6.21387    0.26807  -23.18   <2e-16 ***
## I(2 * log(g) + log(h))  1.00473    0.02835   35.44   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08041 on 29 degrees of freedom
## Multiple R-squared:  0.9774, Adjusted R-squared:  0.9767 
## F-statistic:  1256 on 1 and 29 DF,  p-value: < 2.2e-16
c(AIC(res0),BIC(res0))
## [1] 146.6447 149.5127
c(AIC(res1),BIC(res1))
## [1] -62.71125 -56.97530
c(AIC(res2),BIC(res2))
## [1] -64.37179 -60.06983



실습 과제


  • iris 데이터에서, 종속변수 sl 을 가장 잘 설명하는 선형모형을 구축하시오

    • 어떤 모형이 가장 좋은 모형인지에 대하여 그 생각과 이유를 제시하고

    • 그 기준에 맞는 가장 좋은 모형을 찾으시오


  • iris 데이터에서, 종속변수 sl 이고, 독립변수로 sw, sp 를 사용하는 모형에서,

    • 교호작용이 있는 경우와 없는 경우의 모형을 구성하여 적합하고,

    • 그 결과를 그림으로 표현하고,

    • 어떤 모형이 더 좋은 모형이라고 생각하는지 그 근거를 설명하시오.


  • mtcars 데이터에서, 다음 모형들의 의미와 그 차이를 비교하고

    R 을 이용하여 적합한 결과 어떤 차이가 나타나는지 설명하시오

    • mpg ~ -1 + wt + factor(cyl)

    • mpg ~ wt + cyl

    • mpg ~ -1 + wt + cyl

    • mpg ~ -1 + wt + I(cyl-6)


  • mtcars 데이터에서, mpg 변수를 종속변수로 하는 가장 좋은 모형을 구성하려 한다.

    • gear 변수를 독립변수 중의 하나로 포함하는 모형

    • gear 변수를 ordered 변수로 변환하여 포함하는 모형

    • 위 두 가지 경우의 결과를 비교하여 설명하시오.


  • 케플러의 제3 법칙:

    • 태양계 행성 운행에 관한 케플러의 제3법칙을 회귀분석 관점에서 설명하시오

    • 태양계 행성 운행에 관한 자료를 R 로 분석하여, 케플러의 제3법칙을 확인하시오

    • 참고: [wiki], [YDLee]


  • trees 데이터에서,

    • \(log(v) = \beta_0 + \beta_1 \, \log(g) + \beta_2 \, log(h) +\epsilon,~\) 모형을 가정할 때

    • \(~ \beta_1 = 2 \, \beta_2\) 인지 가설검정하시오


  • trees 데이터에서,

    • 종속변수로 v 를 사용하는 모형 res0 와, log(v) 를 사용하는 res1 모형이 있다

    • 두 모형 res0 와 res1 을 비교하기 위하여 다음 방법을 적용하는 것의 타당성에 대하여 설명하시오

      • anova 에 의한 비교

      • AIC 에 의한 비교

      • BIC 에 의한 비교