[R 통계분석] 기술통계학
> ## 00) data 준비
>
> # 통계청 마이크로데이터 통합서비스(https://mdis.kostat.go.kr/) - 인구 > 인구주택총조사 > 1%_인구사항(제공) > 2010
> data <- read.csv("./data/ch02.csv", na.strings = c("."))
> str(data)
'data.frame': 468284 obs. of 6 variables:
$ 성별코드 : int 1 1 1 1 1 1 1 1 1 1 ...
$ 만연령 : int 0 0 0 0 0 0 0 0 0 0 ...
$ 가구주관계코드 : int 3 3 3 3 3 3 3 3 3 3 ...
$ 교육정도코드 : int 1 1 1 1 1 1 1 1 1 1 ...
$ 총출생아수_남아수: int NA NA NA NA NA NA NA NA NA NA ...
$ 총출생아수_여아수: int NA NA NA NA NA NA NA NA NA NA ...
> unique(data$성별코드)
[1] 1 2
> data$성별코드 <- factor(data$성별코드, levels = c(1, 2),
+ labels = c("남자","여자"))
> sort(unique(data$가구주관계코드))
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14
> data$가구주관계코드 <- factor(data$가구주관계코드, levels=1:14,
+ labels = c("가구주", "가구주의 배우자", "자녀",
+ "자녀의 배우자", "가구주의 부모",
+ "배우자의 부모", "손자녀, 그 배우자",
+ "증손자녀, 그 배우자", "조부모",
+ "형제자매, 그 배우자",
+ "형제자매의 자녀, 그 배우자",
+ "부모의 형제자매, 그 배우자", "기타 친인척",
+ "그외같이사는사람"))
> sort(unique(data$교육정도코드))
[1] 1 2 3 4 5 6 7 8
> data$교육정도코드 <- factor(data$교육정도코드, levels=1:8,
+ labels = c("안 받았음", "초등학교", "중학교",
+ "고등학교", "대학-4년제 미만", "대학-4년제 이상",
+ "석사과정", "박사과정"))
> str(data)
'data.frame': 468284 obs. of 6 variables:
$ 성별코드 : Factor w/ 2 levels "남자","여자": 1 1 1 1 1 1 1 1 1 1 ...
$ 만연령 : int 0 0 0 0 0 0 0 0 0 0 ...
$ 가구주관계코드 : Factor w/ 14 levels "가구주","가구주의 배우자",..: 3 3 3 3 3 3 3 3 3 3 ...
$ 교육정도코드 : Factor w/ 8 levels "안 받았음","초등학교",..: 1 1 1 1 1 1 1 1 1 1 ...
$ 총출생아수_남아수: int NA NA NA NA NA NA NA NA NA NA ...
$ 총출생아수_여아수: int NA NA NA NA NA NA NA NA NA NA ...
> colnames(data)[c(1,3,4)]
[1] "성별코드" "가구주관계코드" "교육정도코드"
> colnames(data)[c(1,3,4)] <- c("성별", "가구주관계", "교육정도")
> str(data)
'data.frame': 468284 obs. of 6 variables:
$ 성별 : Factor w/ 2 levels "남자","여자": 1 1 1 1 1 1 1 1 1 1 ...
$ 만연령 : int 0 0 0 0 0 0 0 0 0 0 ...
$ 가구주관계 : Factor w/ 14 levels "가구주","가구주의 배우자",..: 3 3 3 3 3 3 3 3 3 3 ...
$ 교육정도 : Factor w/ 8 levels "안 받았음","초등학교",..: 1 1 1 1 1 1 1 1 1 1 ...
$ 총출생아수_남아수: int NA NA NA NA NA NA NA NA NA NA ...
$ 총출생아수_여아수: int NA NA NA NA NA NA NA NA NA NA ...
> save.image("data.rda")
> summary(data)
성별 만연령 가구주관계 교육정도 총출생아수_남아수 총출생아수_여아수
남자:226965 Min. : 0.00 가구주 :179293 고등학교 :134246 Min. : 0.00 Min. : 0.00
여자:241319 1st Qu.:22.00 자녀 :145533 대학-4년제 이상: 81110 1st Qu.: 1.00 1st Qu.: 0.00
Median :40.00 가구주의 배우자 :106311 초등학교 : 80710 Median : 1.00 Median : 1.00
Mean :39.34 가구주의 부모 : 12069 중학교 : 55704 Mean : 1.32 Mean : 1.23
3rd Qu.:55.00 손자녀, 그 배우자 : 7832 안 받았음 : 51085 3rd Qu.: 2.00 3rd Qu.: 2.00
Max. :85.00 형제자매, 그 배우자: 5332 대학-4년제 미만: 50753 Max. :12.00 Max. :11.00
(Other) : 11914 (Other) : 14676 NA's :310308 NA's :310308
> head(subset(data, subset = is.na(총출생아수_남아수) & !is.na(총출생아수_여아수)))
[1] 성별 만연령 가구주관계 교육정도 총출생아수_남아수 총출생아수_여아수
<0 행> <또는 row.names의 길이가 0입니다>
> head(subset(data, subset = !is.na(총출생아수_남아수) & is.na(총출생아수_여아수)))
[1] 성별 만연령 가구주관계 교육정도 총출생아수_남아수 총출생아수_여아수
<0 행> <또는 row.names의 길이가 0입니다>
>
Section 01. 그래프
> ## 01) 두 변수 간의 관계를 나타내는 산점도
>
> str(cars)
'data.frame': 50 obs. of 2 variables:
$ speed: num 4 4 7 7 8 9 10 10 10 11 ...
$ dist : num 2 10 4 22 16 10 18 26 34 17 ...
> plot(cars$speed, cars$dist,
+ main = "속도와 제동거리", xlab = "속도(mph)", ylab = "제동거리(ft)",
+ pch = 1, col = "red")
>
> Nile
Time Series:
Start = 1871
End = 1970
Frequency = 1
[1] 1120 1160 963 1210 1160 1160 813 1230 1370 1140 995 935 1110 994 1020 960 1180 799 958 1140 1100 1210 1150
[24] 1250 1260 1220 1030 1100 774 840 874 694 940 833 701 916 692 1020 1050 969 831 726 456 824 702 1120
[47] 1100 832 764 821 768 845 864 862 698 845 744 796 1040 759 781 865 845 944 984 897 822 1010 771
[70] 676 649 846 812 742 801 1040 860 874 848 890 744 749 838 1050 918 986 797 923 975 815 1020 906
[93] 901 1170 912 746 919 718 714 740
> str(Nile)
Time-Series [1:100] from 1871 to 1970: 1120 1160 963 1210 1160 1160 813 1230 1370 1140 ...
> plot(Nile,
+ main = "Nile강의 연도별 유량 변화", xlab = "연도", ylab = "유량")
> plot(Nile, type = "p",
+ main = "Nile강의 연도별 유량 변화", xlab = "연도", ylab = "유량")
>
> ## 02) 막대그래프와 히스토그램
>
> load("data.rda")
> tableBabyM <- table(data$총출생아수_남아수)
> tableBabyM
0 1 2 3 4 5 6 7 8 9 10 11 12
30788 69624 41010 11165 3667 1228 346 104 21 8 4 10 1
> barplot(tableBabyM,
+ main = "출생아(남자)별 빈도", xlab = "출생아수", ylab = "빈도")
>
> tableSexEdu <- table(data$성별, data$교육정도)
> tableSexEdu
안 받았음 초등학교 중학교 고등학교 대학-4년제 미만 대학-4년제 이상 석사과정 박사과정
남자 19161 34214 26588 66548 25673 45530 7107 2144
여자 31924 46496 29116 67698 25080 35580 4634 791
> barplot(tableSexEdu, legend.text = T, col = c("orange", "green"),
+ main = "학력에 따른 성별 인원수", xlab = "학력", ylab = "빈도")
> barplot(tableSexEdu, legend.text = T, col = c("orange", "green"), beside = T,
+ main = "학력에 따른 성별 인원수", xlab = "학력", ylab = "빈도")
>
> hist(data$만연령, main = "연령별 분포", xlab = "연령", ylab = "빈도")
> hist(data$만연령, breaks = c(seq(0, 90, 10)), right = F,
+ main = "연령별 분포", xlab = "연령", ylab = "빈도")
>
> ## 03) 원도표
>
> load("data.rda")
> tableEdt <- table(data$교육정도)
> tableEdt
안 받았음 초등학교 중학교 고등학교 대학-4년제 미만 대학-4년제 이상 석사과정
51085 80710 55704 134246 50753 81110 11741
박사과정
2935
> pie(tableEdt, main = "학력수준별 비중")
>
Section 02. 모수와 통계량
> ranicafe <- read.csv("./data/cafedata.csv", header=T, na.strings="na", stringsAsFactors=FALSE )
> ranicafe <- na.omit(ranicafe)
> str(ranicafe)
'data.frame': 47 obs. of 22 variables:
$ t : int 1 2 3 4 5 6 7 8 9 10 ...
$ Date : chr "2010-01-19" "2010-01-20" "2010-01-21" "2010-01-22" ...
$ Day.Code : int 2 3 4 5 1 2 3 4 5 1 ...
$ Day.of.Week : chr "Tue" "Wed" "Thu" "Fri" ...
$ Bread.Sand.Sold : int 5 6 8 4 3 7 6 0 3 2 ...
$ Bread.Sand.Waste : int 3 8 2 2 0 1 6 0 4 6 ...
$ Wraps.Sold : int 25 7 14 5 10 5 19 7 4 13 ...
$ Wraps.Waste : int 5 17 0 7 0 3 3 0 9 3 ...
$ Muffins.Sold : int 5 3 4 5 8 1 6 6 0 3 ...
$ Muffins.Waste : int 1 5 0 0 0 0 0 1 4 0 ...
$ Cookies.Sold : int 5 1 1 3 3 5 10 0 3 6 ...
$ Cookies.Waste : int 3 6 0 1 0 0 0 0 2 0 ...
$ Fruit.Cup.Sold : int 1 0 0 3 2 2 2 0 1 2 ...
$ Fruit.Cup.Waste : int 4 3 3 0 0 0 0 0 1 0 ...
$ Chips : int 12 0 0 20 0 4 2 20 3 16 ...
$ Juices : int 8 0 13 0 5 4 5 6 4 7 ...
$ Sodas : int 20 13 23 13 13 33 15 27 12 19 ...
$ Coffees : int 41 33 34 27 20 23 32 31 30 27 ...
$ Total.Soda.and.Coffee : int 61 46 57 40 33 56 47 58 42 46 ...
$ Sales : num 200 196 103 163 102 ...
$ Max.Daily.Temperature..F.: int 36 34 39 40 36 26 34 33 20 37 ...
$ Total.Items.Wasted : int 16 39 5 10 0 4 9 1 20 9 ...
- attr(*, "na.action")= 'omit' Named int 34
..- attr(*, "names")= chr "34"
> head(ranicafe)
t Date Day.Code Day.of.Week Bread.Sand.Sold Bread.Sand.Waste Wraps.Sold Wraps.Waste Muffins.Sold Muffins.Waste
1 1 2010-01-19 2 Tue 5 3 25 5 5 1
2 2 2010-01-20 3 Wed 6 8 7 17 3 5
3 3 2010-01-21 4 Thu 8 2 14 0 4 0
4 4 2010-01-22 5 Fri 4 2 5 7 5 0
5 5 2010-01-25 1 Mon 3 0 10 0 8 0
6 6 2010-01-26 2 Tue 7 1 5 3 1 0
Cookies.Sold Cookies.Waste Fruit.Cup.Sold Fruit.Cup.Waste Chips Juices Sodas Coffees Total.Soda.and.Coffee Sales
1 5 3 1 4 12 8 20 41 61 199.95
2 1 6 0 3 0 0 13 33 46 195.74
3 1 0 0 3 0 13 23 34 57 102.68
4 3 1 3 0 20 0 13 27 40 162.88
5 3 0 2 0 0 5 13 20 33 101.76
6 5 0 2 0 4 4 33 23 56 186.94
Max.Daily.Temperature..F. Total.Items.Wasted
1 36 16
2 34 39
3 39 5
4 40 10
5 36 0
6 26 4
>
> library(ggplot2)
> ggplot(ranicafe, aes(Coffees)) +
+ geom_bar(fill="gray") +
+ ggtitle("라니의 카페 커피 판매량") +
+ theme(plot.title = element_text(size = 20, face="bold")) +
+ xlim(0, 50) + xlab("판매량") +
+ ylab("횟수") + scale_y_continuous(breaks=0:10)
>
> ## 04) 최대값과 최소값
>
> sort(ranicafe$Coffees)[1]
[1] 3
> sort(ranicafe$Coffees, decreasing = T)[1]
[1] 48
> min(ranicafe$Coffees)
[1] 3
> max(ranicafe$Coffees)
[1] 48
>
> # 최빈값
> rc <- ranicafe$Coffees
> stem(rc)
The decimal point is 1 digit(s) to the right of the |
0 | 34444
0 | 5688
1 | 01134
1 | 668
2 | 001123344
2 | 55677789
3 | 001112334
3 | 55
4 | 1
4 | 8
> #-> 4가 4번 등장하여 최빈값은 4
> tableRc <- table(rc)
> names(tableRc)[which(tableRc == max(tableRc))]
[1] "4"
>
> ## 05) 평균
> weight <- 1 / length(rc)
> sum(rc * weight)
[1] 21.51064
> mean(rc)
[1] 21.51064
>
> ## 07) 중앙값
> (median.idx <- (length(rc)+1) / 2)
[1] 24
> (rc.srt <- sort(rc))
[1] 3 4 4 4 4 5 6 8 8 10 11 11 13 14 16 16 18 20 20 21 21 22 23 23 24 24 25 25 26 27 27 27 28 29 30 30 31 31 31
[40] 32 33 33 34 35 35 41 48
> rc.srt[median.idx]
[1] 23
> median(rc)
[1] 23
>
> ## 09) 개별 관찰값과 평균과의 차이에 대한 평균
> height <- c(164, 166, 168, 170, 172, 174, 176)
> (height.m <- mean(height))
[1] 170
> (height.dev <- height - height.m)
[1] -6 -4 -2 0 2 4 6
> sum(height.dev)
[1] 0
>
> ## 10) 편차 제곱의 평균 구하기
> (height.dev2 <- height.dev ^ 2)
[1] 36 16 4 0 4 16 36
> sum(height.dev2)
[1] 112
> mean(height.dev2)
[1] 16
>
> ## 11) 표준편차 구하기
> sqrt(mean(height.dev2))
[1] 4
>
> ## 12) 분산과 표준편차 구하기
> var(height)
[1] 18.66667
> sd(height)
[1] 4.320494
> # R은 대부분 표본조사라고 판단하고 갯수(n)가 아닌 자유도(n-1)로 나눈다.
>
> ## 13) '라나의 카페' 커피 판매량에 대한 평균과 표준편차
> rc.m <- mean(rc)
> rc.sd <- sd(rc)
> cat("커피 판매량", round(rc.m, 1), "±", round(rc.sd, 2), "잔")
커피 판매량 21.5 ± 11.08 잔
>
> # 정규분포
> par(mar=c(2, 0, 0, 0))
> x <- seq(-3, 3, by=0.001)
> plot(x, dnorm(x), type="l", lwd=2, xlab="", ylab="", axes=F)
> axis(1)
> lines(c(0, 0), c(0, dnorm(0)), lty=1, lwd=1.5)
> lines(c(-1, -1), c(0, dnorm(-1)), lty=3, lwd=2)
> lines(c(1, 1), c(0, dnorm(1)), lty=3, lwd=2)
>
> # 표준편차가 다른 두 자료
> par(mar=c(2, 0, 0, 0))
> x <- seq(158, 182, by=0.01)
> plot(x, dnorm(x, mean=170, sd=4), type="l", lwd=2, ylim=c(0, 0.2), axes=F)
> lines(x, dnorm(x, mean=170, sd=2), lty=2, lwd=2)
> axis(1)
> legend("topright",
+ lwd=2,
+ lty=c(1,2),
+ c(paste("sd :", 4), paste("sd :", 2))
+ )
>
> ## 14) 변동계수 구하기
> rj <- ranicafe$Juices
> (rj.m <- mean(rj))
[1] 4.93617
> (rj.sd <- sd(rj))
[1] 3.703138
> (rc.cv <- round(rc.sd / rc.m, 3))
[1] 0.515
> (rj.cv <- round(rj.sd / rj.m, 3))
[1] 0.75
>
> ## 15) 사분위수 범위와 상자도표
> (qs <- quantile(rc))
0% 25% 50% 75% 100%
3 12 23 30 48
> (qs[4] - qs[2])
75%
18
> IQR(rc)
[1] 18
> bp <- boxplot(rc, main = "커피 판매량에 대한 상자도표")
> axis(2)
> md <- median(rc)
> arrows(1.4, qs[4], 1.22, qs[4], angle=30, length=0.1)
> text(1.42, qs[4], "Q3", adj=0)
> arrows(1.4, md, 1.22, md, angle=30, length=0.1)
> text(1.42, md, "Q2", adj=0)
> arrows(1.4, qs[2], 1.22, qs[2], angle=30, length=0.1)
> text(1.42, qs[2], "Q1", adj=0)
>
> # 제동거리의 히스토그램
> hist(cars$dist, breaks = seq(0, 120, 10))
>
> ## 16) 이상치 판별
> ( Q <- quantile(cars$dist))
0% 25% 50% 75% 100%
2 26 36 56 120
> (ll <- Q[2] - 1.5*IQR(cars$dist))
25%
-19
> (ul <- Q[4] + 1.5*IQR(cars$dist))
75%
101
> cars$dist[cars$dist < ll]
numeric(0)
> cars$dist[cars$dist > ul]
[1] 120
> boxplot(cars$dist, main = "Boxplot of Distance")
출처 : 이윤환, ⌜제대로 알고 쓰는 R 통계분석⌟, 한빛아카데미, 2020