3.1 연속형 (정량형) 변수

여기서는 연속형(continuous) 또는 정량형(quantitative) 형태의 자료 표현 또는 분석법을 설명한다.

3.1.1 기술 통계량(Descriptive Statistic)

평균(mean), 분산(variance), 표준편차(standard variation), 중앙값(median), 최소값(min), 최대값(max), 범위(range) 이 흔히 쓰이고, 좀 더 사용빈도는 낮지만, 4분위수 범위(interquartile range, IQR), 왜도(skewness), 첨도(kurtosis) 등이 사용된다. 이런 것들을 통칭해서 기술통계량(descriptive statistic)이라고 한다. 통계량(statisstic)의 복수형은 statistics(통계량들)이다.

단일 변수 분석에는 기본적으로 1) 간단한 수치적 계산과 2) 그래프와 도표에 의한 분석방법으로 나뉘어 지며, 간단한 수치적 계산에는 중심위치(location)의 측도와 산포(dispersion)의 측도로 나누어 진다. 중심위치의 측도에는 평균, 중앙값을 비롯하여, 기본적인 수치적 측도인 사분위수, 최대값, 최소값등이 있으며, 산포의 측도에는 분산, 표준편차, 사분위수 범위, 범위 등이 있다. 그래프나 도표에 의한 분석방법에는 선도표, 기둥그래프, 줄기와 잎그림, 상자그림, 산점도, 시계열(time-series) 그림 등이 있다.

R을 설치하면 많은 자료집합(dataset)이 포함되어 있다. 그 중에 ‘lh’ 라는 dataset이 있는데, ?lh 라는 명령으로 그 내용을 살펴볼 수 있다. 그 설명에 따르면 어느 한 여성에서 10분 간격으로 48회 채혈하여 luteinizing hormone을 측정한 것이라고 한다. 

?lh

R에서 descriptive statistics를 위해 가장 많이 사용하는 함수는 summary이다. 다음과 같이 하면 분산과 표준편차는 못보지만 평균과 5개의 기본적인 quantile number를 볼 수 있다. 

summary(lh)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.40    2.00    2.30    2.40    2.75    3.50

분산과 표준편차를 보기위해서는 다음과 같은 함수를 사용한다. 

var(lh)
[1] 0.3043
sd(lh)
[1] 0.5516

sasLM package에는 SAS PROC UNIVARIATE와 유사한 UNIV()함수가 있어서 상당히 많은 기술통계량을 보여준다. 

require(sasLM)
UNIV(lh)
          nAll            nNA        nFinite           Mean       Variance 
      48.00000        0.00000       48.00000        2.40000        0.30426 
            SD             CV            SEM LowerConfLimit UpperConfLimit 
       0.55159       22.98306        0.07962        2.23983        2.56017 
   TrimmedMean            Min             Q1         Median             Q3 
       2.39318        1.40000        2.00000        2.30000        2.75000 
           Max          Range            IQR            MAD       Skewness 
       3.50000        2.10000        0.80000        0.59304        0.29289 
    SkewnessSE       Kurtosis     KurtosisSE  GeometricMean    GeometricCV 
       0.34315       -0.69371        0.67440        2.33783       23.63751

psych라는 package의 describe라는 함수도 summary()보다는 많은 기술통계량을 보여준다.

R에 package가 설치되어 있지 않다면 다음과 같은 명령어로 먼저 설치해준다.

install.packages("psych")

이미 package가 설치되어 있다면, 다음과 같이 package를 load한 후 describe함수를 사용할 수 있다.

require(psych)
describe(lh)
   vars  n mean   sd median trimmed  mad min max range skew kurtosis   se
X1    1 48  2.4 0.55    2.3    2.39 0.59 1.4 3.5   2.1 0.27    -0.84 0.08

자세히 보면 skewness와 kurtosis의 값이 sasLM의 결과와 다른 것을 볼 수 있다. 이는 계산공식이 달라서 그러한데, sasLM은 SAS, SPSS와 동일한 공식을 사용한다.

 

3.1.2 Plotting

입력 이후의 자료 분석 첫 단계는 plotting이다.

plot 함수는 자료의 유형에 따라 다른 그림이 나올 수 있는데, 위 자료는 시계열 자료 형태이어서 거기에 적당하다고 하는 형태로 보여준다.

연속형 변수인 경우에는 boxplot이 많이 사용된다.

oPar = par(mfrow=c(1, 2))
plot(lh, main="plot(lh)", cex.axis=0.8, cex.lab=0.8, cex.main=0.9)
boxplot(lh, main="boxplot(lh)", cex.axis=0.8, cex.lab=0.8, cex.main=0.9)
Time series plot and box plot

Figure 3.1: Time series plot and box plot

Histogram도 연속형 변수를 그림으로 표현할 때 많이 사용한다.

여기에서 구간을 자르는 것이 내부 default algorithm에 의해 결정되는데, SAS, SPSS와 같은 다른 소프트웨어와 다를 수 있다. 만약 수동으로 지정해주고 싶으면 breaks option을 이용할 수 있다.

밀도를 구하는 density 함수 결과를 plot 함수에 입력으로 사용하면 hist와 비교할 만한 graph를 그릴 수 있다. 

oPar = par(mfrow=c(1, 2))
hist(lh, main="hist(lh)", cex.axis=0.8, cex.lab=0.8, cex.main=0.9)
plot(density(lh), main="plot(density(lh))", cex.axis=0.8, cex.lab=0.8, cex.main=0.9)
Histogram vs. density plot

Figure 3.2: Histogram vs. density plot

y scale을 잘 맞추어주면 hist 그림과 겹쳐서 그릴 수도 있다.

oPar = par(mfrow=c(1, 2))
hist(lh, probability=TRUE, cex.axis=0.8, cex.lab=0.8, cex.main=0.9)
lines(density(lh))
Overlay of histogram and density

Figure 3.3: Overlay of histogram and density

많이 사용되지는 않지만, 많은 교재에서 stem-leaf(줄기-잎) 그림을 설명하고 있는데, R에서는 다음과 같이 하면 볼 수 있다. 

stem(lh)

  The decimal point is 1 digit(s) to the left of the |

  14 | 000
  16 | 00
  18 | 000000
  20 | 000000
  22 | 00000000
  24 | 00000
  26 | 000000
  28 | 000
  30 | 000
  32 | 000
  34 | 000

 

stripchart는 유용성에 비해 덜 이용되고 있다. 

stripchart(lh, "stack", pch=16)
Stripchart or stack plot

Figure 3.4: Stripchart or stack plot

Default method가 “overplot”이어서 method를 “stack”이나 “jitter”로 바꾸는 것이 좋다. 여러가지 옵션을 이용하여 찍히는 기호(symbol), 색, 축 표시, 제목 등을 지정해 줄 수 있다.

 

ggplot2는 다양한 plot을 제공해 주고, 고품질의 output을 위해 시도해 볼 만하다.

다음은 유용성은 별로 없지만, 단일 변수일 때 violin plot하는 방법을 보여준다.

require(ggplot2)
ggplot(data.frame(lh), aes(x=1, y=as.numeric(lh))) + geom_violin()
Violin plot of single variable

Figure 3.5: Violin plot of single variable

3.1.3 Tabulation

연속형 변수인데, 구간을 나누어 빈도를 구하고(계수, counting) 싶을 때가 있다. 그럴 때는 cut 함수와 table 함수를 쓸 수 있다.

다음은 구간을 나눈 결과를 lh2에 저장하고, table 함수로 counting한 것이다. 최소값이 1.4이고, 최대값이 3.5이며, observation 수가 48개이므로 1부터 4까지 0.5간격으로 잘라보면 적당할 것 같다.

lh2 = cut(lh, seq(1, 4, by=0.5))
table(lh2)
lh2
(1,1.5] (1.5,2] (2,2.5] (2.5,3] (3,3.5] (3.5,4] 
      3      11      16      11       7       0

결과에서 경계값이 어느 쪽에 속하는지 괄호에 쓰인 기호를 보면 할 수 있다. 결과에서 ‘(’ 와 ‘)’는 경계값을 포함하지 않는다는 뜻이고,’[’ 와 ’]’는 경계값을 포함한다는 뜻이다. 따라서, (1.5, 2]는 1.5초과 2이하를 의미한다. 만약 경계값을 ’이상’과 ’미만’으로 하고 싶으면 cut 함수에서 right=FALSE라는 option을 사용하면 된다.

 

만약에 앞과 같이 미리 구간을 자르고 frequency를 구해 둔 경우에는 barplot을 사용할 수 있다.

barplot(table(lh2))
Barplot of tabulated data

Figure 3.6: Barplot of tabulated data

lh와 lh2는 정확히 순서대로 1:1 match 되므로 다음과 같이 하나의 data.frame으로 만들 수 있다.

lh3 = data.frame(lh, lh2)

또한 lh2는 이제 범주화된 자료이므로 다음에 나오는 범주형 자료 분석 기법을 사용할 수 있다.

3.1.4 통계적 추론 (Statistical Inference)

통계적 추론은 모수 추정(parameter estimation)과 가설 검정(hypothesis test)으로 나뉘어 지는데, 자세한 것은 뒷 장들에서 논의하고, 이 장에서는 기초적인 것만 다룬다.

sasLM package를 설치하였다면 다음과 같이 t 분포를 이용한 신뢰구간을 구할 수 있다. SEM은 Standard Error of Sample Mean으로서 흔히 SE라고 부른다. 하지만, SE의 정의는 standard deviation of a stistic (통계량의 표준편차)이므로 여기에서는 좀 더 엄밀성을 위해서 SEM으로 표현하였다. 

tsum0(lh3, "lh", c("Mean", "SD", "N", "SEM", "LCL", "UCL", "Skewness", "Kurtosis"))
    Mean       SD        N      SEM      LCL      UCL Skewness Kurtosis 
 2.40000  0.55159 48.00000  0.07962  2.23983  2.56017  0.29289 -0.69371 
attr(,"y")
[1] "lh"

tsum0의 3번째 인자(argument)는 요약(summarization, colligative) 함수들이다. tsum0는 아직 data.frame과 따옴표로 둘러싸인 column name을 첫 두 인자로 받기 때문에, tsum0(lh)와 같은 명령어는 받지 못한다.

사실 이것들은 UNIV() 에서도 모두 볼 수 있는 것이다.

위에서 skewness와 kurtosis가 psych::describe의 결과와 다른데, 그 이유는 sasLM package는 SAS 및 SPSS의 방식으로 계산하기 때문이다. (어느 한 가지만 옳은 것은 아니고, 통계에서는 여러가지 합당한 계산식이 있는 경우들이 있다.) psych package의 skew 함수와 kurtosi 함수가 이들을 계산하는데, type=2라는 옵션을 사용해야 SPSS와 같은 값을 산출해준다. R의 기본 함수인 IQR로 SPSS와 다른 결과를 보이는데, 비슷하게 나오려면 type=6 이라는 option을 사용해야 한다.

대부분의 다른 R package (moments 등)에서 계산하는 skewness와 kurtosis가 SAS, SPSS 등과 다르다. 그 이유도 다른 공식을 사용하기 때문인데, 만약 같은 것을 원한다면 확인해 보고 사용해야 한다.

LCL은 Lower Confidence Limit의 약어이고, UCL은 Upper Confidence Limit이다. 여기에서는 신뢰수준(confidence level)을 0.95로 고정한 것이다. LCL은 lower bound of confidence interval이라고도 할 수 있으므로 lower.ci, LCI, LB, LL이라고 약칭할 수 있고, 마찬가지로 UCL은 upper.ci, UCI, UB, UL로 약칭할 수도 있을 것이다.

다음과 같이 formula를 써서 같은 결과를 얻을 수 있다.

tsum(~lh, lh3, e=c("Mean", "SD", "N", "SEM", "LCL", "UCL", "Skewness", "Kurtosis"))

만약 앞의 용어들을 바꾸고 싶다면 다음과 같이 바꾸어서 할 수 있다.

SE = SEM
LL = LCL
UL = UCL
tsum0(lh3, "lh", c("Mean", "SD", "length", "SE", "LL", "UL"))

위의 방식은 사실 SE, LL, UL이라는 함수를 새로 만드는 것이다. 다만, 이전에 만들어둔 다른 함수들과 이름이 충돌하지 않는 것이 좋다. 관측값의 갯수를 알기 위해 R의 기본 함수로 length를 사용하는데, tsum은 n으로 이름을 바뀌어 줄력한다.

 

3.1.5 Single Group t-test

t 분포를 이용하여 단일 표본에 대한 가설 검정도 할 수 있고, 위와는 신뢰 수준이 다른 신뢰구간을 구할 수도 있다. 

t.test(lh, mu=2.6)

    One Sample t-test

data:  lh
t = -2.5, df = 47, p-value = 0.02
alternative hypothesis: true mean is not equal to 2.6
95 percent confidence interval:
 2.24 2.56
sample estimates:
mean of x 
      2.4

위의 경우에 무한모집단의 평균이 2.6이라는 귀무가설(즉, 2.6이 아니라는 대립가설)을 검정한 결과, p-value(유의확률)이 0.05보다 작아서 유의수준(significance level) 0.05에서 귀무가설을 기각한다. 즉, 2.6이 아니라고 통계적인 결론을 내린다. 이게 어떤 의학적 의미를 갖는가 하는 것은 별개의 문제이다. 통계분석을 하고 유의성 검정은 거의 기계적으로 할 수 있지만, 그것의 의미(meaning)를 찾는 것은 그렇지 않다. 이 차후 과정을 interpretation(번역이라기보다 해석)이라고 하고, 이것은 통계의 의미를 잘 아는 해당 분야의 전문가가 해야 한다. 해석의 한자가 두 가지 있다. 解析은 영어로 analysis (이것은 다시 分析으로 번역된다)이고 解釋은 interpretation (이것은 다시 飜譯으로 번역된다)에 가깝다.

Single group t-test는 모집단의 분포가 정규분포인 경우에 sample size와 상관없이 쓸 수 있으며, sample size가 30 이상인 경우에는 모집단의 분포에 상관없이 표본평균들의 분포는 정규분포를 따르므로(중심 극한 정리 central limit theorem) 평균에 대한 추론은 할 수 있게 된다. 따라서, 평균에 대한 추론이면서 sample size가 30 이상인 경우에 다음의 정규성 검정(normality test) 결과는 중요하지 않다.

 

3.1.6 정규성 검정 (Normality Test)

정규성을 검정하는 여러가지 방법이 있지만, R에서 가장 손쉽고 흔히 사용하는 방법은 Shapiro-Wilk test인데, 다음과 같이 한다. 

shapiro.test(lh)

    Shapiro-Wilk normality test

data:  lh
W = 0.97, p-value = 0.3
shapiro.test(log(lh))

    Shapiro-Wilk normality test

data:  log(lh)
W = 0.98, p-value = 0.5

p-value가 0.05보다 크면 정규성을 크게 벗어나지 않는다고 한다. 하지만, 이것이 정규분포를 따름을 입증한 것은 아니다. log변환 후의 값이 약간 더 정규분포에 가까운 것을 알 수 있다.

 

그림을 그리기 위해서는 Q-Q plot (quantile-quantile plot)을 사용하는데, 명령어는 다음과 같다. 

oPar = par(mfrow=c(1, 2))
qqnorm(lh, main="Original scale", cex.axis=0.8, cex.lab=0.8, cex.main=0.9)
qqline(lh)
qqnorm(log(lh), main="Log transformed data", cex.axis=0.8, cex.lab=0.8, cex.main=0.9)
qqline(log(lh))
QQ plot of original and log transformed data

Figure 3.7: QQ plot of original and log transformed data

 

그림으로도 log 변환 후가 좀 더 직선에 가까운 것을 볼 수 있다.

정규분포를 따르는 지는 skewness와 kurtosis가 0 근처임을 보기도 한다.

MASS에 있는 fitdistr 라는 함수도 어떤 분포를 따르는 지 추론하는데 유용하다. 

require(MASS)
fitdistr(lh, "normal")
    mean       sd   
  2.40000   0.54582 
 (0.07878) (0.05571)

위에서 표준편차의 추정량(0.54581743)이 sd 함수로 구한 값(0.5515934)과 다른데, 이는 sqrt(n - 1)로 나누지 않고 sqrt(n)으로 나누는 maximum likelihood(ML) 방법을 썼기 때문이다. ML에 의한 추정값은 MLE라고 하는데, 분산이나 표준편차의 경우에는 biased estimator이기 때문에 평소에는 많이 사용하지 않는다. 다음과 같이 변환해서 확인해 볼 수 있다.

sd(lh)*sqrt((length(lh) - 1)/length(lh)) # 0.5458174

주어진 자료로 다른 분포에 fitting해 볼 수도 있다. (가능한 분포 이름들은 ?명령어로 확인할 수 있다.)