[Tidy Modeling with R] 3. R Modeling Fundamentals

이번 포스팅에서는 3장의 R 모델링 기본원리에 대해 알아보기 위해 온도와 분당 귀뚜라미의 울음 비율 사이의 관계를 살펴보는 것으로 시작합니다.

우선 전체 데이터에 대해 loess 평활 회귀을 적합하여 추가한 그래프와 각 종별 데이터를 사용하여 적합한 선형회귀의 그래프를 그려보았습니다.

crickets %>% ggplot(mapping=aes(x=temp, y=rate)) + 
  geom_point() + geom_smooth(method = "loess")

crickets %>% ggplot(mapping=aes(x=temp, y=rate, color=species)) +
  geom_point() + geom_smooth(method="lm")

우선, data-driven 방법인 loess 평활 회귀에서 보면 울음 빈도와 온도는 상당히 선형적인 경향이 있다는 것을 알 수 있습니다.

또한 각 종별 데이터를 사용하여 적합한 선형회귀에서는 같은 온도일 때 종별로 우는 빈도에 차이가 있는 것처럼 보입니다.

 

R에서는 lm( ) 함수를 사용해 OLS 선형 회귀를 적합할 수 있습니다. 이 함수에서 중요한 파라미터는 formula이며 formula는 심볼릭(Symbolic)한데요.

예를 들면, rate ~ temp라는 formula가 있다고 했을 때, "~" 왼쪽에 있는 rate는 반응변수이며 "~" 오른쪽에 위치한 temp는 설명변수(예측변수)입니다.

이러한 방식으로 여러개의 설명변수를 추가할 수 있으며 설명변수로 범주형 변수를 추가하는 경우 (범주수-1)개의 더미변수로 변환됩니다.

 

또한, 상호작용항(Interaction)을 추가하기 위해 다음과 같은 3개의 formula를 적용할 수 있습니다.

• $y \sim x_1 + x_2 + x_1:x_2$
• $y \sim (x_1 + x_2)^2$
• $y \sim x_1*x_2$

 

추가적으로 formula는 표시변수(Indicate Variable)을 자동으로 생성할 뿐만 아니라 아래와 같은 편의성을 제공합니다.

• In-Line 함수를 formula에 사용할 수 있습니다. (Eg. $ y \sim log(x) $ )
• R은 formula 내부에 사용할 수 있는 다양한 함수를 제공하고 있습니다. (Eg. poly, ns function)
• 많은 변수가 있는 경우 dot(.)을 사용한 shortcut이 가능해 반응 변수를 제외한 모든 변수들을 지칭할 수 있습니다.

 

처음 예시로 돌아가서 이번에는 상호작용항을 추가하여 적합을 해보겠습니다.

interaction_fit1 <- lm(rate ~ (temp + species)^2, data = crickets)
interaction_fit2 <- lm(rate ~ temp * species, data = crickets)
interaction_fit3 <- lm(rate ~ temp + species + temp:species, data=crickets)
summary(interaction_fit1)
summary(interaction_fit2)
summary(interaction_fit3)

interaction fit1, fit2, fit3의 결과과 동일합니다.

• temp의 회귀계수의 유의확률이 0.05보다 작으므로 temp의 회귀계수는 0이라고 할 수 없습니다.
• species와 temp와 species의 상호작용항의 회귀계수의 유의확률이 0.05보다 크므로 회귀계수가 0이 아니라고 할 수 없습니다. 
• 위에서 언급한 것 처럼, 범주형 자료인 species는 2가지의 종이 있으므로 1개의 더미변수가 생깁니다. 

 

 

이제 상호작용항이 있는 모델과, 상호작용항이 없는 주효과만 있는 모델 중 어느것을 선택해야 할까요?

한가지 방법은 가능도비 검정(likelihood ratio test)을 이용하여 복잡한 모델과 단순한 모델의 가능도비를 계산하는 것으로 귀무가설은 "단순한 모델이 복잡한 모델처럼 데이터를 잘 설명한다" 입니다. 따라서 귀무가설을 채택하면 모델을 복잡한 모델에서 단순한 모델로 축약할 수 있으며, 귀무가설을 기각하면 복잡한 모델에서 단순한 모델로 축약할 수 있는 증거가 없으므로 복잡한 모델을 쓰는 것이 맞습니다. 다만 해당 내용은 통계파트에서 다루는것이 적절하다고 생각해서 자세한 내용은 생략하겠습니다.

R에서는 anova함수 또는 lrtest 함수를 사용해서 검정을 수행할 수 있습니다. 

main_effect_fit <- lm(rate~temp+species, data=crickets)

anova(interaction_fit1, main_effect_fit, test="LRT")
lmtest::lrtest(main_effect_fit, interaction_fit1)

(* 결과가 다른 이유는 https://stats.stackexchange.com/questions/155474/why-does-lrtest-not-match-anovatest-lrt   참고)

각각 0.2441, 0.2176의 유의확률을 가지므로 0.05 유의수준하에서 귀무가설을 채택할 수 있습니다. 따라서 단순한 모델인 교호작용항이 없는 모델이 복잡한 모델인 교호작용항이 있는 모델처럼 데이터를 잘 적합하므로 상호작용항이 모델에 유의하지 않음을 알 수 있습니다.

 

따라서, 상호작용항이 없는 주효과만 포함된 모델에서 회귀계수와 통계량을 summary 함수를 통해 살펴보면 아래와 같습니다.

summary(main_effect_fit)

• temp와 species의 유의확률이 0.05 이하이므로 분당 울음 비율과 두 변수과 선형적인 관계가 있다고 할 수 있습니다.
• 또한 temp의 회귀계수가 3.6으로 temp(온도)가 1만큼 증가하면 분당 울음 비율이 3.6회 증가하는 관계를 가집니다.

 

적합한 모델을 이용해서 예측을 진행할 수 있는데, R에는 predict( ) 함수를 사용해서 미지의 값을 예측할 수 있습니다.

예를들면, 종이 O. niveus이고 온도가 20이라고 해봅시다. 그러면 아래의 결과와 같이 분당 54.778 번의 비율로 운다는 것을 알 수 있습니다.

predict(main_effect_fit, data.frame(species="O. niveus", temp=20))

 

R formula는 많은 모델링 패키지에서 사용되며 일반적으로 아래와 같은 동작을 수행합니다.

• formula는 모델이 사용하는 열을 정의합니다.
• formula를 사용하여 열을 적절한 형식으로 변환(인코딩)합니다.
• 열의 역할은 formula에 의해 정의됩니다.

 

3 A Review of R Modeling Fundamentals | Tidy Modeling with R