9.6 모형 진단
9.6 모형 진단¶
회귀·모형적합에서 lm()·glm() 등으로 모형을 적합했다고 해서 끝이 아닙니다. "이 모형이 데이터를 얼마나 잘 설명하는가?", "변수를 하나 더 추가한 모형이 정말로 더 나은가?", "잔차에 아직 남아 있는 패턴은 없는가?"와 같은 질문에 답해야 비로소 그 모형을 신뢰하고 사용할 수 있습니다.
R에 내장된 trees 데이터셋(나무 31그루의 지름·높이·부피 실측치)으로 지름(Girth)·높이(Height)에서 부피(Volume)를 예측하는 회귀모형을 적합해 보겠습니다.
data(trees)
fit <- lm(Volume ~ Girth + Height, data = trees)
summary(fit)
#>
#> Call:
#> lm(formula = Volume ~ Girth + Height, 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 ***
#> Girth 4.7082 0.2643 17.816 < 2e-16 ***
#> Height 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(fit)을 실행하면 계수·표준오차·R²·F-통계량 등 다양한 정보가 텍스트로 한꺼번에 출력됩니다. 하지만 이 텍스트에서 특정 숫자(예: 회귀계수만, 잔차만)를 프로그램적으로 꺼내 쓰거나, 반복문으로 여러 후보 모형을 적합한 뒤 특정 지표만 모아 비교표를 만들려면 텍스트를 그대로 읽는 것만으로는 불편합니다.
이 절에서 다루는 함수들은 크게 네 가지 역할로 나뉩니다.
- 모형 객체(
lm,glm,aov등)에서 특정 구성요소를 그대로 꺼내는 추출 함수:coef(),fitted(),residuals() - 모형이 데이터에 얼마나 잘 들어맞는지를 숫자 하나로 요약하는 적합도·정보기준 함수:
deviance(),logLik(),AIC(),BIC() - 범주형 변수(그룹)에 따른 평균 차이를 검정하거나, 여러 모형을 통계적으로 비교하는 분산분석·모형비교 함수:
aov(),anova() - 잔차의 분포 형태를 살펴보는 분포 확인 함수:
density()
lm()·glm()·aov()처럼 모형을 적합하는 함수는 계수·잔차·적합값·자유도 등 수십 개의 요소를 담은 리스트(list) 형태의 S3 객체를 반환합니다. 이 절의 함수는 대부분 제네릭(generic) 함수여서, 같은 함수 이름이라도 넘겨받은 객체의 클래스(lm, glm, aov, nls 등)에 따라 내부적으로 서로 다른 메서드(residuals.lm(), residuals.glm() 등)가 호출되어 그 객체에 맞는 방식으로 계산됩니다. summary()(9.2절)가 벡터·데이터프레임·팩터마다 다르게 동작했던 것과 같은 원리입니다.
아래 표는 이 절에서 다루는 함수를 정리한 것입니다.
| 분류 | 함수 | 반환값 |
|---|---|---|
| 모형 요소 추출 | fitted() |
적합값(예측값) |
residuals() |
잔차(실제값 − 적합값) | |
coef() |
회귀계수(모수 추정치) | |
| 적합도·정보기준 | deviance() |
이탈도(선형회귀에서는 잔차제곱합) |
logLik() |
로그가능도 | |
AIC() |
아카이케 정보기준(작을수록 좋음) | |
BIC() |
베이즈 정보기준(작을수록 좋음, 변수 개수에 더 엄격) | |
| 분산분석·모형비교 | aov() |
분산분석 전용 모형 객체 |
anova() |
분산분석표 또는 중첩모형 비교표 | |
| 분포 확인 | density() |
커널 밀도추정 객체 |
적합값과 잔차 추출: fitted(), residuals()¶
모형이 예측한 값(적합값)과 실제 관측값의 차이(잔차)를 확인하지 않으면, 모형이 특정 구간에서만 잘 맞고 다른 구간에서는 체계적으로 틀리는 문제를 놓치기 쉽습니다.
trees 데이터로 적합한 모형에서 첫 6그루의 예측 부피와 실제 부피의 차이를 비교해 보겠습니다.
data(trees)
fit <- lm(Volume ~ Girth + Height, data = trees)
head(fitted(fit)) # 모형이 예측한 부피
#> 1 2 3 4 5 6
#> 4.837660 4.553852 4.816981 15.874115 19.869008 21.018327
head(trees$Volume) # 실제 부피
#> [1] 10.3 10.3 10.2 16.4 18.8 19.7
head(residuals(fit)) # 실제값 - 적합값
#> 1 2 3 4 5 6
#> 5.4623403 5.7461484 5.3830187 0.5258848 -1.0690084 -1.3183270
fitted()는 모형 수식에 원래 데이터를 대입해 얻은 예측값을, residuals()는 "실제값 − 적합값"을 반환합니다. 잔차가 0 근처에 무작위로 흩어져 있지 않고 특정 패턴(곡선 모양, 점점 커지는 폭 등)을 보인다면 선형성·등분산성 같은 회귀분석의 전제조건을 의심해 봐야 합니다.
i번째 관측치에 대해 적합값의 수식은 ŷᵢ = Xᵢβ̂(X는 설계행렬, β̂는 추정된 계수), 잔차 수식은 eᵢ = yᵢ − ŷᵢ입니다.
fitted()¶
fitted(object, ...)는 모형 객체 object가 각 관측치에 대해 예측한 적합값(fitted value)을 반환합니다. fitted.values()는 완전히 동일한 결과를 주는 별칭(alias)입니다.
object:lm(),glm(),aov(),nls()등으로 적합한 모형 객체...: 클래스별 메서드에 추가로 전달할 인자(대부분의 경우 사용하지 않음)
data(trees)
fit <- lm(Volume ~ Girth + Height, data = trees)
head(fitted(fit))
#> 1 2 3 4 5 6
#> 4.837660 4.553852 4.816981 15.874115 19.869008 21.018327
head(fitted.values(fit)) # fitted()와 완전히 동일
#> 1 2 3 4 5 6
#> 4.837660 4.553852 4.816981 15.874115 19.869008 21.018327
💡
fitted()와predict()는 어떻게 다른가요?fitted()는 모형을 적합할 때 사용했던 원래 데이터에 대한 예측값만 돌려줍니다. 모형에 없던 새로운 데이터에 대한 예측이 필요하다면predict()를 사용해야 합니다. 인자 없이predict(fit)을 호출하면fitted(fit)과 같은 결과가 나오는 것도 이 때문입니다.
residuals()¶
residuals(object, ...)는 모형 object의 잔차(residual)를 반환합니다. resid()는 완전히 동일한 결과를 주는 별칭입니다.
object: 적합된 모형 객체type: (glm등 일부 클래스에서) 잔차의 정의 방식을 지정. 예를 들어glm객체에는"deviance"(기본값)·"pearson"·"working"·"response"등 서로 다른 정의의 잔차가 있습니다....: 클래스별 메서드에 추가로 전달할 인자
data(trees)
fit <- lm(Volume ~ Girth + Height, data = trees)
head(residuals(fit))
#> 1 2 3 4 5 6
#> 5.4623403 5.7461484 5.3830187 0.5258848 -1.0690084 -1.3183270
head(resid(fit)) # residuals()와 완전히 동일
#> 1 2 3 4 5 6
#> 5.4623403 5.7461484 5.3830187 0.5258848 -1.0690084 -1.3183270
# 직접 계산으로 검증: 실제값 - 적합값
head(trees$Volume - fitted(fit))
#> 1 2 3 4 5 6
#> 5.4623403 5.7461484 5.3830187 0.5258848 -1.0690084 -1.3183270
회귀계수 추출: coef()¶
summary(fit)의 출력에서 계수 숫자만 골라 다른 계산에 쓰거나, 반복문으로 여러 모형을 적합한 뒤 계수만 모아 비교표를 만들어야 할 때, 텍스트 출력을 눈으로 보고 옮겨 적는 방식은 비효율적이고 오류가 나기 쉽습니다.
trees 데이터 모형에서 Girth(지름)의 계수만 뽑아 "지름이 1인치 늘어날 때 부피가 몇 세제곱피트 늘어나는가"를 확인해 보겠습니다.
data(trees)
fit <- lm(Volume ~ Girth + Height, data = trees)
coef(fit)["Girth"]
#> Girth
#> 4.708161
coef()는 모형 객체 안에 저장된 계수 벡터를 이름(변수명)이 붙은 숫자 벡터 형태로 그대로 꺼내 옵니다. 이렇게 이름이 붙어 있기 때문에 coef(fit)["Girth"]처럼 변수 이름으로 특정 계수에 바로 접근할 수 있습니다.
coef()¶
coef(object, ...)는 모형 object의 계수(coefficient)를 반환합니다. coefficients()는 완전히 동일한 결과를 주는 별칭입니다.
object: 적합된 모형 객체...: 클래스별 메서드에 추가로 전달할 인자
data(trees)
fit <- lm(Volume ~ Girth + Height, data = trees)
coef(fit)
#> (Intercept) Girth Height
#> -57.9876589 4.7081605 0.3392512
모형 적합도와 정보기준¶
변수를 하나 더 추가하면 잔차제곱합(RSS)은 거의 항상 줄어듭니다. 그렇다고 변수를 무한정 추가하는 것이 좋은 모형은 아니며(과적합), 단순히 RSS가 더 작다는 이유만으로 변수가 많은 모형이 항상 우월하다고 볼 수는 없습니다. 변수 개수가 다른 모형끼리 공정하게 비교하려면 "설명력"과 "복잡도(변수 개수)"를 함께 고려하는 지표가 필요합니다.
trees 데이터에서 Girth만 쓰는 단순회귀모형과, Girth·Height를 모두 쓰는 다중회귀모형을 비교해 보겠습니다.
data(trees)
fit1 <- lm(Volume ~ Girth, data = trees) # 변수 1개
fit2 <- lm(Volume ~ Girth + Height, data = trees) # 변수 2개
AIC(fit1, fit2)
#> df AIC
#> fit1 3 181.6447
#> fit2 4 176.9100
Height를 추가한 fit2의 AIC(176.91)가 fit1의 AIC(181.64)보다 작으므로, 변수 개수가 늘어난 것을 감안하더라도 fit2가 더 나은 모형이라고 판단할 수 있습니다.
deviance()는 모형이 설명하지 못하고 남긴 정도(선형회귀에서는 잔차제곱합과 같음)를, logLik()은 관측된 데이터가 이 모형에서 나왔을 가능성(가능도)을 로그 스케일로 나타냅니다. AIC()·BIC()는 로그가능도에 변수 개수(모형의 복잡도)에 대한 벌점을 더해, 변수 개수가 다른 모형끼리도 공정하게 비교할 수 있는 단일 지표로 만든 것입니다. 두 지표 모두 값이 작을수록 더 좋은 모형을 의미합니다.
로그가능도를 logL, 모형의 추정 모수 개수를 k, 관측치 개수를 n이라 하면 AIC = -2 × logL + 2k, BIC = -2 × logL + k × log(n)입니다. n이 커질수록 log(n)이 2보다 커지므로, BIC는 AIC보다 변수 개수가 많은(복잡한) 모형에 더 엄격한 벌점을 부과합니다.
deviance()¶
deviance(object, ...)는 모형 object의 이탈도(deviance)를 반환합니다. 선형회귀(lm)에서는 잔차제곱합(residual sum of squares, RSS)과 정확히 같은 값이며, 로지스틱 회귀 등 glm() 모형에서는 가능도 기반의 일반화된 적합도 부족 지표입니다.
object: 적합된 모형 객체...: 클래스별 메서드에 추가로 전달할 인자
data(trees)
fit <- lm(Volume ~ Girth + Height, data = trees)
deviance(fit)
#> [1] 421.9214
# 직접 계산으로 검증: 잔차제곱합
sum(residuals(fit)^2)
#> [1] 421.9214
logLik()¶
logLik(object, ...)는 모형 object의 로그가능도(log-likelihood)를 반환합니다.
object: 적합된 모형 객체...: 클래스별 메서드에 추가로 전달할 인자
data(trees)
fit <- lm(Volume ~ Girth + Height, data = trees)
logLik(fit)
#> 'log Lik.' -84.45499 (df=4)
결과 옆의 df=4는 추정된 모수의 개수(절편·Girth·Height 계수 3개 + 오차의 분산 1개)를 뜻하며, 이 값이 바로 AIC()·BIC() 계산에 쓰이는 k입니다.
AIC() / BIC()¶
AIC(object, ..., k = 2)는 아카이케 정보기준(Akaike Information Criterion)을, BIC(object, ...)는 베이즈 정보기준(Bayesian Information Criterion, AIC(object, k = log(n))과 동일)을 반환합니다.
object: 적합된 모형 객체(여러 개를 쉼표로 나열하면 한 번에 비교표를 얻을 수 있음)k: (AIC 전용) 모수 개수에 곱할 벌점 계수. 기본값 2가 일반적인 AIC이며,k = log(n)으로 지정하면 BIC와 동일한 결과를 얻습니다....: 비교할 모형을 추가로 나열하거나, 클래스별 메서드에 전달할 인자
data(trees)
fit1 <- lm(Volume ~ Girth, data = trees)
fit2 <- lm(Volume ~ Girth + Height, data = trees)
AIC(fit1); AIC(fit2)
#> [1] 181.6447
#> [1] 176.91
BIC(fit1); BIC(fit2)
#> [1] 185.9467
#> [1] 182.6459
💡 이 함수들은
lm()모형에만 적용되나요? 아닙니다.deviance()·logLik()·AIC()·BIC()는 모두 제네릭 함수이므로glm()으로 적합한 로지스틱 회귀 등에도 동일하게 적용됩니다. 다만glm()에서는deviance()가 잔차제곱합이 아니라 가능도 기반의 이탈도를 반환한다는 점이 다릅니다.data(mtcars) # 변속기 종류(am: 0=자동, 1=수동)를 마력(hp)·무게(wt)로 예측하는 로지스틱 회귀 glm_fit <- glm(am ~ hp + wt, data = mtcars, family = binomial) deviance(glm_fit) # 잔차이탈도(residual deviance) #> [1] 10.05911 glm_fit$null.deviance # 설명변수 없이 절편만 있는 모형의 이탈도 #> [1] 43.22973 AIC(glm_fit) #> [1] 16.05911절편만 있는 모형(
null.deviance, 43.23)보다 hp·wt를 포함한 모형의 이탈도(deviance, 10.06)가 훨씬 작으므로, 두 설명변수가 변속기 종류를 예측하는 데 실질적으로 기여하고 있다고 판단할 수 있습니다. 이 두 이탈도의 차이를 카이제곱분포로 검정하는 것이anova(glm_fit, test = "Chisq")이며, 기본적인 사고방식은 9.6.5절의anova()모형 비교와 동일합니다.⚠️ 주의할 점:
AIC()·BIC()는 같은 데이터, 같은 종속변수로 적합한 모형끼리 비교할 때만 의미가 있습니다. 종속변수를 변환(예:log(Volume)vsVolume)했거나 결측값 제거 등으로 실제 사용된 행 개수가 다른 모형끼리는 비교할 수 없습니다.
분산분석: aov()¶
두 집단의 평균 차이는 t.test()로 검정할 수 있지만, 세 집단 이상의 평균을 동시에 비교해야 할 때 두 집단씩 짝지어 t.test()를 여러 번 반복하면 검정을 반복할수록 우연히 유의한 결과가 나올 확률(1종 오류)이 누적되어 커집니다.
R에 내장된 PlantGrowth 데이터셋(식물 30그루를 대조군·처치군1·처치군2 세 그룹에 배정해 무게를 측정한 실험 데이터)으로, 그룹 간 무게 차이가 있는지 한 번에 검정해 보겠습니다.
data(PlantGrowth)
table(PlantGrowth$group) # 그룹별 표본 수
#> ctrl trt1 trt2
#> 10 10 10
aov_fit <- aov(weight ~ group, data = PlantGrowth)
summary(aov_fit)
#> Df Sum Sq Mean Sq F value Pr(>F)
#> group 2 3.766 1.8832 4.846 0.0159 *
#> Residuals 27 10.492 0.3886
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p-값(0.0159)이 0.05보다 작으므로, 세 그룹의 평균 무게가 모두 같다는 귀무가설을 기각하고 "적어도 한 그룹은 다른 그룹과 평균이 다르다"고 결론지을 수 있습니다.
aov()는 "분산분석(Analysis of Variance)"이라는 이름 그대로, 범주형 변수(그룹)가 연속형 반응변수의 평균에 영향을 주는지를 검정하기 위한 모형을 적합합니다. 내부적으로는 lm()과 동일한 최소제곱 계산을 사용하지만, 결과를 요약할 때 개별 계수보다 그룹 간 분산 비교(F-검정) 표 형태에 특화되어 있다는 점이 다릅니다.
aov()¶
aov(formula, data, ...)는 분산분석 모형을 적합합니다. 내부적으로 그룹 간 제곱합(between-group)과 그룹 내 제곱합(within-group)의 비율로 F-통계량을 계산합니다.
formula:반응변수 ~ 그룹변수형태의 모형 수식. 그룹변수가 여러 개면weight ~ group1 + group2(교호작용 없음) 또는weight ~ group1 * group2(교호작용 포함) 형태로 확장할 수 있습니다.data: formula에 사용된 변수가 들어 있는 데이터프레임...:lm()과 동일하게weights(가중치),subset(일부 행만 사용) 등을 지정할 수 있습니다.
data(PlantGrowth)
aov_fit <- aov(weight ~ group, data = PlantGrowth)
aov_fit # print() 결과 - 제곱합·자유도 요약
#> Call:
#> aov(formula = weight ~ group, data = PlantGrowth)
#>
#> Terms:
#> group Residuals
#> Sum of Squares 3.76634 10.49209
#> Deg. of Freedom 2 27
#>
#> Residual standard error: 0.6233746
#> Estimated effects may be unbalanced
aov()의 반환값도 결국 lm 계열 객체이므로, 앞서 다룬 coef()·fitted()·residuals()를 그대로 적용할 수 있습니다.
💡
coef(aov_fit)은 왜 그룹 이름을 그대로 보여주지 않을까요? 결과에group대신grouptrt1·grouptrt2라는 이름이 나오는 것은, 범주형 변수(팩터)를 회귀모형에 넣을 때 기준 그룹(여기서는 알파벳순으로 가장 앞선 "ctrl")을 절편(Intercept)에 흡수시키고 나머지 그룹은 기준 그룹과의 "차이"로 표현하기 때문입니다(가변수(dummy variable) 처리). 그래서 분산분석의 목적에는 계수 하나하나보다summary(aov_fit)처럼 그룹 전체에 대한 F-검정 결과를 보는 쪽이 더 적합합니다.
모형 비교: anova()¶
하나의 모형 안에서 어떤 변수가 반응변수의 분산을 얼마나 설명하는지 보고 싶을 때도 있고, 변수 하나를 추가하기 전과 후의 모형이 통계적으로 유의하게 다른지 직접 검정하고 싶을 때도 있습니다.
trees 데이터에서 Height를 추가하는 것이 통계적으로 의미가 있는지 anova()로 직접 검정해 보겠습니다.
data(trees)
fit1 <- lm(Volume ~ Girth, data = trees) # 축소모형
fit2 <- lm(Volume ~ Girth + Height, data = trees) # 완전모형
anova(fit1, fit2)
#> Analysis of Variance Table
#>
#> Model 1: Volume ~ Girth
#> Model 2: Volume ~ Girth + Height
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 29 524.30
#> 2 28 421.92 1 102.38 6.7943 0.01449 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p-값(0.01449)이 0.05보다 작으므로, Height를 추가한 완전모형이 축소모형보다 통계적으로 유의하게 더 잘 맞는다고 판단할 수 있습니다. 앞서 AIC(fit1, fit2) 비교에서 fit2가 더 낫다고 나왔던 결론과 방향이 일치합니다.
anova()는 넘겨받은 모형이 하나인지 여러 개인지에 따라 동작이 달라지는 제네릭 함수입니다. 모형을 하나만 넣으면 그 모형에 포함된 각 변수가 잔차제곱합을 얼마나 줄이는지 순서대로 분해한 분산분석표를 보여주고, 서로 중첩(nested)된 모형 두 개 이상을 넣으면 두 모형의 잔차제곱합 차이를 F-검정으로 비교합니다.
변수가 적은 모형(M1)과 변수가 많은 모형(M2)을 비교할 때, F = ((RSS₁ − RSS₂) / (df₁ − df₂)) / (RSS₂ / df₂)이며, 이 F 값이 F(df₁−df₂, df₂) 분포에서 얼마나 극단적인지로 p-값을 계산합니다. 단, 이 방식은 M1이 M2에 완전히 포함되는(중첩된) 관계일 때만 유효합니다.
anova()¶
anova(object, ...)는 모형 object의 분산분석표를 반환하거나, 여러 모형을 인자로 넘기면 중첩모형 간 F-검정 결과를 반환합니다.
object: 적합된 모형 객체. 이 하나만 넣으면 "단일모형 분산분석표", 중첩된 모형 여러 개를 쉼표로 이어서 넣으면 "모형 비교"로 동작합니다.test: (glm등 일부 클래스에서) 검정 방식을 지정합니다. 예:"Chisq"(카이제곱검정),"F"(F검정)...: 비교할 추가 모형을 나열하거나, 클래스별 메서드에 전달할 인자
data(trees)
fit2 <- lm(Volume ~ Girth + Height, data = trees)
# 인자 1개 - 단일모형 분산분석표: 각 변수가 설명하는 제곱합을 순서대로 분해
anova(fit2)
#> Analysis of Variance Table
#>
#> Response: Volume
#> Df Sum Sq Mean Sq F value Pr(>F)
#> Girth 1 7581.8 7581.8 503.1503 < 2e-16 ***
#> Height 1 102.4 102.4 6.7943 0.01449 *
#> Residuals 28 421.9 15.1
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
⚠️
anova(fit2)의 결과는 변수를 수식에 적은 순서에 따라 달라질 수 있습니다. 위 표는 Girth를 먼저 넣어 설명한 뒤, 남은 잔차에서 Height가 추가로 얼마나 설명하는지를 보여주는 순차적(sequential, Type I) 분해 방식입니다.Height + Girth순서로 수식을 바꾸면 각 행의 제곱합이 달라질 수 있습니다. 변수 순서와 무관하게 "다른 변수를 모두 통제한 상태에서 이 변수 하나만의 순수한 기여도"를 보고 싶다면summary(fit2)의 계수별 t-검정 결과를 참고하거나,car패키지의Anova(type = "II"/"III")(Base R 범위 밖) 같은 외부 도구를 검토해야 합니다.
잔차의 분포 확인: density()¶
선형회귀의 여러 통계적 추론(계수의 p-값 등)은 잔차가 정규분포를 따른다는 가정 위에 서 있습니다. 히스토그램(18.2절)으로도 분포 모양을 볼 수 있지만 막대 구간(bin) 폭에 따라 인상이 달라지는 단점이 있어, 더 부드러운 곡선으로 분포를 확인하고 싶을 때가 있습니다.
trees 데이터 모형의 잔차가 대략 종 모양(정규분포와 비슷한 모양)을 띠는지 확인해 보겠습니다.
data(trees)
fit <- lm(Volume ~ Girth + Height, data = trees)
r <- residuals(fit)
d <- density(r) # 커널 밀도추정 계산
d
#>
#> Call:
#> density.default(x = r)
#>
#> Data: r (31 obs.); Bandwidth 'bw' = 1.639
#>
#> x y
#> Min. :-11.323 Min. :8.849e-05
#> 1st Qu.: -5.142 1st Qu.:6.114e-03
#> Median : 1.039 Median :3.770e-02
#> Mean : 1.039 Mean :4.040e-02
#> 3rd Qu.: 7.220 3rd Qu.:6.831e-02
#> Max. : 13.402 Max. :1.018e-01
density()는 관측치 하나하나에 정규분포 모양의 작은 봉우리(커널)를 씌운 뒤 모두 합산해, 이산적인 데이터로부터 부드럽게 이어진 확률밀도 곡선을 추정합니다. 반환값은 숫자 벡터가 아니라 x좌표·y좌표 등을 담은 density 클래스 객체이며, 이 객체를 plot()에 넘기면 곡선을 그래프로 그릴 수 있습니다.
density()가 반환하는 리스트에는 다음 요소가 들어 있습니다.
x: 밀도를 계산한 지점들(기본값 512개)y: 각 x 지점에서의 추정 밀도값bw: 실제로 사용된 대역폭(bandwidth) — 커널 하나하나의 폭을 조절하는 값으로, 클수록 곡선이 더 뭉툭(과소적합)해지고 작을수록 더 들쭉날쭉(과적합)해집니다.n: 밀도 추정에 사용된 관측치 개수
density()¶
density(x, bw = "nrd0", kernel = "gaussian", n = 512, na.rm = FALSE, ...)는 데이터 x의 커널 밀도추정(kernel density estimation) 결과를 반환합니다.
x: 밀도를 추정할 숫자 벡터bw: 대역폭 또는 대역폭을 계산하는 방법. 기본값"nrd0"(Silverman의 경험적 규칙을 변형한 방식) 외에"SJ"(Sheather-Jones, 대체로 더 정교하다고 알려짐) 등을 지정하거나,0.5처럼 숫자를 직접 지정할 수도 있습니다.kernel: 커널의 모양. 기본값은"gaussian"(정규분포 모양)이며"epanechnikov"·"rectangular"등도 가능합니다.n: 밀도를 계산할 지점의 개수(기본값 512). 그래프를 부드럽게 그리기 위한 해상도 개념이며, 실제 관측치 개수와는 무관합니다.na.rm:TRUE이면 결측값 제외
data(trees)
fit <- lm(Volume ~ Girth + Height, data = trees)
r <- residuals(fit)
d <- density(r)
str(d)
#> List of 7
#> $ x : num [1:512] -11.3 -11.3 -11.2 -11.2 -11.1 ...
#> $ y : num [1:512] 0.000116 0.000128 0.00014 0.000153 0.000168 ...
#> $ bw : num 1.64
#> $ n : int 31
#> $ call : language density.default(x = r)
#> $ data.name: chr "r"
#> $ has.na : logi FALSE
#> - attr(*, "class")= chr "density"
아래 코드를 직접 실행하면 잔차의 밀도 곡선과, 같은 평균·표준편차를 갖는 정규분포 곡선을 겹쳐 그려서 잔차가 정규분포에서 얼마나 벗어나는지 눈으로 비교할 수 있습니다(18장 R 그래프 참고. 그래프 자체는 이 문서에서 직접 생성하지 않으니, 아래 코드를 R에서 실행해 확인해 보시기 바랍니다).
plot(density(r), main = "잔차의 커널 밀도추정", xlab = "잔차")
curve(dnorm(x, mean = mean(r), sd = sd(r)), add = TRUE, lty = 2, col = "blue")
legend("topright", legend = c("잔차 밀도", "정규분포"), lty = c(1, 2),
col = c("black", "blue"))
더 알아보기: 회귀진단의 확장 도구 — 이상치와 영향점 찾기
잔차(
residuals())를 그대로 살펴보는 것만으로는 "관측치 하나가 회귀선 전체를 얼마나 심하게 왜곡시키고 있는가"까지는 알기 어렵습니다. Base R의stats패키지는 이런 개별 관측치의 영향력을 진단하는 함수를 별도로 제공합니다.
hatvalues(model): 각 관측치의 레버리지(leverage, 지렛값)를 반환합니다. 레버리지가 크다는 것은 그 관측치의 설명변수(x) 값이 다른 관측치들과 동떨어져 있어, 회귀선의 기울기를 크게 좌우할 잠재력이 있다는 뜻입니다.rstandard(model): 잔차를 표준화(각 잔차를 그에 해당하는 표준오차로 나눔)한 표준화잔차를 반환합니다.rstudent(model): 해당 관측치를 제외하고 다시 적합한 모형을 기준으로 표준화한 스튜던트화잔차를 반환하며, 이상치 탐지에rstandard()보다 더 민감합니다.cooks.distance(model): 쿡의 거리(Cook's distance)를 반환합니다. 레버리지와 잔차 크기를 함께 반영해, 그 관측치 하나를 제거했을 때 회귀계수 전체가 얼마나 달라지는지를 하나의 숫자로 요약한 값입니다. 흔히4/n(n은 관측치 개수)이나 1을 넘는 값을 주의 깊게 살펴볼 기준으로 삼습니다.influence.measures(model): 위 진단값들(dfbetas,dffits,cov.ratio,cook.d,hat등)을 표 하나로 한꺼번에 정리해 반환합니다.data(trees) fit <- lm(Volume ~ Girth + Height, data = trees) # 쿡의 거리가 가장 큰 관측치 3개 sort(cooks.distance(fit), decreasing = TRUE)[1:3] #> 31 18 3 #> 0.6052326 0.1775359 0.167319231번째 나무의 쿡의 거리(0.61)가 다른 관측치보다 두드러지게 크므로, 이 나무 하나가 회귀계수 전체에 상당한 영향을 주고 있다고 의심해 볼 수 있습니다.
나아가
lm객체를plot()에 직접 넘기면, 이 진단값들을 활용한 4가지 표준 진단 그래프(잔차 대 적합값, 정규 Q-Q, 표준화잔차의 절대값 제곱근, 잔차 대 레버리지)를 한 번에 확인할 수 있습니다. 아래 코드 역시 그래프이므로 직접 실행해서 확인해 보시기 바랍니다.par(mfrow = c(2, 2)) # 그래프 4개를 2x2로 배치(18.6절 그래프 요소 제어 참고) plot(fit) par(mfrow = c(1, 1)) # 배치를 원래대로 되돌림
