[R] 기술 통계 (평균, 표준편차, 표준오차, 최댓값, 최솟값, 중위수, 분위수 등)
1,000명으로 어떤 연구를 했다고 하자. 그들의 키, 몸무게 등 지표들은 서로 다를 것이다. 논문의 저자가 이 모든 것을 독자들에게 보여주고자 한다면 행이 1,000인 표를 제시해야 할 것이다. 그렇게 큰 표를 실어줄 저널이 없기도 하거니와, 독자들이 보기에도 한눈에 들어오지 않는다. 그 대신 키의 '평균', 몸무게의 '평균'을 제시하면 한눈에 들어오니 보기가 좋다. 연속 변수는 평균, 표준편차 등으로 요약을 하여 보여주고, 범주형 자료 (흡연 여부, 음주 여부 등)는 도수분포표 혹은 분할표로 제시하게 된다. 분할표를 작성하는 방법은 다음 링크에서 확인할 수 있다.
2022.08.31 - [기술 통계/R] - [R] 도수분포표 (Frequency table), 분할표 (Contingency table) 만들기 - CrossTable()
어떤 지표로 요약해줄 것인가?
1) 정규성을 따를 때: 평균 및 표준편차
어떤 변수가 정규 분포를 따른다고 할 수 있다면, 평균과 표준편차만 알면 된다. 단 두 개의 지표만 있으면 전체 분포를 알아낼 수 있기 때문이다. (정규성을 따르는지는 정규성 검정으로 확인할 수 있으며 정규성 검정을 하는 방법은 다음 링크에서 확인할 수 있다.)
2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (1) : Q-Q plot - qqnorm(), qqline()
2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (2) : 히스토그램 - hist(), dnorm()
2022.08.12 - [기술 통계/R] - [R] 정규성 검정 (4) : 정량적 검정 (Lilliefors test) - lillie.test()
2) 정규성을 따르지 않을 때: 중위수, 최댓값, 최솟값, 분위수, 사분위 범위 등
정규성을 따르지 않는다면 평균과 표준편차를 안다고 해도 전체 분포를 알아낼 수는 없다. 따라서 분포에 대한 직접적인 정보를 주는데, 예를 들어 '하위 25%에 위치하는 사람의 ALT값은 얼마인가?' 등을 제시하는 것이다. 그런 지표로는 중위수, 최댓값, 최솟값, 분 위수, 사분위 범위 등이 있다.
이번 포스팅에서는 이 모든 지표들 (평균, 표준편차, 표준오차, 중위수, 최댓값, 최솟값, 분위수, 사분위 범위 등)을 구하는 법에 대해 소개할 것이다.
*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.
2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.08.29)
분석용 데이터 (update 22.08.29)
2022년 08월 29일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 통계 프로그램 사용 방법 1) 엑셀 파일 2) CSV 파일 3) 코드북
medistat.tistory.com
코드를 보여드리기에 앞서 워킹 디렉토리부터 지정하겠다.
워킹 디렉토리에 관한 설명은 다음 링크된 포스트에서 볼 수 있다.
2022.08.05 - [통계 프로그램 사용 방법/SAS] - [SAS] 라이브러리 만들기 - LIBNAME
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")
코드
이번 포스팅을 일일이 이해하기는 귀찮고 코드만 가져가고 싶은 사람은 다음 코드를 사용하면 된다. :
#기술통계 함수 작성
univariate<-function(x){
n<-length(x)-sum(is.na(x))
missing<-sum(is.na(x))
mean<-mean(x)
sd<-sd(x)
var<-var(x)
se<-sqrt(var/n)
min<-min(x)
max<-max(x)
median<-median(x)
mode <- function(x) {
nodup <- unique(x)
if (length(which(tabulate(match(x, nodup))==max(tabulate(match(x, nodup)))))>=2 )
{NA}
else{nodup[which.max(tabulate(match(x, nodup)))]}
}
mode<-mode(x)
lclm<-mean-qt(0.975, n-1)*se
uclm<-mean+qt(0.975, n-1)*se
P25<-data.frame(quantile(x)[2])[1,1]
P50<-data.frame(quantile(x)[3])[1,1]
P75<-data.frame(quantile(x)[4])[1,1]
qrange<-data.frame(quantile(x)[4])[1,1]-data.frame(quantile(x)[2])[1,1]
name<-c("N", "N of missing","Mean", "Standard Deviation", "Variation", "Standard Error", "Minimum",
"Maximum", "Median", "Mode", "95% Lower limit of mean (two-sided)", "95% Upper limit of mean(two-sided)", "P25", "P50", "P75", "Quantile Range")
value<-c(n, missing, mean, sd, var, se, min, max, median, mode, lclm, uclm, P25, P50, P75, qrange)
result<-cbind(name, value)
print(result)
}
#기술 통계량 계산: "df$ALT"자리에 원하는 변수를 넣으면 된다.
univariate(df$ALT)
결과
name value
[1,] "N" "1000"
[2,] "N of missing" "0"
[3,] "Mean" "35.11593830339"
[4,] "Standard Deviation" "5.00643005419849"
[5,] "Variation" "25.0643418875819"
[6,] "Standard Error" "0.158317219175875"
[7,] "Minimum" "19.61781834"
[8,] "Maximum" "49.03881459"
[9,] "Median" "35.13945912"
[10,] "Mode" NA
[11,] "95% Lower limit of mean (two-sided)" "34.8052658601898"
[12,] "95% Upper limit of mean(two-sided)" "35.4266107465902"
[13,] "P25" "31.7559160525"
[14,] "P50" "35.13945912"
[15,] "P75" "38.5879822725"
[16,] "Quantile Range" "6.83206622""N" : 분석 대상의 수
"N of missing" : 결측치 수
"Mean" : 평균
"Standard Deviation" : 표준 편차
"Variation" : 분산
"Standard Error" : 표준 오차
"Minimum" : 최솟값
"Maximum" : 최댓값
"Median" : 중위수
"Mode" : 최빈값 (최빈값이 여러개인 경우 NA로 표기됨)
"95% Lower limit of mean (two-sided)" : 평균의 95% 신뢰구간의 하한
"95% Upper limit of mean(two-sided)" : 평균의 95% 신뢰구간의 상한
"P25" : 1사분위수
"P50" : 2사분위수 (중위수)
"P75" : 3사분위수
"Quantile Range" : 사분위 범위
코드 (분석 대상의 수, 결측치 수, 평균, 표준 편차, 분산, 표준 오차)
코드를 하나씩 뜯어 살펴보면 다음과 같다.
univariate<-function(x){
n<-length(x)-sum(is.na(x))
missing<-sum(is.na(x))
mean<-mean(x)
sd<-sd(x)
var<-var(x)
se<-sqrt(var/n)
}univariate<-function(x){ :입력되는 값이 x인 함수를 만들건데 함수의 내용은 { }안에 들어있으니 살펴봐라. 그리고 함수의 이름은 univariate이다.
n<-length(x)-sum(is.na(x)) : 분석 대상의 수[n]는 전체 개체의 수[length(x)]에서 결측치의 개수 [sum(is.na(x))]를 뺀 값이다.
missing<-sum(is.na(x)) : 결측된 것을 TRUE로, 결측되지 않은 것을 FALSE로 반환한다 [is.na(x)]. TRUE인 것을 모두 합친다 [sum(is.na(x))]. 그 값이 결측치의 개수이다.
mean<-mean(x) : 평균은 mean() 함수를 이용한다.
sd<-sd(x) : 표준편차는 sd() 함수를 이용한다.
var<-var(x) : 분산은 var() 함수를 이용한다.
se<-sqrt(var/n) : 표준오차(se)는 분산(var)을 분석 대상의 수(n)로 나눈 뒤 루트를 씌운 값이다.
min<-min(x)
max<-max(x)
median<-median(x)
mode <- function(x) {
nodup <- unique(x)
if (length(which(tabulate(match(x, nodup))==max(tabulate(match(x, nodup)))))>=2 )
{NA}
else{nodup[which.max(tabulate(match(x, nodup)))]}
}
mode<-mode(x) min<-min(x) : 최솟값은 min() 함수를 이용한다.
max<-max(x) : 최댓값은 max() 함수를 이용한다.
median<-median(x) : 중앙값은 median() 함수를 이용한다.
mode <- function(x) { : 최빈값은 구하는 함수가 R에 내장되어있지 않으므로 함수를 따로 만들어야 한다. 데이터(x)를 입력받아 최빈값을 구하는 함수를 만들 것이고 그 함수의 이름은 mode이다.
nodup <- unique(x) : 중복되는 값을 모두 삭제하고 한 개씩만 대표적으로 남긴다. 대표들로 이루어진 벡터를 nodup이라고 하자.
if (length(which(tabulate(match(x, nodup))==max(tabulate(match(x, nodup)))))>=2 )
{NA}
else{nodup[which.max(tabulate(match(x, nodup)))]}
여기까지의 내용이 복잡하니 하나씩 코드를 살펴보기로 한다.
(1) 큰 구성은 if (A) {B} else {C} 이다. A라는 조건이 만족하면 B를 반환하고, A를 만족하지 않으면 C를 반환한다는 것이다.
(2) match(x, nodup): x가 nodup에서 몇 번째 데이터와 일치하는지를 반환하는 함수다.
Example) 예를 들어 x가 (1,5,2,3,3,4,2)라고 하자.
nodup은 중복되는 값을 삭제했을 터이니 (1,5,2,3,4)일 것이다.
그렇다면 match(x,nodup)은 (1,2,3,4,4,5,3)이 된다.
(3) tabulate(match(x,nodup)): tabulate는 각 정수가 나타난 횟수를 세준다. (1,2,3,4,4,5,3)에서 1은 1개, 2는 1개, 3은 2개, 4는 2개, 5는 1개다. 따라서 tabulate(match(x,nodup))은 (1,1,2,2,1)이 된다.
최빈값은 가장 흔한 단 한 개의 값이어야 하는데, 여기에서 3과 4는 2개, 나머지는 1개이므로 이 데이터의 최빈값은 없다. 즉, tabulate(match(x,nodup))인 (1,1,2,2,1)의 최댓값(2)이 tabulate(match(x,nodup))인 (1,1,2,2,1)에 2개 이상 존재한다면 최빈값이 없다고 반환해야 한다.
(4) tabulate(match(x, nodup))==max(tabulate(match(x, nodup))): tabulate(match(x,nodup))인 (1,1,2,2,1)이 tabulate(match(x,nodup))인 (1,1,2,2,1)의 최댓값(2)과 같은지를 확인해주는 코드다 결과는 (FALSE, FALSE, TRUE, TRUE, FALSE)이다.
(5) which(tabulate(match(x, nodup))==max(tabulate(match(x, nodup))))는 TRUE인 위치를 반환하므로 (3,4)를 보여준다.
(6) length(which(tabulate(match(x, nodup))==max(tabulate(match(x, nodup)))))는 (3,4)에 있는 숫자의 개수, 즉 2를 반환한다.
이렇게 함으로써 최빈값'같은' 값이 2개 이상 존재할 때 NA를 반환하도록 하였다.
그렇지 않을 때 (즉, 단 한 개의 최빈값이 존재할 때)는 다음과 같은 로직에 의해 구해진다고 생각하면 된다.
Example) x가 (1,5,2,2,3)이라고 하자.
nodup은 (1,5,2,3)이 된다.
match(x,nodup)은 (1,2,3,3,4)이 된다.
tabulate(match(x,nodup))는 (1,1,2,1)이 된다. 여기에서 횟수가 가장 많아 '2'라고 표현된 값의 근원을 찾아가야 한다.
which.max(tabulate(match(x, nodup)))는 (1,1,2,1)에서 최댓값이 존재하는 위치를 반환하므로 3이 된다.
(1,1,2,1)의 2는 match(x,nodup)인 (1,2,3,3,4)에서 온 것이다. 1이 1개, 2가 1개, 3이 2개, 4가 1개. 즉, tabulate는 개수를 세어주므로 중복된 값을 빼는 효과가 있다. 즉 2가 3번째에 위치한 이유는 nodup (1,5,2,3)의 2가 nodup에서 3번째에 위치해있기 때문이다. 따라서 nodup의 3번째 값을 반환해야 최빈값을 구할 수 있다. 그러므로 nodup[which.max(tabulate(match(x, nodup)))]이라는 코드를 쓰게 된다.
}
mode<-mode(x)
lclm<-mean-qt(0.975, n-1)*se
uclm<-mean+qt(0.975, n-1)*se
P25<-data.frame(quantile(x)[2])[1,1]
P50<-data.frame(quantile(x)[3])[1,1]
P75<-data.frame(quantile(x)[4])[1,1]
qrange<-data.frame(quantile(x)[4])[1,1]-data.frame(quantile(x)[2])[1,1]
name<-c("N", "N of missing","Mean", "Standard Deviation", "Variation", "Standard Error", "Minimum",
"Maximum", "Median", "Mode", "95% Lower limit of mean (two-sided)", "95% Upper limit of mean(two-sided)", "P25", "P50", "P75", "Quantile Range")
value<-c(n, missing, mean, sd, var, se, min, max, median, mode, lclm, uclm, P25, P50, P75, qrange)
result<-cbind(name, value)
print(result)
}
lclm<-mean-qt(0.975, n-1)*se : t 분포에서 자유도(n-1)에 걸맞게 2.5%에 해당하는 값을 찾고, 표준오차를 곱하여 평균에서 빼면 95% 신뢰구간의 하한이 구해진다.
uclm<-mean+qt(0.975, n-1)*se : t 분포에서 자유도(n-1)에 걸맞게 2.5%에 해당하는 값을 찾고, 표준오차를 곱하여 평균에서 더하면 95% 신뢰구간의 상한이 구해진다.
P25<-data.frame(quantile(x)[2])[1,1] : quantile() 함수의 2번째 값은 1사분위수를 의미한다. 데이터프레임으로 바꾸고 값만을 취한다.
P50<-data.frame(quantile(x)[3])[1,1] : quantile() 함수의 2번째 값은 2사분위수를 의미한다. 데이터프레임으로 바꾸고 값만을 취한다.
P75<-data.frame(quantile(x)[4])[1,1] : quantile() 함수의 2번째 값은 3사분위수를 의미한다. 데이터프레임으로 바꾸고 값만을 취한다.
qrange<-data.frame(quantile(x)[4])[1,1]-data.frame(quantile(x)[2])[1,1] : 3사분위수에서 1사분위수를 빼어 값을 구한다.
name<-c("N", "N of missing","Mean", "Standard Deviation", "Variation", "Standard Error", "Minimum",
"Maximum", "Median", "Mode", "95% Lower limit of mean (two-sided)", "95% Upper limit of mean(two-sided)", "P25", "P50", "P75", "Quantile Range") : 각 값의 이름을 지정하여 name이라는 벡터에 저장한다.
value<-c(n, missing, mean, sd, var, se, min, max, median, mode, lclm, uclm, P25, P50, P75, qrange) : 위에서 구한 값들로 벡터를 만들고 value라는 이름을 붙인다.
result<-cbind(name, value) : name과 value를 열 별로 합치고 result에 저장한다.
print(result) : result를 출력한다.
}
만약 25%, 50%, 75%가 아닌 임의의 n% 백분위수를 구하고 싶다면 Pn을 사용하면 된다.
64%백분위수를 보고싶다면 P64를 사용하면 된다.
사분위수는 SAS, SPSS, R의 결과가 서로 다를 수 있다. 왜냐하면 각 프로그램에서 사분위수를 구하는 방법이 다를 수 있기 때문이다. 각 프로그램에는 사분위수를 구하는 여러가지 방법이 내장되어 있으며, 골라서 사용할 수도 있다.
SAS의 사분위수 확인하기:2022.09.23 - [기술 통계/SAS] - [SAS] 기술 통계 (평균, 표준편차, 표준오차, 최댓값, 최솟값, 중위수, 분위수 등) - PROC UNIVARIATE, PROC MEANS
SPSS의 사분위수 확인하기: 2022.09.29 - [기술 통계/SPSS] - [SPSS] 기술 통계 (평균, 표준편차, 표준오차, 최댓값, 최솟값, 중위수, 분위수 등)
R 기술 통계량 정복 완료!
작성일: 2022.09.27.
최종 수정일: 2022.09.29.
이용 프로그램: R 4.1.3
RStudio v1.4.1717
RStudio 2021.09.1+372 "Ghost Orchid" Release
운영체제: Windows 10, Mac OS 10.15.7