4.3 모형식 표현
R에서 통계적 모형 표현 방법
- 지금까지 별다른 설명 없이
~가 들어간 수식표현을 특정함수(예:lm(),t.test(), 심지어 그래프 생성에 필요한 함수 등)의 인수로 사용함. - R은 (통계적) 모형을 표현하기 위해 formula 표현을 사용 \(\rightarrow\) 일반적으로
좌변 ~ 우변형태로 표시 - 보통은 특정 함수 내에서 호출되며 데이터에 포함되어 있는 변수를 평가하지 않고 해당 함수에서 해석할 수 있도록 변수값을 불러올 수 있음.
- formula는 “language” 객체의 일종이며 “formula” 클래스를 속성으로 갖는 평가되지 않은 표현식(unevaluated expression)
typeof(quote(x + 10)) # 객체의 형태가 "language"[1] "language"class(quote(x + 10)) # 객체의 클래스가 "call"[1] "call"- R에서 formula을 특정하는
~의 의미는 “즉시 평가(evaluate)하지 않고 이 코드의 의미를 전달(캡쳐)” \(\rightarrow\) 인용(quote) 연산자로 볼 수 있는 이유임
# 수식 표현
a <- y ~ x
b <- y ~ x + b
c <- ~ x + y + z
typeof(c); class(c); attributes(c)[1] "language"[1] "formula"$class
[1] "formula"
$.Environment
<environment: R_GlobalEnv>- 가장 기본적인 formula의 형태는 아래와 같음
반응변수(response variable) ~ 독립변수(independent variables)~는 “(좌변)은 (우변)의 함수로 나타낸 모형” 으로 해석됨.- 우변과 좌변 모두 일반적으로 여러 개의 변수들이 있을 수 있으며, 해당 변수의 추가는
+로 표시됨 - 좌변은 반응변수, 우변은 설명변수를 의미
일반적으로 좌변에 \(y\)로 표현되는 반응변수는 학문 분야에 따라 종속변수(dependent variable), 표적변수(target variable), 결과변수(outcome variable), 레이블(label, \(y\)가 범주형일 경우) 등으로 명칭되며, 우변에 \(y\)를 설명하기 위해 사용하는 변수(\(x\))를 마찬가지로 분야와 성격에 따라 독립변수(independent variable), 설명변수(exploratory variable), 예측변수(predictor variable), 위험 인자(risk factor), 공변량(covariate) 등으로 명칭된다.
- formula를 구성하는 항으로 벡터 객체를 직접 사용할 수 있으나, 데이터 프레임의 변수명을 formula의 항으로 구성할 수 있음.
- formula의 항으로 지정된 벡터 또는 변수들의 값은 formula를 호출해 사용할 때 까지 연결되지 않은 “언어”로써만 존재
- formula는 양면수식(two-sided formula, 좌변과 우변 모두에 하나 이상의 항이 존재)과 단면수식(one-sided formula, 우변에만 하나 이상의 항 존재)
set.seed(1)
x1 <- rnorm(100, 2, 4)
x2 <- runif(100, 2, 4)
x3 <- rpois(100, 3)
y <- 2*x1 + 3*x2 + 0.5*x3 + 5 + rnorm(100, 0, 4)
f1 <- y ~ x # y는 x의 함수
f2 <- y ~ x1 + x2 # y는 x1과 x2의 함수를 지칭하는 모형
typeof(f2); class(f2); attributes(f2)[1] "language"[1] "formula"$class
[1] "formula"
$.Environment
<environment: R_GlobalEnv>names(iris)[1] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width" "Species" f3 <- Species ~ Sepal.Length + Sepal.Width + Petal.Length
# 붓꽃의 종을 꽃받침 길이, 너비, 꽃잎의 길이에 대한 선형 함수로 표현
f4 <- ~ x1 + x2
f5 <- y ~ x1 + x2 + x3
length(f3); length(f4)[1] 3[1] 2f4[[1]]; f4[[2]]`~`x1 + x2f5[[1]]; f5[[2]]; f5[[3]]`~`yx1 + x2 + x3수식 사용 이유
- 변수 간의 관계를 알기 쉽게 표현
- 복잡한 변수의 관계를 표현해 함수 내에서 손쉽게 해당 항에 대응하는 데이터 값에 접근 가능
수식표현 방법
- 위에서 기술한 바와 같이
좌변항 ~ 우변항으로 표현 formula()또는as.formula()함수를 통해 텍스트를 formula 형태로 생성 가능
f6 <- "y ~ x1 + x2 + x3"
h <- as.formula(f6)
hy ~ x1 + x2 + x3h <- formula(f6)
hy ~ x1 + x2 + x3fs <- c(f1, f2, f6) # formula 객체를 concatenate 할 경우 list 객체
fl <- lapply(fs, as.formula)formula로 표현한 모형의 항에 대응하는 값으로 데이터 행렬 및 데이터 프레임 생성
model.frame(): formula 객체에 표현된 항에 대응하는 데이터 값으로 이루어진 데이터 프레임 반환model.matrix(): 디자인 행렬을 생성하는 함수
d1 <- model.frame(y ~ x1 + x2) # 벡터값을 데이터 프레임으로 반환
head(d1)# formula를 구성하고 있는 변수명에 대응하는 변수를 데이터 프레임에서 추출
f3Species ~ Sepal.Length + Sepal.Width + Petal.Lengthd2 <- model.frame(f3, iris)
head(d2)# model.matrix()에서는 디자인 행렬만 반환
# y 값은 포함하지 않고 우변에 해당하는 항에 대응하는 데이터 반환
X1 <- model.matrix(y ~ x1 + x2)
head(X1) (Intercept) x1 x2
1 1 -0.5058152 2.535016
2 1 2.7345733 2.437291
3 1 -1.3425144 3.033594
4 1 8.3811232 2.537901
5 1 3.3180311 2.362337
6 1 -1.2818735 3.037152formula 관련 함수
terms(): formula 객체에 포함되어 있는 항의 구조 파악all.vars(): formula에 포함되어 있는 변수명 추출update(): formula를 구성하는 항을 수정
terms(f2)y ~ x1 + x2
attr(,"variables")
list(y, x1, x2)
attr(,"factors")
x1 x2
y 0 0
x1 1 0
x2 0 1
attr(,"term.labels")
[1] "x1" "x2"
attr(,"order")
[1] 1 1
attr(,"intercept")
[1] 1
attr(,"response")
[1] 1
attr(,".Environment")
<environment: R_GlobalEnv>all.vars(f2)[1] "y" "x1" "x2"update(f2, ~ . + x3)y ~ x1 + x2 + x3formula에 사용되는 연산자와 의미
| Symbol | Example | Meaning |
|---|---|---|
+ |
X + Z |
변수 항 포함 |
- |
X + Z - X |
변수 항 제거 |
: |
X:W |
X와 W의 교호작용 |
%in% |
X%in%W |
지분(nesting) |
* |
X*Z |
X + Z + X:Z |
^ |
(X+W+Z)^3 |
3차 교호작용항까지 모든 항을 포함 |
I |
I(X^2) |
as is: X^2을 포함 |
1 |
X - 1 |
절편 제거 |
. |
Y ~ . |
모든 변수 포함(X + W + Z) |
일반 연산 시 A %in% B의 의미는 A가 B의 원소를 포함하는지에 대한 논리값을 반환해 주지만, formula에서 %in%은 중첩 또는 지분(nesting)을 내포함.
R의 리스트 객체는 중첩 및 지분 구조의 대표적 형태임. 예를 들어 리스트에 포함된 한 원소에 대응하는 데이터의 형태 및 값은 동일 리스트의 다른 원소에 대응한
데이터의 형태 및 값이 다름. 즉, 리스트 객체는 한 객체에 여러 형태의 데이터 구조를 가질 수 있고 이를 중접된 구조라고 함.
또한 실험계획법(experimental design)에서 지분 설계(nested design) 방법이 있는데, 두 개 요인(factor) A와 B가 존재할 때 A의 수준에 따른 B의 수준이 모두 다른 경우, 즉 교호작용이 존재하지 않는 형태의 실험설계법을 지칭함. 예를 들어 A사와 B사의 오랜지 주스 당도에 차이가 있는지를 알기 위해 각 회사에서 생산하고 있는 주스 3종을 랜덤하게 선택했다고 가정해 보자. 여기서 주 관심사는 두 회사에서 생산한 주스의 당도의 동질성이지 오랜지 주스 간 당도 차이는 관삼사항이 아님. 이 경우, 주 관심요인은 회사(C)이고, 요인 C는 회사 A, B라는 두 개의 수준(level)을 갖고 있음. 오랜지 주스(O)는 각 회사 별로 3개의 수준을 갖고 있는데, 각 회사에서 생산하는 오랜지 주스는 생산 공정에 차이가 있기 떄문에 각 회사에 지분되어 있음. 즉, 회사 A에서 생산한 오랜지 주스 O1, O2, O3은 회사 B에서 생산한 O1, O2, O3과 다름.
set.seed(10)
x <- runif(30)
z <- runif(30, 2, 3)
w <- sample(0:1, 30, replace = TRUE)
y <- 3 + 2.5*x + x^2 + 1.5*z + 0.5 * w + 2*w*x + rnorm(30, 0, 2)
head(model.matrix(~ x)) (Intercept) x
1 1 0.50747820
2 1 0.30676851
3 1 0.42690767
4 1 0.69310208
5 1 0.08513597
6 1 0.22543662head(model.matrix(~ x + z)) # x + z (Intercept) x z
1 1 0.50747820 2.535597
2 1 0.30676851 2.093088
3 1 0.42690767 2.169803
4 1 0.69310208 2.899832
5 1 0.08513597 2.422638
6 1 0.22543662 2.747746head(model.matrix(~ x + z - x)) # x + z에서 z 항 제거 (Intercept) z
1 1 2.535597
2 1 2.093088
3 1 2.169803
4 1 2.899832
5 1 2.422638
6 1 2.747746head(model.matrix(~ x:w)) # 교호작용항만 포함 (Intercept) x:w
1 1 0.5074782
2 1 0.3067685
3 1 0.0000000
4 1 0.0000000
5 1 0.0000000
6 1 0.0000000head(model.matrix(~ x*w)) # x + w + x:w (Intercept) x w x:w
1 1 0.50747820 1 0.5074782
2 1 0.30676851 1 0.3067685
3 1 0.42690767 0 0.0000000
4 1 0.69310208 0 0.0000000
5 1 0.08513597 0 0.0000000
6 1 0.22543662 0 0.0000000head(model.matrix(~ (x + z + w)^3)) # x + z + w + x:z + z:w + x:w + x:w:z (Intercept) x z w x:z x:w z:w x:z:w
1 1 0.50747820 2.535597 1 1.2867602 0.5074782 2.535597 1.2867602
2 1 0.30676851 2.093088 1 0.6420935 0.3067685 2.093088 0.6420935
3 1 0.42690767 2.169803 0 0.9263056 0.0000000 0.000000 0.0000000
4 1 0.69310208 2.899832 0 2.0098799 0.0000000 0.000000 0.0000000
5 1 0.08513597 2.422638 0 0.2062536 0.0000000 0.000000 0.0000000
6 1 0.22543662 2.747746 0 0.6194427 0.0000000 0.000000 0.0000000head(model.matrix(~ x + I(x^2))) # x + x^2 (Intercept) x I(x^2)
1 1 0.50747820 0.257534127
2 1 0.30676851 0.094106916
3 1 0.42690767 0.182250156
4 1 0.69310208 0.480390494
5 1 0.08513597 0.007248133
6 1 0.22543662 0.050821668# head(model.matrix(~ x + x^2))
head(model.matrix(~ x - 1)) x
1 0.50747820
2 0.30676851
3 0.42690767
4 0.69310208
5 0.08513597
6 0.22543662dat <- data.frame(y, x, w, z)
head(model.matrix(y ~ ., data = dat)) (Intercept) x w z
1 1 0.50747820 1 2.535597
2 1 0.30676851 1 2.093088
3 1 0.42690767 0 2.169803
4 1 0.69310208 0 2.899832
5 1 0.08513597 0 2.422638
6 1 0.22543662 0 2.747746head(model.matrix(y ~ .^2, data = dat)) (Intercept) x w z x:w x:z w:z
1 1 0.50747820 1 2.535597 0.5074782 1.2867602 2.535597
2 1 0.30676851 1 2.093088 0.3067685 0.6420935 2.093088
3 1 0.42690767 0 2.169803 0.0000000 0.9263056 0.000000
4 1 0.69310208 0 2.899832 0.0000000 2.0098799 0.000000
5 1 0.08513597 0 2.422638 0.0000000 0.2062536 0.000000
6 1 0.22543662 0 2.747746 0.0000000 0.6194427 0.000000# nested 구조
dat2 <- expand.grid(C = c("A", "B"), O = paste0("O", 1:3),
y = runif(3, 1, 2))
dat2 <- dat2[order(dat2$C, dat2$O), ]
model.matrix(y ~ C + O %in% C, data = dat2) (Intercept) CB CA:OO2 CB:OO2 CA:OO3 CB:OO3
1 1 0 0 0 0 0
7 1 0 0 0 0 0
13 1 0 0 0 0 0
3 1 0 1 0 0 0
9 1 0 1 0 0 0
15 1 0 1 0 0 0
5 1 0 0 0 1 0
11 1 0 0 0 1 0
17 1 0 0 0 1 0
2 1 1 0 0 0 0
8 1 1 0 0 0 0
14 1 1 0 0 0 0
4 1 1 0 1 0 0
10 1 1 0 1 0 0
16 1 1 0 1 0 0
6 1 1 0 0 0 1
12 1 1 0 0 0 1
18 1 1 0 0 0 1
attr(,"assign")
[1] 0 1 2 2 2 2
attr(,"contrasts")
attr(,"contrasts")$C
[1] "contr.treatment"
attr(,"contrasts")$O
[1] "contr.treatment"