반응형

[이론] 콕스 회귀 모델에서 생존 함수의 표현

 

콕스 회귀 모델에서 생존 함수는 다음과 같이 표현하곤 한다. 

$$S(t,\mathbf{X}) = \left[ S_0 (t) \right] ^{e^{\sum_{i=1}^{p} \beta_i X_i}} $$

 

왜 이렇게 표현되는지를 이해해야 추후에 나오는 log-log plot 등을 이해할 수 있기에 여기에서 소개한다.

 

1. Notation

먼저 notation에 대해 소개하고자 한다.

1)$\mathbf{X}$

$\mathbf{X}$는 $X$들의 집합체이며, $X$는 모델을 설명하는 변수를 의미한다. 예를 들면, 특정 위험 요소에 노출 여부, 성별, 나이, 인종, 음주 여부, 흡연 여부 등이 있을 수 있다. 이 모든 것들을 적기에는 귀찮고 공간 낭비이므로 $\mathbf{X}$로 표현한다.

 

2)$S(t,\mathbf{X})$

$S(t,\mathbf{X})$는 특정 $\mathbf{X}$을 가진 사람의 시간 $t$에서의 생존 확률이다. 예를 들면, 남자, 53세, 백인, 음주자, 비흡연자의 12개월에서의 생존 확률인 것이다.

 

3) $S_0(t)$

$S_0 (t)$는 모델을 설명하는 변수들의 값이 기본값(reference)인 사람의 시간 $t$에서의 생존 확률을 의미한다. 예를 들어 모델의 기본값을 남성, 0세, 흑인, 비음주자, 비흡연자로 잡았다면, 이런 사람의 12개월에서의 생존 확률을 의미한다.

 

4) $e^{\sum_{i=1}^{p} \beta_i X_i}$

이 식에서 모델을 설명하는 변수는 $p$개 였음을 알 수 있고, $\beta_i$는 각 변수에 대한 회귀계수를 의미한다.

 

 

2. 생존 함수와 위험 함수의 관계

이전 포스팅에서 다음과 같은 관계를 얻어냈다. 2022.11.10 - [통계 이론] - [이론] 생존 함수와 위험 함수의 관계

$$S(t) = \exp \left[ - \int_{0}^{t} h(u) du  \right] = e^{- \int_{0}^{t} h(u) du}$$

 

 

3. 콕스 회귀 모델에서 생존 함수의 표현

추가적인 notation을 두 개만 다루고 넘아가겠다.

1)$h_{\mathbf{X}}(t)$

특정 $\mathbf{X}$을 가진 사람의 시간 $t$에서의 위험 함수를 나타낸다.

 

2)1)$h_0(t)$

모델을 설명하는 변수들의 값이 기본값(reference)인 사람의 시간 $t$에서의 위험 함수를 나타낸다.

 

그렇다면 다음을 알 수 있다.

 

$$\begin{align} S(t,\mathbf{X})&=\exp\left[-\int_{0}^{t}h_{\mathbf{X}}(u)du\right]\\&=\exp\left[-\int_{0}^{t}h_{0}(u)e^{\sum_{i=1}^{p}\beta_{i}X_i}du\right]\\&=\exp\left[-\int_{0}^{t}h_{0}(u)du\times e^{\sum_{i=1}^{p}\beta_{i}X_i}\right] \\&=\left[S_{0}(t)\right]^{\exp\left(\sum_{i=1}^{p}\beta_{i}X_i \right)} \\&& \end{align}$$

 

 

 

[이론] 콕스 회귀 모델에서 생존 함수의 표현 정복 완료!

 

작성일: 2022.11.10.

최종 수정일: 2022.11.10.

 

 

 

반응형
반응형

생존 함수와 위험 함수의 관계

 

 생존 분석에 대해 공부를 하다 보면 생존 함수, 위험 함수가 나오게 되고, 그들의 관계를 나타내는 식이 등장한다.

 

1. 생존 함수와 위험 함수의 정의

먼저 생존 함수($S(t)$), 위험 함수($h(t)$)에 대해 설명하고 넘어가고자 한다.

 

1) 생존 함수

생존 함수 $S(t)$는 시간 $t$일 때, 살아있는 사람의 분율 (proportion)이다. 즉, 시간 $t$까지 살아남을 확률을 의미한다.

- 처음 ($t=0$)에는 모든 사람이 살아있으므로 $S(0)=1$이다.

- 시간이 무한히 흐르고 나면 모든 사람이 죽으므로 $S(\infty)=0$이다.

 

2) 위험 함수

위험 함수 $h(t)$는 그 개념이 조금 더 복잡하다. 거칠게 설명하면 다음과 같다.

"시간 $t$까지 살아남았을 때, 그 순간에 단위시간당 사망할 확률"

이 설명은 이해하기 어려우므로 예시를 들어 이해해보고자 한다.

 

(0) 상황

아침 9시에 200명으로 연구를 시작했다.

하지만 12시가 되었을 때 100명이 살아남았다.

 

(1) 1분간 관찰해보기

12시부터 12시 1분까지 10명이 죽었다고 하자. 즉, 사망할 확률은 10%(소수점으로 나타내면 0.1)이다

그렇다면, 60초동안 사망할 확률은 0.1이므로 위험 함수의 값은 다음과 같이 정해진다.

$$h(12시)= \frac {0.1} {60} = 0.00166667$$

 

그런데 위에서 위험 함수는 "그 순간에 단위시간당 사망할 확률"이라고 하였으므로 60초라는 간격은 너무 길다. 60초가 아니라 10초 동안 관찰해보자

 

(2) 10초간 관찰해보기

12시부터 12시 0분 10초까지 5명이 죽었다고 하자. 즉, 사망할 확률은 5%(소수점으로 나타내면 0.05)이다

그렇다면, 10초동안 사망할 확률은 0.05이므로 위험 함수의 값은 다음과 같이 정해진다.

$$h(12시)= \frac {0.05} {10} = 0.005$$

 

이렇게 시간 간격을 무한히 줄이기 위해 0에 가까워지는 극한을 적용하면 위험 함수의 값을 구할 수 있다.

$t$이상 살 확률을 시각 변수 $T$에 대해 $P(T\geq t)$라고 표현한다면, 위험 함수 $h(t)$는 다음과 같이 나타낼 수 있다.

$$h(t) = \lim_{\Delta t \rightarrow 0} \frac {P (t \leq T \leq t+ \Delta t \vert T \geq t)} {\Delta t} $$

 

2. 생존 함수와 위험 함수의 관계

보통 다음과 같은 관계로 나타난다.

$$h(t) = -\frac{S'(t)} {S(t)} = - \frac {dS(t)/dt} {S(t)}$$

또한 이는 일종의 미분방정식이므로 양변을 적분하고 정리하면 다음을 얻을 수 있다.

$$S(t) = \exp \left[ - \int_{0}^{t} h(u) du  \right] = e^{- \int_{0}^{t} h(u) du}$$

 

어떻게 이런 관계가 도출되는지 설명하고자 한다.

 

시간 $[0,1]$ 사이를 굉장히 큰 수 $n$에 대해 $n$개의 구간으로 나눈다.

처음 ($t=0$)에는 모든 사람이 살아있으므로 $S(0)=1$이다. 이번에는 $S \left( \frac {1} {n} \right)$을 구해볼 것이다.

그렇다면 우리가 관심을 갖는 시간의 구간은 $\left[ 0, \frac{1}{n} \right]$이다.

 

 위험 함수는 "그 순간에 단위시간당 사망할 확률"이므로 매 순간마다 다른 값을 가지는 것이 당연하지만, $n$이 매우 큰 수이므로 구간 $\left[ 0, \frac{1}{n} \right]$은 매우 짧은 찰나의 순간일 것이고, 그 사이에는 값이 거의 변하지 않는다고 봐도 무방하다. 이 구간의 $h(t)$는 $h \left( \frac {1} {n} \right)$으로 고정되어 있다고 가정하자.

 그런데 위험함수는 단위 시간당 사망할 확률이므로 기준이 구간 $[0,1]$이다. 그런데 우리가 관심을 갖는 구간은 $\left[ 0, \frac{1}{n} \right]$이므로 이 구간에서 사망할 확률은 $h \left( \frac {1} {n} \right)$에 $\frac{1}{n}$을 곱해주어야 한다.

$$\left[ 0, \frac{1}{n} \right]에~사망할~확률=h \left( \frac {1} {n} \right)\frac{1}{n}$$

이 구간에서 생존할 확률은 이 값을 1에서 빼주면 된다.

 

$$\left[ 0, \frac{1}{n} \right]에~생존할~확률=1-h \left( \frac {1} {n} \right)\frac{1}{n}$$

 

그렇다면 $S \left( \frac {1} {n} \right)$은 처음 ($t=0$)에 생존한 사람의 분율인 $S(0)$에 $\left( 1-h \left( \frac {1} {n} \right)\frac{1}{n} \right)$을 곱하여 구할 수 있다.

 

$$S \left( \frac {1} {n} \right) =S(0) \times \left( 1-h \left( \frac {1} {n} \right)\frac{1}{n} \right)=1-h \left( \frac {1} {n} \right)\frac{1}{n} $$

 

$\left[ 0, \frac{1}{n} \right]$에서의 생존할 확률을 구한 것과 같은 논리로 이를 일반화하면 $\left[ \frac{j-1}{n}, \frac{j}{n} \right]$에 생존할 확률은 다음과 같다.

$$ \left[ \frac{j-1}{n}, \frac{j}{n} \right]에~생존할~확률=1-h \left( \frac {j} {n} \right)\frac{1}{n} $$

 

한편 어떤 수 $k$에 대해  $\frac {k} {n}$ 까지 살아남을 확률을 구하는 것은

$\left[0,\frac{1}{n}\right]$에 생존한 사람이

$\left[\frac{1}{n},\frac{2}{n}\right]$에 생존한 사람이

$\left[\frac{2}{n},\frac{3}{n}\right]$에 생존한 사람이

$...$

$\left[\frac{k-1}{n},\frac{k}{n}\right]$에 생존할 확률을 구하는 것과 같다.

 

따라서 $\frac {k} {n}$ 까지 살아남을 확률인 $S \left( \frac{k}{n}\right)$은 다음과 같이 구할 수 있다.

$$\begin{align}S \left( \frac{k}{n}\right) &= S(0) \times\left(1-h \left( \frac {1} {n} \right)\frac{1}{n} \right)\times\left(1-h \left( \frac {2} {n} \right)\frac{1}{n} \right)\cdots\times\left(1-h \left( \frac {k} {n} \right)\frac{1}{n} \right)\\&=\prod_{j=1}^{k} \left( 1-h \left( \frac{j}{n} \right) \frac{1}{n} \right)\\&& \end{align}$$

 

그러므로 다음 등식이 성립한다.

$$S \left( \frac{k}{n} \right) \div S \left( \frac{k-1}{n} \right) = 1-h \left( \frac{k}{n} \right) \frac{1}{n}$$

 

이 식을 정리하면 다음을 얻게 되고

$$h \left( \frac{k}{n} \right) = - \frac{S \left( \frac{k}{n} \right) - S \left( \frac{k-1}{n} \right)}{ \frac{1}{n} S \left( \frac{k}{n} \right)}$$

 

$\lim_{ n \rightarrow \infty, \frac{k}{n} \rightarrow t}$의 극한을 씌우게 되면 다음 등식이 성립한다.

$$h(t) = -\frac{S'(t)} {S(t)} = - \frac {dS(t)/dt} {S(t)}$$

 

따라서 위 수식을 얻게 된다.

 

 

 

 

[이론] 생존 함수와 위험 함수의 관계 정복 완료!

 

작성일: 2022.11.10.

최종 수정일: 2022.11.10.

 

 

반응형
반응형

[R] 피셔 정확 검정에서 workspace 부족 에러 해결 방법

 

피셔 정확 검정 (Fisher's exact test)을 R로 구현하다 보면 다음과 같은 에러가 뜰 때가 있다.

 

1)

  FEXACT error 6.  LDKEY=592 is too small for this problem,   (ii := key2[itp=867] = 369304672, ldstp=17760)

Try increasing the size of the workspace and possibly 'mult'

 

2)

 FEXACT error 40. Out of workspace.

 

이럴 때에는 workspace를 늘려주어야 분석이 가능해진다.

Fisher's exact test를 시행하는 R 코드에서 다음과 같이 (workspace=2e8) 혹은 (workspace=2e9), ...과 같은 argument를 추가해주면 구동될 '수도' 있다. Fisher's exact test를 시행하는 R코드는 다음 링크에서 확인할 수 있다. 2022.09.02 - [범주형 자료 분석/R] - [R] 피셔 정확 검정 - fisher.test()

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.10.11)

 

분석용 데이터 (update 22.10.11)

2022년 10월 11일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 범주형 자료 분석 - 모평균 검정 - 반복 측정 자료 분석 - 통계

medistat.tistory.com

 

 

코드를 보여드리기에 앞서 워킹 디렉토리부터 지정하고 데이터를 불러와 df에 객체로 저장하겠다.

워킹 디렉토리에 관한 설명은 다음 링크된 포스트에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 작업 디렉토리 (Working Directory) 지정 - getwd(), setwd()

 

데이터 불러오는 방법은 다음 링크에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : EXCEL - read_excel(), read.xlsx()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 저장하기 : CSV 파일 - write.csv(), write_csv()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : SAS file (.sas7bdat) - read.sas7bdat(), read_sas()

 

 

(workspace=2e8) 혹은 (workspace=2e9), ...과 같은 argument를 추가하여 코드를 실행해보겠다.

#작업 디렉토리 지정
setwd("C:/Users/user/Documents/Tistory_blog")

#데이터 불러오기
install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

#피셔 정확 검정 시행하기
fisher.test(df$SEX, df$RH, workspace=2e8)

 

물론 우리의 데이터는 그렇게 크지 않아 workspace가 부족하지 않으며, workspace를 늘려줄 필요도 없다. 하지만 만약 여러분의 데이터로 피셔 정확 검정을 시행하던 차에 저런 문제가 생기는 경우 workspace를 늘려보기를 권한다.

 

[R] 피셔 정확 검정 workspace 에러 해결 방법 정복 완료!

 

작성일: 2022.11.10.

최종 수정일: 2022.11.10.

이용 프로그램: R 4.1.3

RStudio v1.4.1717

RStudio 2021.09.1+372 "Ghost Orchid" Release 

운영체제: Windows 10, Mac OS 10.15.7

반응형
반응형

[R] 일표본 T검정 (One-sample T-test) - t.test()

  세상에 존재하는 모든 사람을 대상으로 연구를 하면 참 좋겠지만, 시간적인 이유로, 그리고 경제적인 이유로 일부를 뽑아서 연구를 진행할 수밖에 없다. 모든 사람을 모집단이라고 하고, 뽑힌 일부를 표본이라고 한다. 우리는 표본으로 시행한 연구로 모집단에 대한 결론을 도출해내고자 할 것이다. 

 1000명에게 피검사를 시행하였고, 간 기능 검사의 일환으로 ALT 수치를 모았다. 이 데이터를 기반으로 1000명이 기원한 모집단 인구에서의 ALT평균이 어떻게 될지 예측하는 것이 T-test이다. T-test는 크게 세 가지로 나눌 수 있다.

 

 1) 일표본 T검정 (One sample T-test)

 : 모집단의 평균이 특정 값이라고 할 수 있는가?

 예) 모집단의 ALT 평균이 50이라고 할 수 있는가?

 

 2) 독립 표본 T검정 (Independent samples T-test): 2022.11.12 - [모평균 검정/R] - [R] 독립 표본 T검정 (Independent samples T-test) - t.test(), var.test(), levene.test()

 

 : 두 모집단의 평균이 다르다고 할 수 있는가?

 예) 고혈압 환자와 일반인의 수축기 혈압 평균이 서로 다르다고 할 수 있는가?

 

 3) 대응표본 T검정 (Paired samples T-test): 2022.11.25 - [반복 측정 자료 분석/R] - [R] 대응 표본 T검정 (Paired samples T-test) - t.test()

 : 모집단의 짝지어진 변수들의 평균이 다르다고 할 수 있는가?

예) 간기능 개선제 복용 전 ALT 평균은 간기는 개선제 복용 후 ALT 평균과 다르다고 할 수 있는가?

 

이번 포스팅에서는 일표본 T검정 (One sample T-test)에 대해 알아볼 것이다.

 

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.10.11)

 

분석용 데이터 (update 22.10.11)

2022년 10월 11일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 범주형 자료 분석 - 모평균 검정 - 반복 측정 자료 분석 - 통계

medistat.tistory.com

 

코드를 보여드리기에 앞서 워킹 디렉토리부터 지정하겠다.

워킹 디렉토리에 관한 설명은 다음 링크된 포스트에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 작업 디렉토리 (Working Directory) 지정 - getwd(), setwd()

 

setwd("C:/Users/user/Documents/Tistory_blog")

 

데이터를 불러와 df에 객체로 저장하겠다.

데이터 불러오는 방법은 다음 링크에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : EXCEL - read_excel(), read.xlsx()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 저장하기 : CSV 파일 - write.csv(), write_csv()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : SAS file (.sas7bdat) - read.sas7bdat(), read_sas()

 

install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

 

목표: 모집단의 ALT 평균이 50이라고 할 수 있는가?

 이번 포스팅의 목적은 1000명의 데이터를 가지고, 이 1000명이 기원한 모집단의 ALT평균이 50이라고 할 수 있는지 판단하는 것이다.

 

전제조건 (정규성) - 코드

일표본 (One sample) T 검정의 전제조건은 해당 분포가 정규성을 따른다는 것이다. 정규성 검정을 하도록 하겠다. 정규성 검정에 관한 분석 내용은 다음 글에서 살펴볼 수 있다. 

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (1) : Q-Q plot - qqnorm(), qqline()

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (2) : 히스토그램 - hist(), dnorm()

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (3) : 정량적 검정 (Shapiro-Wilk, Kolmogorov-Smirnov) - shapiro.test(), ks.test()

2022.08.12 - [기술 통계/R] - [R] 정규성 검정 (4) : 정량적 검정 (Lilliefors test) - lillie.test()

2022.08.16 - [기술 통계/R] - [R] 고급 Q-Q Plot - Van der Waerden, Rankit, Tukey, Blom

 

# 1) Q-Q plot 그리기
qqnorm(df$ALT)
qqline(df$ALT)

# 2) 히스토그램 그리기
hist(df$ALT, prob=TRUE)
ALTrange<-seq(min(df$ALT),max(df$ALT),length=max(max(df$ALT)-min(df$ALT),100))
ND<-dnorm(ALTrange,mean=mean(df$ALT),sd=sd(df$ALT))
lines(ALTrange, ND, lwd=2)

# 3) Shapiro-Wilk test 시행하기
shapiro.test(df$ALT)

 

 

정규성 검정 - 결과

 

 

	Shapiro-Wilk normality test

data:  df$ALT
W = 0.9982, p-value = 0.3749

 

 N수가 2,000개 미만이므로 Shapiro-Wilk 통계량의 p-value를 보면 0.05 이상이며, Q-Q Plot은 대부분의 데이터가 선상에 있고, 히스토그램에서도 정규성을 따르는 것처럼 볼 수 있다. 따라서 ALT로는 일표본 T검정 (One-sample T-test)를 시행할 수 있다. 

 

 

일표본 T검정 (One sample T test)의 귀무가설과 대립 가설

 

이번 일표본 T검정 (One sample T test)의 귀무가설은 하나다.

귀무가설: $H_0=$ 모집단의 ALT 평균은 50이다.

 

하지만 대립 가설은 3개가 된다.

1) $H_1=$ 모집단의 평균은 50보다 크다. (단측 검정)

2) $H_1=$ 모집단의 평균은 50보다 작다. (단측 검정)

3) $H_1=$ 모집단의 평균은 50보다 크거나 작다. (양측 검정)

 

대립 가설 별로 코드를 보고자 한다.

 

1)  $H_1=$ 모집단의 평균은 50보다 크다. (단측 검정)

t.test(df$ALT, mu=50, alternative="greater")

t.test(df$ALTmu=50alternative="greater") : T-test를 시행한다. 검정할 변수는 df 데이터의 ALT변수다. ALT의 모평균이 50이라고 할 수 있는지 검정할 것이다. 대립 가설은 "모평균이 50보다 크다"이다.

결과

	One Sample t-test

data:  df$ALT
t = -94.014, df = 999, p-value = 1
alternative hypothesis: true mean is greater than 50
95 percent confidence interval:
 34.85529      Inf
sample estimates:
mean of x 
 35.11594

위의 결과 중 다음 내용에 주목해야 한다.

One Sample t-test
data:  df$ALT
t = -94.014, df = 999, p-value = 1
alternative hypothesis: true mean is greater than 50
95 percent confidence interval:
 34.85529      Inf
 sample estimates:
mean of x  35.11594

 

P-value1.0000이다. 즉, 대립 가설 (모평균의 참값은 50보다 크다.)을 선택하는 것이 말이 안 된다는 것이다. 평균이 50보다 크다고 보는 것보다는 50이라고 보는 것이 더 합리적이라는 뜻이다. 실제로 평균은 35.1159라고 적혀있으므로 합당한 통계분석임을 알 수 있다.

 

2)  $H_1=$ 모집단의 평균은 50보다 작다. (단측 검정)

t.test(df$ALT, mu=50, alternative="less")

t.test(df$ALTmu=50alternative="less") : T-test를 시행한다. 검정할 변수는 df 데이터의 ALT변수다. ALT의 모평균이 50이라고 할 수 있는지 검정할 것이다. 대립 가설은 "모평균이 50보다 작다"이다.

 

결과

	One Sample t-test

data:  df$ALT
t = -94.014, df = 999, p-value < 2.2e-16
alternative hypothesis: true mean is less than 50
95 percent confidence interval:
     -Inf 35.37659
sample estimates:
mean of x 
 35.11594

위의 결과 중 다음 내용에 주목해야 한다.

 

One Sample t-test

data:  df$ALT t = -94.014, df = 999, p-value < 2.2e-16

alternative hypothesis: true mean is less than 50

95 percent confidence interval:

     -Inf 35.37659

sample estimates:

mean of x   35.11594 

 

P-value는 $2.2 \times 10^{-16}$보다 작다. 즉, 대립 가설 (모평균의 참값은 50보다 작다.)이 합리적이라는 말이다. 평균이 50이라고 보는 것보다는 50보다 작다고 보는 것이 더 합리적이라는 뜻이다. 실제로 평균은 35.1159라고 적혀있으므로 합당한 통계분석임을 알 수 있다.

 

 

3)  $H_1=$ 모집단의 평균은 50보다 작거나 크다. (양측 검정)

t.test(df$ALT, mu=50, alternative="two.sided")

t.test(df$ALT, mu=50)

위 두 개의 코드는 같은 결과를 내준다. 왜냐하면 양측 검정 (alternative="two.sided")가 기본 값이어서 alternative 옵션을 지정해주지 않으면 alternative="two.sided"라고 인식하고 코드가 돌아가기 때문이다.

 

 

t.test(df$ALTmu=50alternative="two.sided") : T-test를 시행한다. 검정할 변수는 df 데이터의 ALT변수다. ALT의 모평균이 50이라고 할 수 있는지 검정할 것이다. 대립 가설은 "모평균은 50이 아니다"이다.

 

결과

One Sample t-test

data:  df$ALT
t = -94.014, df = 999, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 50
95 percent confidence interval:
 34.80527 35.42661
sample estimates:
mean of x 
 35.11594

 

위의 결과 중 다음 내용에 주목해야 한다.

 

One Sample t-test 

data:  df$ALT t = -94.014, df = 999, p-value < 2.2e-16 

alternative hypothesis: true mean is not equal to 50 

95 percent confidence interval:

     34.80527 35.42661 

sample estimates: 

mean of x   35.11594 

 

P-value는 $2.2 \times 10^{-16}$보다 작다. 즉, 대립 가설 (모평균의 참값은 50dㅣ 아니다.)이 합리적이라는 말이다. 평균이 50이라고 보는 것보다는 50보다 작거나 크다고 보는 것이 더 합리적이라는 뜻이다. 실제로 평균은 35.1159라고 적혀있으므로 합당한 통계분석임을 알 수 있다. 실제로 평균은 35.1159라고 적혀있으므로 합당한 통계분석임을 알 수 있다.

 

주로 양측 검정이 이용되나, 사실 일표본 T검정은 거의 시행되지 않고 다음 글로 나올 독립 표본 T검정이 더 많이 사용된다. 

 

코드 정리

#워킹 디렉토리 지정
setwd("C:/Users/user/Documents/Tistory_blog")

#데이터 불러오기
install.packages("readr")
library("readr")
df<-read_csv("Data.csv")



##정규성 검정##
# 1) Q-Q plot 그리기
qqnorm(df$ALT)
qqline(df$ALT)

# 2) 히스토그램 그리기
hist(df$ALT, prob=TRUE)
ALTrange<-seq(min(df$ALT),max(df$ALT),length=max(max(df$ALT)-min(df$ALT),100))
ND<-dnorm(ALTrange,mean=mean(df$ALT),sd=sd(df$ALT))
lines(ALTrange, ND, lwd=2)

# 3) Shapiro-Wilk test 시행하기
shapiro.test(df$ALT)


##T-test 1
# 단측 (대립가설=크다) 검정 t-test
t.test(df$ALT, mu=50, alternative="greater")

##T-test 2
# 단측 (대립가설=작다) 검정 t-test
t.test(df$ALT, mu=50, alternative="less")

##T-test 3
# 양측 검정 t-test
t.test(df$ALT, mu=50, alternative="two.sided")
t.test(df$ALT, mu=50)

 

R 일표본 T검정 (One-sample T-test) 정복 완료!

 

작성일: 2022.11.03.

최종 수정일: 2022.11.29.

이용 프로그램: R 4.1.3

RStudio v1.4.1717

RStudio 2021.09.1+372 "Ghost Orchid" Release 

운영체제: Windows 10, Mac OS 10.15.7

반응형
반응형

[SAS] 일원 배치 분산 분석 (One-way ANOVA, ANalysis Of VAriance) - PROC GLM, PROC ANOVA

 이전 포스팅에서 음주 여부에 따른 ALT 평균에 차이가 있는지를 알아보기 위해서는 T test를 사용한다는 내용을 알아보았다. 2022.10.04 - [모평균 검정/SAS] - [SAS] 독립 표본 T검정 (Independent samples T-test) - PROC TTEST

 음주 여부는 "1) 음주자", "2) 비음주자"로 나뉘는 이분형 변수다. 즉, 음주자의 ALT 평균과 비음주자의 ALT 평균을 비교한 것이었다. 

 

 만약, 두 개 이상의 범주로 나뉘는 변수에 대해 모평균의 검정을 하고자 한다면 어떻게 해야 할까? 예를 들어, 흡연 상태를 "1) 비흡연자", "2) 과거 흡연자", "3) 현재 흡연자"로 나누었고, 각 그룹의 폐기능 검사 중 하나인 FVC (Forced Vital Capacity)의 평균에 차이가 있는지 알아보고자 한다. 

 

 이때 사용할 수 있는 방법이 일원 배치 분산 분석인 One-way ANOVA (Analysis of Variance)이며 이번 포스팅에서 다뤄볼 주된 내용이다.

 

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.10.11)

 

분석용 데이터 (update 22.10.11)

2022년 10월 11일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 범주형 자료 분석 - 모평균 검정 - 반복 측정 자료 분석 - 통계

medistat.tistory.com

 

시작하기 위해 라이브러리를 만들고, 파일을 불러온다.

라이브러리 만드는 방법: 2022.08.05 - [통계 프로그램 사용 방법/SAS] - [SAS] 라이브러리 만들기 - LIBNAME

파일 불러오는 방법: 2022.08.05 - [통계 프로그램 사용 방법/SAS] - [SAS] 데이터 불러오기 및 저장하기 - PROC IMPORT, PROC EXPORT

 

*라이브러리 지정하기;
LIBNAME hong "C:/Users/User/Documents/Tistory_blog";

*파일 불러오기;
PROC IMPORT
DATAFILE="C:\Users\user\Documents\Tistory_blog\Data.xlsx"
DBMS=EXCEL
OUT=hong.df
REPLACE;
RUN;

 

목표: 흡연 상태에 따라 FVC 평균은 서로 다르다고 모집단 수준에서 말할 수 있는가?

 

전제조건 (정규성, 등분산성)

 ANOVA에는 두 가지 전제조건이 필요하다. 본 포스팅의 예시에 맞추어 설명하면 다음과 같다.

 1) 정규성: 현재 흡연자, 과거 흡연자, 비흡연자 별로 FVC의 분포를 보았을 때 각각의 분포는 정규성을 따른다. 

 2) 등분산성: 현재 흡연자, 과거 흡연자, 비흡연자 별로 FVC의 분포를 보았을 때 분산은 서로 같다.

 

따라서 각각을 검정해야 한다.

 

1) 정규성 

정규성 검정에 관한 내용은 다음 링크에서 확인할 수 있다.

2022.08.12 - [기술 통계/SAS] - [SAS] 정규성 검정 - PROC UNIVARIATE

 

 SMOK: 흡연 상태

 FVC_pPRED: FVC의 백분위 예측치

PROC UNIVARIATE DATA=hong.df NORMAL PLOT;
CLASS SMOK;
VAR FVC_pPRED;
HISTOGRAM FVC_pPRED/ NORMAL (MU=EST SIGMA=EST);
RUN;

 

결과 - 1) 비흡연자 (SMOK=0)

.

 표본의 수가 2,000이 되지 않으므로 Shapiro-Wilk 통계량의 p-value를 보아야 하고, 이는 0.05보다 크다. Q-Q plot, 히스토그램에서 정규성을 띠지 않는다고 할만한 근거가 없으므로 정규성을 따른다고 결론을 짓는다.

 

결과 - 2) 과거 흡연자 (SMOK=1)

 표본의 수가 2,000이 되지 않으므로 Shapiro-Wilk 통계량의 p-value를 보아야 하고, 이는 0.05보다 크다. Q-Q plot, 히스토그램에서 정규성을 띠지 않는다고 할만한 근거가 없으므로 정규성을 따른다고 결론을 짓는다.

 

결과 - 3) 현재 흡연자 (SMOK=2)

 표본의 수가 2,000이 되지 않으므로 Shapiro-Wilk 통계량의 p-value를 보아야 하고, 이는 0.05보다 크다. Q-Q plot, 히스토그램에서 정규성을 띠지 않는다고 할만한 근거가 없으므로 정규성을 따른다고 결론을 짓는다.

 

 따라서 정규성 전제조건은 만족한다고 할 수 있다.

 

2) 등분산성 

등분산성 검정에 관한 내용은 다음 링크에서 확인할 수 있다.

2022.10.04 - [모평균 검정/SAS] - [SAS] 독립 표본 T검정 (Independent samples T-test) - PROC TTEST

2022.10.05 - [모평균 검정/SAS] - [SAS] 등분산성 검정 (Homogeneity of variance) - PROC GLM

 

위 두 개 글에서 제시한 네 가지 방법 모두를 사용하여 등분산 검정을 시행해 보겠다.

*1) Levene;
PROC GLM DATA=hong.df;
CLASS SMOK;
MODEL FVC_pPRED=SMOK;
MEANS SMOK/ HOVTEST=LEVENE(TYPE=ABS);
RUN;

*2) O'Brien;
PROC GLM DATA=hong.df;
CLASS SMOK;
MODEL FVC_pPRED=SMOK;
MEANS SMOK/ HOVTEST=OBRIEN;
RUN;

*3) Brwon and Forsythe;
PROC GLM DATA=hong.df;
CLASS SMOK;
MODEL FVC_pPRED=SMOK;
MEANS SMOK/ HOVTEST=BF;
RUN;

*4) Bartlett;
PROC GLM DATA=hong.df;
CLASS SMOK;
MODEL FVC_pPRED=SMOK;
MEANS SMOK/ HOVTEST=BARTLETT;
RUN;
QUIT;

 

결과

1) Levene

2) O'Brien

3) Brown and Forsythe

4) Bartlett

 어떤 방법을 사용하더라도 p-value>0.05로 등분산성이라는 전제조건을 만족한다는 것을 알 수 있다. 따라서 일원 배치 분산 분석 (One-way ANOVA)를 시행할 수 있다.

 

 

일원 배치 분산 분석 (One way ANOVA) 코드

*1)사후분석- 본페로니;
PROC GLM DATA=hong.df;
CLASS SMOK;
MODEL FVC_pPRED=SMOK;
MEANS SMOK/BON CLDIFF;
RUN;
QUIT;

*2)사후분석- 본페로니;
PROC GLM DATA=hong.df;
CLASS SMOK;
MODEL FVC_pPRED=SMOK;
LSMEANS SMOK/ADJUST=BON PDIFF;
RUN;
QUIT;

*3)사후분석- 던칸;
PROC GLM DATA=hong.df;
CLASS SMOK;
MODEL FVC_pPRED=SMOK;
MEANS SMOK/DUNCAN;
RUN;
QUIT;

세 가지 방법으로 코드를 소개했지만, 모두 대동소이하다. 첫 세 줄이 ANOVA 코드로 모두 동일하다. 네 번째 줄은 사후 분석 및 출력 방법에 대한 코드다. 사후 분석에 대해서는 아래에 설명하도록 하겠다.

 

1) 사후 분석: 본페로니(BONferroni) 방법을 사용하고, 차이의 신뢰구간 (Confidence Limits for DIFFerence)을 나타내라.
PROC GLM DATA=hong.df; : 일반화 선형 모델 (GLM)을 시행하겠다. 데이터는 hong라이브러리의 df파일을 사용하라.
CLASS SMOK; : 흡연 상태에 따라 분산 분석 (ANOVA)를 시행하라.
MODEL FVC_pPRED=SMOK; : 흡연 상태(SMOK)에 따라 FVC % 예측치(FVC_pPRED)의 평균에 차이가 있는지 나타내라.
MEANS SMOK/BON CLDIFF; : 흡연 상태(SMOK)에 따른 FVC % 예측치(FVC_pPRED)의 평균을 보여달라. 사후 분석은 본페로니(BONferroni, BON) 방법을 사용하고, 차이의 신뢰구간 (Confidence Limits for DIFFerence, CLDIFF)을 보여달라

 - CLDIFF를 사용하기 위해서는 "MEANS"로 시작하는 구문을 사용해야 한다.

2) 사후 분석: 본페로니(BONferroni) 방법을 사용하고, 차이의 p-value (P-value for DIFFerence)를 나타내라.
PROC GLM DATA=hong.df; : 일반화 선형 모델 (GLM)을 시행하겠다. 데이터는 hong라이브러리의 df파일을 사용하라.
CLASS SMOK; : 흡연 상태에 따라 분산 분석 (ANOVA)를 시행하라.
MODEL FVC_pPRED=SMOK; : 흡연 상태(SMOK)에 따라 FVC % 예측치(FVC_pPRED)의 평균에 차이가 있는지 나타내라.
LSMEANS SMOK/ADJUST=BON PDIFF; : 흡연 상태(SMOK)에 따른 FVC % 예측치(FVC_pPRED)의 최소 제곱 평균(Least Squares MEANS, LSMEANS)을 보여달라. 사후 분석은 본페로니(BONferroni, BON) 방법을 사용하고, 차이의 p-value (P-value for DIFFerence)를 보여달라

 - PDIFF를 사용하기 위해서는 "LSMEANS"로 시작하는 구문을 사용해야 한다.

3) 사후 분석: 던칸(DUNCAN) 방법을 사용하라.
PROC GLM DATA=hong.df; : 일반화 선형 모델 (GLM)을 시행하겠다. 데이터는 hong라이브러리의 df파일을 사용하라.
CLASS SMOK; : 흡연 상태에 따라 분산 분석 (ANOVA)를 시행하라.
MODEL FVC_pPRED=SMOK; : 흡연 상태(SMOK)에 따라 FVC % 예측치(FVC_pPRED)의 평균에 차이가 있는지 나타내라.
MEANS SMOK/DUNCAN; : 흡연 상태(SMOK)에 따른 FVC % 예측치(FVC_pPRED)의 평균을 보여달라. 사후 분석은 던칸(DUNCAN) 방법을 사용하라.

결과

 위 세 개의 코드에서 공통적인 부분 (ANOVA)의 결과는 다음과 같다.

 

P-value가 <0.0001이므로 귀무가설을 기각하고 대립 가설을 채택한다. 그렇다면 여기에서 귀무가설 및 대립 가설은 무엇이었는가?

 

 귀무가설:  $H0=$ 세 집단의 모평균은 "모두" 동일하다.

 대립 가설:  $H1=$ 세 집단의 모평균이 모두 동일한 것은 아니다.

 

우리는 대립 가설을 채택해야 하므로 "세 집단의 모평균이 모두 동일한 것은 아니다."라고 결론 내릴 것이다.

이 말은, 세 집단 중 두 개씩 골라 비교했을 때, 적어도 한 쌍에서는 차이가 난다는 것이다. 따라서 세 집단 중 두 개씩 골라 비교를 해보아야 하며, 이를 사후 분석 (post hoc analysis)이라고 한다. 세 집단에서 두 개씩 고르므로 가능한 경우의 수는 $_3 C_2=3$이다.

 (1) 비흡연자 vs 과거 흡연자

 (2) 비흡연자 vs 현재 흡연자

 (3) 과거 흡연자 vs 현재 흡연자

 

사후 분석 결과 부분을 보겠다.

 

 

1) 사후 분석으로 본페로니(BONferroni) 방법을 사용하고, 차이의 신뢰구간 (Confidence Limits for DIFFerence)을 나타내라.

 빨간 박스를 먼저 설명하면, 비흡연자 (SMOK=0)가 과거 흡연자 (SMOK=1)보다 FVC가 9.5328%p 더 크다. 그 차이의 95% 신뢰구간은 7.2168에서 11.8499이다. 즉 신뢰구간에 0이 포함되어있지 않으므로 유의미한 차이가 있다고 할 수 있다. 그런 경우 노란 박스의 "***" 표시를 해준다. 밑의 결과도 해석은 같다. 

 즉, 모든 사후 분석 결과에서 유의미한 차이가 있음을 알 수 있다.

 

 

2) 사후 분석: 본페로니(BONferroni) 방법을 사용하고, 차이의 p-value (P-value for DIFFerence)를 나타내라.

 위의 방법 1)에서는 95% 신뢰구간만을 보여주고, p-value가 0.05보다 작은지만 보여줄 뿐, 정확한 p-value의 값은 알려주지 않는다. 정확한 p-value값을 알고 싶은 경우 LSMEANS구문에 PDIFF옵션을 사용하면 된다.

 예를 들어 위의 빨간 박스는 그룹 1 (비흡연자)와 그룹 2(과거 흡연자) 사이의 평균의 차이에 대한 p-value를 보여준다. 위와 다르게 p-value를 보여준다는 장점은 있지만, 그 차이가 얼마인지, 신뢰구간은 얼마인지 보여주지 못하는 단점이 있다. 둘 중 어떤 방법으로 출력을 할지는 개인이 선택을 하면 된다.

 

3) 사후 분석: 던칸(DUNCAN) 방법을 사용하라.

 

이 시점에서 이런 의문이 들 수 있다.

"각각의 그룹별로 평균을 비교하면 되지, 굳이 왜 ANOVA라는 방법까지 사용하는 것인가?"

 아주 논리적인 의문점이다. 하지만, 반드시 ANOVA를 사용해야 한다. 그 이유는 다음과 같다. 본 사례는 흡연 상태에 따른 조합 가능한 경우의 수가 3인데, 각각 유의성의 기준을 $\alpha=0.05$로 잡아보자. 이때 세 번의 비교에서 모두 귀무가설이 학문적인 진실인데(평균에 차이가 없다.), 세 번 모두 귀무가설을 선택할 확률은 $(1-0.95)^3 \approx 0.86$이다. (이해가 어려우면 p-value에 대한 설명 글을 읽고 오길 바란다. 2022.09.05 - [통계 이론] - [이론] p-value에 관한 고찰)

 그런데, ANOVA의 귀무가설은 "모든 집단의 평균이 같다."이다. 따라서 모든 집단의 평균이 같은 것이 학문적 진실일 때, 적어도 한 번이라도 대립 가설을 선택하게 될 확률은 $1-0.86=0.14$가 된다. 즉, 유의성의 기준이 올라가게 되어, 덜 보수적인 결정을 내리게 되고, 다시 말하면 대립 가설을 잘 선택하는 쪽으로 편향되게 된다. 학문적으로는 '다중 검정 (multiple testing)을 시행하면 1종 오류가 발생할 확률이 증가하게 된다.'라고 표현한다.

 따라서, 각각을 비교해보는 것이 아니라 한꺼번에 비교하는 ANOVA를 시행해야 함이 마땅하다. 

 

 여기에서 한 번 더 의문이 들 수 있다.

 "사후 분석을 할 때에는 1종 오류가 발생하지 않는가?"

그렇다. 1종 오류가 발생할 확률이 있으므로, p-value의 기준을 더 엄격하게 (0.05보다 더 작게) 잡아야 한다. P-value를 보정하는 방법은 여러 가지가 나와있는데, 가장 간단하고 제일 보수적인 방법이 위에 예시에서 사용한 본페로니 방법이다. 그 외에도 많은 방법이 있고 다른 방법을 사용하기 위해서는 BON의 자리에 다음 중 하나를 입력하면 된다.

 

코드 내용 코드 내용
BON Bonferroni t test SCHEFFE Scheffe's multiple comparison
DUNCAN Duncan's multiple range test SIDAK Pairwise t test (Sidak's inequality)
DUNNETT Dunnett's two-tailed t test SMM (혹은 GT2) Pairwise comparisons (studentized maximum modulus, Sidak's uncorrelated t-inequality)
DUNNETTL Dunnett's one-tailed t test (less) SNK Student-Newman-Keuls
DUNNELLU Dunnett's one-tailed t test (greater) T (혹은 LSD) Pairwise t test
GABRIEL Gabriel's multiple comparison TUKEY Tukey's studentized range test
REGWQ Ryan-Einot-Habriel-Welsch WALLER Waller-Duncan k-ratio  t test

 

굉장히 많은 방법이 있지만, Bonferroni, Duncan, Scheffe, Tukey가 가장 많이 쓰이는 듯하다. 이 중에 Duncan, Tukey는 그룹별로 표본의 수가 같을 때만 사용해야 하므로 유의하도록 한다.

 

그중에 DUNCAN옵션을 사용하여 산출한 결과는 다음과 같다.

막대기가 연결되어 있으면 평균이 서로 다르다고 할 수 없는 것이고, 끊어져 있으면 서로 다르다고 할 수 있는 것이다. 위의 경우 비흡연자 (SMOK=0)은 파란 막대기가 있고, 과거 흡연자 (SMOK=1)는 빨간 막대기로 파란 막대기와 분리되어 있다. 따라서 과거 흡연자의 FVC는 비흡연자의 FVC와 유의하게 다르다고 할 수 있다. 현재 흡연자 (SMOK=2)는 초록 막대기로 파랑 및 빨간 막대기와 분리되어 있다. 따라서 현재 흡연자의 FVC는 비흡연자 및 과거 흡연자의 FVC와 유의하게 다르다고 할 수 있다. 즉 모든 그룹에서 FVC의 평균은 다르다고 할 수 있다.

 

어떤 사후 분석을 쓸 것인가

 이 논의에 대해 정답이 따로 있는 것은 아니다. 적절한 방법을 사용하여 논문에 제시하면 되고, 어떤 것이 정답이라고 콕 집어 이야기할 수는 없다. 다만, 사후 분석 방법이 여러 가지가 있다는 것은 '사후 분석 방법에 따라 산출되는 결과가 달라질 수 있다.'는 것을 의미하고, 심지어는 '어떤 사후 분석 방법을 채택하냐에 따라 유의성 여부가 달라질 수도 있다.'는 것을 의미한다. 심지어, ANOVA에서는 유의한 결과가 나왔는데, 사후 분석을 해보니 유의한 차이를 보이는 경우가 없을 수도 있다. 따라서 어떤 사후 분석 방법 결과에 따른 결과인지 유의하여 해석할 필요가 있다.

 

PROC ANOVA

 PROC ANOVA에 대한 설명을 듣고자 본 포스팅에 들어온 독자도 있을 것이다. 하지만 이는 다루지 않을 것이다. 왜냐하면 PROC ANOVA는 그룹별로 표본의 수가 같아야만 사용할 수 있는, 아주 제한적인 상황에서만 쓸 수 있는 코드이기 때문이다. 따라서 웬만하면 PROC GLM을 사용하길 권한다. 만약 사용하고자 한다면, PROC GLM을 PROC ANOVA로만 바꾸어 다음과 같이 사용하면 된다.

PROC ANOVA DATA=hong.df;
CLASS SMOK;
MODEL FVC_pPRED=SMOK;
MEANS SMOK;
RUN;

SAS 일원 배치 분산 분석 (One-way ANOVA, ANalysis Of VAriance) 정복 완료!

작성일: 2022.10.12.

최종 수정일: 2022.11.30.

이용 프로그램: SAS v9.4

운영체제: Windows 10

반응형
반응형

 

[SAS] 대응 표본 T검정 (Paired samples T-test) - PROC TTEST

 

  세상에 존재하는 모든 사람을 대상으로 연구를 하면 참 좋겠지만, 시간적인 이유로, 그리고 경제적인 이유로 일부를 뽑아서 연구를 진행할 수밖에 없다. 모든 사람을 모집단이라고 하고, 뽑힌 일부를 표본이라고 한다. 우리는 표본으로 시행한 연구로 모집단에 대한 결론을 도출해내고자 할 것이다. 

 1000명에게 피검사를 시행하였고, 간 기능 검사의 일환으로 ALT 수치를 모았다. 그리고, 간 기능 개선제를 복용시킨 뒤 ALT 검사를 다시 시행하였다. 이 데이터를 기반으로 1000명이 기원한 모집단 인구에서 간 기능 개선제 복용은 ALT 수치를 개선시키는지 여부를 확인하는 것이 T-test이다. T-test는 크게 세 가지로 나눌 수 있다.

 

 1) 일표본 T검정 (One sample T-test) : 2022.09.30 - [모평균 검정/SAS] - [SAS] 일표본 T검정 (One-sample T-test) - PROC TTEST

 : 모집단의 평균이 특정 값이라고 할 수 있는가?

 예) 모집단의 ALT 평균이 50이라고 할 수 있는가?

 

 

 2) 독립 표본 T검정 (Independent samples T-test, Two samples T-test) : 2022.10.04 - [모평균 검정/SAS] - [SAS] 독립 표본 T검정 (Independent samples T-test) - PROC TTEST

 : 두 모집단의 평균이 다르다고 할 수 있는가?

 예) 고혈압 환자와 일반인의 수축기 혈압 평균이 서로 다르다고 할 수 있는가?

 

 3) 대응 표본 T검정 (Paired samples T-test) 

 : 모집단의 짝지어진 변수들의 평균이 다르다고 할 수 있는가?

예) 간 기능 개선제 복용 전 ALT 평균은 간기는 개선제 복용 후 ALT 평균과 다르다고 할 수 있는가?

 

이번 포스팅에서는 대응 표본 T검정 (Paired samples T-test)에 대해 알아볼 것이다.

 

 

왜 대응 표본 T 검정이 필요할까?

 독립 표본 T 검정을 공부한 사람이라면 대응 표본 T 검정의 필요성에 대해 의구심을 가질 수 있다. 간 기능 개선제 복용 전 ALT와 간 기능 개선제 복용 후 ALT에 대해 독립 표본 T 검정을 돌리면 될 것이라고 생각할 수 있기 때문이다. 하지만 결론적으로 그렇게 하면 안 된다. 예를 들어 다음 사례를 봐보자.

환자 번호 시험 공부 전 시험 점수 시험 공부 후 시험 점수
1번 0.1 0.2
2번 0.2 0.3
... $x$ $x+0.1$
999번 99.9 100.0
1000번 100.0 100.0

1000명이 시험을 보았고, 원래 100점을 받았던 1000번을 제외한 모든 사람이 시험공부로 0.1점이 올랐다. 1000명 중 999명의 점수가 올랐으므로 공부는 시험 점수를 올리는데 유의미한 영향력을 미칠 수 있다고 결론이 나야 합당할 것이다. 하지만 공부 전 후 시험 점수의 평균은 둘 다 50점에 가깝고 큰 변화는 없는 것처럼 보인다. 따라서 독립 표본 T 검정으로 내린 결론은 잘못된 것이다. 즉, 대응되는 (paired) 경우 그 값 자체보다는 둘 사이의 차이에 집중해야 한다.

 

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.10.06)

 

분석용 데이터 (update 22.10.06)

2022년 08월 29일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 통계 프로그램 사용 방법 1) 엑셀 파일 2) CSV 파일 3) 코드북

medistat.tistory.com

 

시작하기 위해 라이브러리를 만들고, 파일을 불러온다.

라이브러리 만드는 방법: 2022.08.05 - [통계 프로그램 사용 방법/SAS] - [SAS] 라이브러리 만들기 - LIBNAME

파일 불러오는 방법: 2022.08.05 - [통계 프로그램 사용 방법/SAS] - [SAS] 데이터 불러오기 및 저장하기 - PROC IMPORT, PROC EXPORT

 

*라이브러리 지정하기;
LIBNAME hong "C:/Users/User/Documents/Tistory_blog";

*파일 불러오기;
PROC IMPORT
DATAFILE="C:\Users\user\Documents\Tistory_blog\Data.xlsx"
DBMS=EXCEL
OUT=hong.df
REPLACE;
RUN;

 

전제: 정규성 검정 (차이)

 대응 표본 T 검정 (Paired samples T test)의 전제 조건은 검정하고자 하는 두 변수의 차이가 정규분포를 따른다는 것이다. 따라서 ALT와 간 기능 개선제 복용 후 ALT(ALT_POSTMED)의 차이를 구하고 정규성 검정을 시행한다. 

차이를 구하는 방법: 2022.10.06 - [통계 프로그램 사용 방법/SAS] - [SAS] 변수 계산 (산술 연산) - DATA, SET

정규성 검정을 하는 방법: 2022.08.12 - [기술 통계/SAS] - [SAS] 정규성 검정 - PROC UNIVARIATE

 

DATA hong.df;
SET hong.df;
*차이를 나타내는 변수는 ALT_DIF;
ALT_DIF=ALT-ALT_POSTMED;
RUN;

PROC UNIVARIATE DATA=hong.df NORMAL PLOT;
VAR ALT_DIF;
HISTOGRAM ALT_DIF/ NORMAL (MU=EST SIGMA=EST);
RUN;

 

결과

N 수는 1000으로 2000 미만이므로 Shapiro-Wilk 통계량의 p-value는 0.05 이상이다. Q-Q plot상에서 대부분의 데이터가 직선상에 있으며, 히스토그램 또한 정규분포를 따르는 것처럼 보이므로 두 변수의 차는 정규분포를 따른다고 할 수 있다. 따라서 대응 표본 T 검정 (Paired samples T test)를 시행할 수 있다.

 

대응 표본  T 검정 코드

PROC TTEST DATA=hong.df;
PAIRED ALT*ALT_POSTMED;
RUN;

 

PROC TTEST DATA=hong.df; : T-test를 시행할 것이고, 데이터는 hong 라이브러리의 df 데이터를 이용하라.
PAIRED ALT*ALT_POSTMED; : 대응 표본 T 검정(PAIRED)을 시행할 것이고, 각각의 변수는 ALT와 ALT_POSTMED이다.

 

결과

ALT의 평균이 ALT_POSTMED의 평균에 비해 2.8363만큼 더 크다. 그리고 그 차이의 p-value는 <0.0001로 유의미하다. 따라서 귀무가설을 기각할 수 있다. 본 검정의 귀무가설은 다음과 같다.

 귀무가설:  $H_0$=간 기능

이를 기각하므로 동일하지 않다고 할 수 있고, 간 기능 개선제 복용 후 ALT가 2.8363만큼 낮아져 간 기능에 호전이 있다고 평가할 수 있다.  

 

$$ALT = ALT\_POSTMED$$

$$ \therefore ALT\_DIF=ALT-ALT\_POSTMED=0$$

 

즉, 대응 표본 T 검정을 시행하는 것은 ALT_DIF변수의 평균이 0과 같은지 검정하는 일표본 T검정 (One sample T test)와 같은 것이다. 일표본 T 검정은 다음 링크에서 확인할 수 있다. 2022.09.30 - [모평균 검정/SAS] - [SAS] 일표본 T검정 (One-sample T-test) - PROC TTEST

 

일표본 T 검정으로 대응 표본 T 검정 시행하기 - 코드

PROC TTEST DATA=hong.df H0=0;
VAR ALT_DIF;
RUN;

 

결과

위와 같은 결과를 내는 것을 알 수 있다.

 

코드 정리

*라이브러리 지정하기;
LIBNAME hong "C:/Users/User/Documents/Tistory_blog";

*파일 불러오기;
PROC IMPORT
DATAFILE="C:\Users\user\Documents\Tistory_blog\Data.xlsx"
DBMS=EXCEL
OUT=hong.df
REPLACE;
RUN;

*차이를 나타내는 변수 만들기;
DATA hong.df;
SET hong.df;
*차이를 나타내는 변수는 ALT_DIF;
ALT_DIF=ALT-ALT_POSTMED;
RUN;

*정규성 검정 시행하기;
PROC UNIVARIATE DATA=hong.df NORMAL PLOT;
VAR ALT_DIF;
HISTOGRAM ALT_DIF/ NORMAL (MU=EST SIGMA=EST);
RUN;

*대응 표본 T 검정 시행하기 (1);
PROC TTEST DATA=hong.df;
PAIRED ALT*ALT_POSTMED;
RUN;

*대응 표본 T 검정 시행하기 (2);
PROC TTEST DATA=hong.df H0=0;
VAR ALT_DIF;
RUN;

 

SAS 대응 표본 T검정 (Paired samples T-test) 정복 완료!

작성일: 2022.10.07.

최종 수정일: 2022.10.07.

이용 프로그램: SAS v9.4

운영체제: Windows 10

 

반응형
반응형

[SAS] 변수 계산 (비교 연산) - DATA, SET, GE, LE, NE, IN, MIN, MAX

조건에 따라 변수를 바꾸고 싶을 때가 있다. 예를 들어, 인구를 1) AST 40 이상, 2) AST 20 이상 40 미만, 3) AST 20 미만으로 나누고 싶을 수 있다. 이런 경우에 사용하는 비교 연산에 대해 알아보고자 한다. 내용은 다음과 같다.

 

1) >, <, >=, <=, =, ^=

2) MAX

3) MIN

4) IN, NOT IN

 

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.10.06)

 

분석용 데이터 (update 22.10.06)

2022년 08월 29일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 통계 프로그램 사용 방법 1) 엑셀 파일 2) CSV 파일 3) 코드북

medistat.tistory.com

시작하기 위해 라이브러리를 만들고, 파일을 불러온다.

라이브러리 만드는 방법: 2022.08.05 - [통계 프로그램 사용 방법/SAS] - [SAS] 라이브러리 만들기 - LIBNAME

파일 불러오는 방법: 2022.08.05 - [통계 프로그램 사용 방법/SAS] - [SAS] 데이터 불러오기 및 저장하기 - PROC IMPORT, PROC EXPORT

 

*라이브러리 지정하기;
LIBNAME hong "C:/Users/User/Documents/Tistory_blog";

*파일 불러오기;
PROC IMPORT
DATAFILE="C:\Users\user\Documents\Tistory_blog\Data.xlsx"
DBMS=EXCEL
OUT=hong.df
REPLACE;
RUN;

 

상황: 전체 인구를 1) AST 40 이상, 2) AST 20 이상 40 미만, 3) AST 20 미만으로 나누고 싶을 때

 AST는 연속 변수로 코딩되어 있음을 코드북에서 확인할 수 있다. 우리는 ALT_S3이라는 변수를 새로 만들어 다음과 같이 코딩되기를 바란다.

AST_S2 0 1 2
AST 40 이상 20 이상 40 미만 20 미만

 이를 실행하는 방법은 여러 가지가 있는데, 코드를 하나씩 살펴보면 다음과 같다. 변수를 바꿀 때 쓰는 DATA-SET구문에 대한 설명은 다음 링크를 확인하길 바란다. 2022.10.06 - [통계 프로그램 사용 방법/SAS] - [SAS] 변수 계산 (산술 연산) - DATA, SET

 

가장 기초적인 방법 (하지만 잘못된 코드다)

DATA hong.df_new;
SET hong.df;

*AST를 기준으로 셋으로 나누기;
IF AST>=40 THEN AST_S3=0;
	ELSE IF AST<40 AND AST>=20 THEN AST_S3=1;
	ELSE IF AST<20 THEN AST_S3=2;

RUN;

조건문 (IF, ELSE IF)와 AND 연산자에 대한 내용은 다음 링크를 확인하길 바란다. 2022.10.06 - [통계 프로그램 사용 방법/SAS] - [SAS] 변수 계산 (논리 연산) - DATA, SET, IF, ELSE IF, ELSE, AND, OR

 

IF AST>=40 THEN AST_S3=0; : AST가 40 이상이면 AST_S3의 값은 0을 부여한다.
     ELSE IF AST<40 AND AST>=20 THEN AST_S3=1; : AST가 40 미만이면서 20 이상이면 AST_S3의 값은 1을 부여한다.
     ELSE IF AST<20 THEN AST_S3=2; : AST가 20 미만이면 AST_S3의 값은 2를 부여한다.

 

 

문제점: 결측치

얼핏 보면 이상할 것이 없어 보이는 이 코드가 왜 잘못되었다고 하는 것일까?

문제는 결측치다.

 

다음 코드로 AST의 분포를 확인해보면 8개의 결측치가 있음을 알 수 있다. PROC MEANS에 관한 내용은 다음 링크에서 확인할 수 있다. 2022.09.23 - [기술 통계/SAS] - [SAS] 기술 통계 (평균, 표준편차, 표준오차, 최댓값, 최솟값, 중위수, 분위수 등) - PROC UNIVARIATE, PROC MEANS

 

결측치 확인 코드

PROC MEANS DATA=hong.df NMISS; 
VAR AST; 
RUN;

 

결과

 

결측치 확인 코드

그런데, AST_S3에는 결측치가 없음을 다음 코드에서 확인할 수 있다.

PROC MEANS DATA=hong.df_new NMISS; 
VAR AST_S3; 
RUN;

 

결과

 

결측치 향방 확인하는 코드

그렇다면 결측치는 모두 어디에 가있는 것일까? SAS에서는 결측치를 "."으로 표현한다. 결측치를 1로, 결측치가 아닌 경우를 0으로 코딩한 변수 AST_MISS를 만들고, AST_S3과의 분할표를 만들어 확인해보자. 분할표에 대한 내용은 다음 링크를 확인하길 바란다. 2022.08.18 - [기술 통계/SAS] - [SAS] 도수분포표 (Frequency table), 분할표 (Contingency table) 만들기 - PROC FREQ

DATA hong.df_new;
SET hong.df;

*AST를 기준으로 셋으로 나누기;
IF AST>=40 THEN AST_S3=0;
	ELSE IF AST<40 AND AST>=20 THEN AST_S3=1;
	ELSE IF AST<20 THEN AST_S3=2;

*AST 결측 여부 확인 변수 만들기;
IF AST=. THEN AST_MISS=1;
	ELSE AST_MISS=0;

RUN;

*분할표 만들기;
PROC FREQ DATA=hong.df_new; 
TABLE AST_S3 * AST_MISS / NOROW NOCOL NOPERCENT; 
RUN;

 

결과

결측값 8개에 AST_S3의 값으로 2가 부여되어 있다. 이는 명백히 "ELSE IF AST<20 THEN AST_S3=2;" 코드에서 기인한 문제다. 즉, SAS는 AST의 결측치에 대해 "AST<20"가 옳다고 (TRUE) 본 것이다. 이를 해결하기 위해서는 코드를 다음과 같이 작성해야 한다. (두 가지 방법 모두 옳다.)

 

코드

*방법 1: "AST>=0"을 추가하기;
DATA hong.df_new;
SET hong.df;

*AST를 기준으로 셋으로 나누기;
IF AST>=40 THEN AST_S3=0;
	ELSE IF AST<40 AND AST>=20 THEN AST_S3=1;
	ELSE IF AST<20 AND AST>=0 THEN AST_S3=2;

RUN;

*방법 2: "AST^=."을 추가하기;
DATA hong.df_new;
SET hong.df;

*AST를 기준으로 셋으로 나누기;
IF AST>=40 THEN AST_S3=0;
	ELSE IF AST<40 AND AST>=20 THEN AST_S3=1;
	ELSE IF AST<20 AND AST^=. THEN AST_S3=2;

RUN;

방법 1은 "AST>=0"을 추가하여 결측치는 해당 조건을 만족하지 않게 조정한 것이다.

방법 2는 "AST^=."을 추가하였다.  "^="는 같지 않다는 뜻이므로, AST가 결측이 아니라는 뜻을 의미한다. 

 

각 방법에는 세 가지 조건이 존재한다. 방법 2를 기준으로 한다면 

조건 1: AST>=40

조건 2: AST<40 AND AST>=20

조건 3: AST<20 AND AST^=.

인데, 결측치는 세 조건 모두 부합하지 않았다. 따라서 결측치는 THEN 뒤로 넘어가 본 적이 없다. THEN 뒤에서는 AST_S3의 값을 배정받으므로 결측치는 AST_S3값을 배정받아본 적이 없는 것이다. 따라서 AST가 결측인 경우에는 AST_S3 또한 결측치가 된다. 

 

>, <, >=, <=등은 실제 코딩을 하다 보면 매우 입력하기 귀찮아지는데, 다음과 같은 영문 부호로 이를 대체할 수도 있다.

연산 부호 영문 부호
= EQ 같다
~=
^=
NE 같지 않다
> GT 크다
< LT 작다
>= GE 크거나 같다
<= LE 작거나 같다
^>
~>
NG 크지 않다
^<
~<
NL 작지 않다

 

 

 

MAX, MIN

최댓값을 구하는 MAX, 최솟값을 구하는 MIN 함수다. 간 기능 지표 (ALT, AST) 중 큰 것을 고른 LIVER_MAX변수와, 작은 것을 고른 LIVER_MIN변수를 새로 만든다고 하자. 그럼 코드는 다음과 같다.

DATA hong.df_new;
SET hong.df;

*최댓값;
LIVER_MAX=MAX(AST, ALT);
*최솟값;
LIVER_MIN=MIN(AST, ALT);
RUN;

LIVER_MAX=MAX(AST, ALT); : AST와 ALT 중 더 큰 값을 LIVER_MAX에 저장하라.
LIVER_MIN=MIN(AST, ALT); : AST와 ALT 중 더 작은 값을 LIVER_MIN에 저장하라.

 

여기에서는 결측치를 어떻게 처리할까?

AST, ALT 둘 중 하나가 결측치인 경우 MIN, MAX모두 결측치가 아닌 값을 반환한다. 

만약, AST, ALT 둘 다 결측치인 경우 MIN, MAX은 결측치를 반환한다.

 

 

상황: 전체 인구를 1) 과거 흡연자, 2) 비흡연자 및 현재 흡연자로 나누고 싶을 때

 이전 포스팅에서 이런 경우 조건문(IF)을 이용하여 나눌 수 있다고 하였다. 2022.10.06 - [통계 프로그램 사용 방법/SAS] - [SAS] 변수 계산 (논리 연산) - DATA, SET, IF, ELSE IF, ELSE, AND, OR

위 글에서도 간단하게 설명했지만 IN이라는 연산자를 이용하면 매우 빠르게 할 수 있다.

 

IN은 다음과 같이 사용할 수 있다. 

 

1) IN (0 2) : 0 혹은 2에 해당하는 경우

2) IN (0:2) : 0 이상 2 이하 정수에 해당하는 경우

 

IN을 NOT IN으로 바꾸면 "해당하지 않는 경우"로 바뀐다.

*결측치가 없음을 확인하기;
PROC MEANS DATA=hong.df NMISS;
VAR SMOK;
RUN;

DATA hong.df_new;
SET hong.df;

*과거 흡연자 // 현재 흡연자 + 비흡연자로 나누기;
IF SMOK IN (0 2) THEN SMOK_S2=0;
	ELSE SMOK_S2=1;

RUN;

IF SMOK IN (0 2) THEN SMOK_S2=0; : SMOK가 0 혹은 2인 경우 SMOK_S2에는 0을 부여한다.
ELSE SMOK_S2=1; : 그 외에는 SMOK_S2에 1을 부여한다.

 

SAS 변수 계산 (비교 연산) 정복 완료!

작성일: 2022.10.06.

최종 수정일: 2022.10.06.

이용 프로그램: SAS v9.4

운영체제: Windows 10

반응형
반응형

[SAS] 변수 계산 (논리 연산) - DATA, SET, IF, ELSE IF, ELSE, AND, OR

 조건에 따라 변수를 바꾸고 싶을 때가 있다. 예를 들어, 인구를 1) 음주 & 흡연자, 2) 음주 & 과거 흡연자, 3) 그 외로 나누고 싶을 수 있다. 이런 경우에 사용하는 논리 연산에 대해 알아보고자 한다. 내용은 다음과 같다.

 

1) 조건문 (IF, ELSE IF)

2) AND

3) OR

 

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.10.06)

 

분석용 데이터 (update 22.10.06)

2022년 08월 29일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 통계 프로그램 사용 방법 1) 엑셀 파일 2) CSV 파일 3) 코드북

medistat.tistory.com

 

시작하기 위해 라이브러리를 만들고, 파일을 불러온다.

라이브러리 만드는 방법: 2022.08.05 - [통계 프로그램 사용 방법/SAS] - [SAS] 라이브러리 만들기 - LIBNAME

파일 불러오는 방법: 2022.08.05 - [통계 프로그램 사용 방법/SAS] - [SAS] 데이터 불러오기 및 저장하기 - PROC IMPORT, PROC EXPORT

 

*라이브러리 지정하기;
LIBNAME hong "C:/Users/User/Documents/Tistory_blog";

*파일 불러오기;
PROC IMPORT
DATAFILE="C:\Users\user\Documents\Tistory_blog\Data.xlsx"
DBMS=EXCEL
OUT=hong.df
REPLACE;
RUN;

 

몇 가지 예시를 통해 논리 연산에 대해 배워보고자 한다.

 

상황 : 전체 인구를 현재 흡연자와 그 외(과거 흡연자와 비흡연자)로 나누려고 할 때

 데이터의 흡연 변수는 다음과 같이 코딩되어 있음을 코드북에서 확인할 수 있다.

SMOK 0 1 2
내용 비흡연자 과거 흡연자 현재 흡연자

 

우리는 SMOK_S2라는 변수를 새로 만들어 다음과 같이 코딩되기를 바란다.

SMOK_S2 0 1
내용 과거 흡연자와 비흡연자 현재 흡연자

 

즉, SMOK를 SMOK_S2로 바꿀 때, 다음과 같이 바뀌어야 한다.

  SMOK SMOK_S2
비흡연자 0 0
과거 흡연자 1 0
현재 흡연자 2 1

 

SMOK의 0과 1은 SMOK_S2에서 0이어야 하고, SMOK의 2는 SMOK_S2에서 1이어야 한다. 이를 실행하는 방법은 여러 가지가 있는데, 코드를 하나씩 살펴보면 다음과 같다. 변수를 바꿀 때 쓰는 DATA-SET구문에 대한 설명은 다음 링크를 확인하길 바란다. 2022.10.06 - [통계 프로그램 사용 방법/SAS] - [SAS] 변수 계산 (산술 연산) - DATA, SET

 

1) 방법 1 (가장 기초적)

DATA hong.df_new;
SET hong.df;

*현재 흡연자 // 과거 흡연자 + 비흡연자로 나누기 (1);
IF SMOK=0 THEN SMOK_S2=0;
	ELSE IF SMOK=1 THEN SMOK_S2=0;
	ELSE IF SMOK=2 THEN SMOK_S2=1;

RUN;

IF구문은 다음과 같은 구조를 띤다.

IF {조건1} THEN {실행1}                  :"조건 1"을 만족하면 "실행1"을 실시한다.

    ELSE IF {조건2} THEN {실행2}    :"조건 1"을 만족하지 않는 데이터 중 "조건 2"를 만족하면 "실행 2"를 실시한다.

    ...

    ELSE {마지막 실행}                     :"조건1", "조건 2", ... 를 만족하지 않는 데이터는 "마지막 실행"을 실시한다.

 

이를 적용하면 다음과 같이 해석할 수 있다.

 

IF SMOK=0 THEN SMOK_S2=0; : SMOK가 0이면 SMOK_S2는 0이다.
     ELSE IF SMOK=1 THEN SMOK_S2=0; : SMOK가 0이 아닌 사람 중 SMOK가 1인 사람의 SMOK_S2는 0이다.
     ELSE IF SMOK=2 THEN SMOK_S2=1; : SMOK가 0, 1이 아닌 사람 중 SMOK가 1인 사람의 SMOK_S2는 0이다.

 

여기에는 ELSE구문이 없으므로, 모든 조건에 부합하지 못한 데이터가 만약 있다면 그 사람의 SMOK_S2는 결측치로 대체된다.

2) 방법 2 

*결측치가 없음을 확인하기;
PROC MEANS DATA=hong.df NMISS;
VAR SMOK;
RUN;

DATA hong.df_new;
SET hong.df;

*현재 흡연자 // 과거 흡연자 + 비흡연자로 나누기 (2);
IF SMOK=0 THEN SMOK_S2=0;
	ELSE IF SMOK=1 THEN SMOK_S2=0;
	ELSE SMOK_S2=1;

RUN;

이번에는 결측치가 없음을 확인하고 "ELSE IF"가 아닌 "ELSE"를 사용할 것이다. PROC MEANS를 이용하여 결측치가 없음을 우선 확인하자. 이에 관한 내용은 다음 링크에서 확인할 수 있다. 2022.09.23 - [기술 통계/SAS] - [SAS] 기술 통계 (평균, 표준편차, 표준오차, 최댓값, 최솟값, 중위수, 분위수 등) - PROC UNIVARIATE, PROC MEANS

 

결측치가 없음을 확인했다면, 다음의 사고 과정은 매우 논리적이다.

 

1) 1000명 전원은 {비흡연자, 과거 흡연자, 현재 흡연자} 셋 중의 하나에 반드시 속한다.

2) 비흡연자에게 SMOK_S2는 0의 값을 부여하고, 과거 흡연자에게 SMOK_S2는 1의 값을 부여한다.

3) 아직 SMOK_S2값을 부여받지 못한 사람은 모두 현재 흡연자이므로 SMOK_S2는 1의 값을 부여한다.

 위의 3) 내용을 코드로 옮기면 "ELSE SMOK_S2=1"가 된다. 

 

만약 결측치가 있다면 흡연 정보를 모르는 사람 또한 SMOK_S2의 값은 1이 부여되어 현재 흡연자로 평가받기 때문에 결측치가 있을 때에는 이 코드를 사용하면 안 된다.

 

3) 방법 3

*결측치가 없음을 확인하기;
PROC MEANS DATA=hong.df NMISS;
VAR SMOK;
RUN;

DATA hong.df_new;
SET hong.df;

*현재 흡연자 // 과거 흡연자 + 비흡연자로 나누기 (3);
IF SMOK=0 OR SMOK=1 THEN SMOK_S2=0;
	ELSE SMOK_S2=1;

RUN;

OR 연산자를 사용하면 코드가 조금 더 간단해진다.

OR 연산자는 OR 앞뒤로 있는 조건 중 어떤 것이라도 만족하면 옳다고(TRUE)로 인식한다. 우리는 과거 흡연자, 비흡연자 모두 SMOK_S2에서는 0이라는 값을 부여할 것이므로 SMOK변수는 0이든 1이든 둘 중 하나이기만 하면 된다. 

IF SMOK=0 OR SMOK=1 THEN SMOK_S2=0;

그래서 SMOK=0, SMOK=1 중 하나라도 만족하면 SMOK_S2에는 0이라는 값을 부여한다는 위 코드를 사용하였다. 

 

4) 방법 4

다음 링크에서 비교 연산자에 대해 공부하고 온다면 다음과 같이 쓸 수 있음도 알 수 있다. 2022.10.06 - [통계 프로그램 사용 방법/SAS] - [SAS] 변수 계산 (비교 연산) - DATA, SET, GE, LE, NE, IN, MIN, MAX

*결측치가 없음을 확인하기;
PROC MEANS DATA=hong.df NMISS;
VAR SMOK;
RUN;

DATA hong.df_new;
SET hong.df;

*현재 흡연자 // 과거 흡연자 + 비흡연자로 나누기 (4);
IF SMOK<=1 THEN SMOK_S2=0;
	ELSE SMOK_S2=1;

RUN;

과거 흡연자와 비흡연자의 SMOK값은 각각 0,1이므로 "SMOK<=1"이라고 쓸 수도 있다.

만약 SMOK에 결측치가 있다면 비교 연산자중 일부는 결측치를 0으로 인식하므로, "SMOK<=1"을 만족한 것으로 본다. 따라서 결측치가 없음을 미리 확인해야 하며, 이에 관한 내용은 다음 링크에서 확인할 수 있다. 2022.10.06 - [통계 프로그램 사용 방법/SAS] - [SAS] 변수 계산 (비교 연산) - DATA, SET, GE, LE, NE, IN, MIN, MAX

 

5) 방법 5, 6

비교 연산자 중 IN을 사용하여 다음과 같이 작성할 수도 있다.

*결측치가 없음을 확인하기;
PROC MEANS DATA=hong.df NMISS;
VAR SMOK;
RUN;

DATA hong.df_new;
SET hong.df;

*현재 흡연자 // 과거 흡연자 + 비흡연자로 나누기 (5);
IF SMOK IN (0 1) THEN SMOK_S2=0;
	ELSE SMOK_S2=1;

RUN;

DATA hong.df_new;
SET hong.df;

*현재 흡연자 // 과거 흡연자 + 비흡연자로 나누기 (6);
IF SMOK IN (0:1) THEN SMOK_S2=0;
	ELSE SMOK_S2=1;

RUN;

방법 5의 IF SMOK IN (0 1) THEN SMOK_S2=0;

- SMOK가 0, 1인 데이터에게만 SMOK_S2의 값으로 0을 부여한다는 뜻이다.

방법 6의 IF SMOK IN (0:1) THEN SMOK_S2=0;

- SMOK가 0 이상 1 이하의 정수인 데이터에게만 SMOK_S2의 값으로 0을 부여한다는 뜻이다.

 

6) 방법 7

방법 6은 사실 "SMOK가 0 이상이면서 1 이하인 정수 데이터"에게 SMOK_S2의 값으로 0을 부여한 것이다. 이는 AND연산자를 사용한 것과 다름없다. AND 연산자는 AND 앞뒤로 있는 조건을 모두 만족해야 TRUE를 반환하는데, 이 경우 SMOK>=0 AND SMOK<=1 이라는 문구를 사용하면 되는 것이다.

*결측치가 없음을 확인하기;
PROC MEANS DATA=hong.df NMISS;
VAR SMOK;
RUN;

DATA hong.df_new;
SET hong.df;

*현재 흡연자 // 과거 흡연자 + 비흡연자로 나누기 (7);
IF SMOK>=0 AND SMOK<=1 THEN SMOK_S2=0;
	ELSE SMOK_S2=1;

RUN;

 

AND 연산자는 OR 연산자보다 우선하므로, 괄호를 적절히 치며 사용해야 한다. (사칙연산에서 곱셈이 덧셈보다 우선하는 것과 동일하다)

 

SAS 변수 계산 (논리 연산) 정복 완료!

작성일: 2022.10.06.

최종 수정일: 2022.11.03.

이용 프로그램: SAS v9.4

운영체제: Windows 10

 

반응형

+ Recent posts