brunch

You can make anything
by writing

C.S.Lewis

by 쩐시 Aug 20. 2021

R 프로그래밍 - 2021 서울대기오염 분석 4(최종)

데이터 분석

앞서 데이터 전처리까지 끝난 seoulair 데이터를 여러가지 방면으로 분석해보도록 한다.

참고로, soeulair 데이터는 아래와 같다.



최종 분석 결과는 " 미세먼지와 초미세먼지는 서로 영향력이 있으며, 비례관계를 가진다 "


<분석 목록>

1. 연간 미세먼지 평균

2. 미세먼지가 가장 심했던 날짜

3. 구별 미세먼지 평균 비교

4. 계절별 분석

5. 미세먼지 등급 분석

6. 구별 미세먼지 등급 비교

7. 1년간 미세먼지 추이를 그래프로 그리기

8. 계절별 미세먼지 등급 비율 그래프 그리기

9. 25개 구의 연간 초미세먼지 평균 누운 막대그래프로 그리기.

10. 연간 초미세먼지 평균이 가장 높은 구와 가장 낮은 구의 평균 t-test+cor.test,linear reggresion 검정





1. 연간 미세먼지 평균


현재 테이블에 결측치를 제거하지 않은 상태기 때문에, 아래와 같이 입력해준다. 


mean(seoulair$pm10, na.rm = T)  "38.67256, 보통 수준" 

mean(!is.na(seoulair$pm10)) -> 이건 시험삼아 해본건데 error가 나니 조심!


뒤에 더 상세한 평균 활용 분석은 group_by&summarise를 통해 진행한다.


2. 미세먼지가 가장 심했던 날짜


seoulair %>% 

  filter(!is.na(pm10)) %>% 

  filter(pm10 == max(pm10)) %>% 

  select(date, district, pm10) 

금년도 5월 강동구의 미세먼지 수준은 경악할 정도로 높다!

코드가 잘못된 건가 싶어 홈페이지에서 직접 찾아보니 사실... 무슨일이 있었던 건지...


3. 구별 미세먼지 평균 비교


서울시 25개 구의 1년간 미세먼지 평균을 구한 후, 평균이 가장 낮은 순& 가장 높은 순 5개 구의 이름과 평균만 출력한다.


1) 평균이 가장 낮은 5개 구(오름차순)

seoulair %>% 

  filter(!is.na(pm10)) %>% 

  group_by(district) %>% 

  summarise(m = mean(pm10)) %>% 

  arrange(m) %>% 

  head(5)

상당히 의외였다, 서대문구 같은 경우는 주변에 중구, 동대문 다 있어서 교통이 굉장히 혼잡한데, 미센먼지의 평균이 가장 낮았다!

뭐... 물론 평균이지만, 대학가가 대부분이라 그런걸까 싶다. 코로나로 인해 비대면 개강을 해서 그런가...


1) 평균이 가장 높은 5개 구(내림차순)

seoulair %>% 

  filter(!is.na(pm10)) %>% 

  group_by(district) %>% 

  summarise(m = mean(pm10)) %>% 

  arrange(desc(m)) %>% 

  head(5)

예상했던대로, 미세먼지 평균이 가장 안좋은 곳은 너무나 교통이 혼잡한 중구였다.

언제가도 하루종일 직장인들로 혼잡한 그곳... 코로나와 상관없이 출근률이 매번 80%, 서울의 중심이니 당연한것일까 생각이 든다.


4. 계절별 분석


미세먼지와 초미세먼지 모두 결측치를 제거한 값을 구하도록 한다. 계절별 각 먼지들의 평균값을 추출하고, 미세먼지를 기준으로 정렬한다.


seoulair %>% 

  filter(!is.na(pm10)&!is.na(pm2.5)) %>% 

  group_by(season) %>% 

  summarise(pm10_m = mean(pm10),

            pm2.5_m = mean(pm2.5)) %>% 

  arrange(pm10_m)

결과치를 보면, 우리 모두가 알다시피 겨울,봄 황사의 계절에 가장 높았고, 여름과 가을에 가장 낮은 것을 확인 할 수 있었다. 겨울은 여름에 비해 무려 2배이상이 높다...

코로나 이전에 봄에 항상 황사마스크를 꼈던것도 불편했는데, 이제는 여름에 kf94를 껴도 거뜬하다! 역시 인간은 환경에 적응한다고 했던가...ㅠㅠ


5. 미세먼지 등급 분석

https://cleanair.seoul.go.kr/

데이터를 제공한 서울특별시 대기환경정보에 의하면, 서울시는 미세먼지에 각 등급을 매겨, 관리를 한다고 한다. 아래 등급표를 보고 미세먼지(pm10)의 등급 파생변수를 만들어 등급, 빈도, 백분율을 표시하고, 빈도가 많은 등급부터 출력하도록 한다.


seoulair %>% 

  filter(!is.na(pm10)) %>% 

  mutate(pm_grade = ifelse(pm10<=30, 'good',

                           ifelse(pm10<=81, 'normal',

                                  ifelse(pm10<=150, 'bad','worse')))) %>% 

  group_by(pm_grade) %>% 

  summarise(cnt = n()) %>%  

  mutate(total = sum(cnt), 

         pct = round(cnt/total*100,1)) %>% 

  select(pm_grade, cnt, pct) %>% 

  arrange(desc(cnt))


'좋음'이 50.4%, '보통'dl 44.3%로 대체적으로 서울시의 미세먼지 수준은 좋다고 할 수 있다. 역시 계절에 따른 빈부격차(?)가 극심하다는 것을 알 수 있다. 물론 타의적인 미세먼지 요인이 크지만... 크흠


6. 구별 미세먼지 등급 비교

구별로 미세먼지 등급을 분류한 후, good의 빈도와 비율을 구하고 비율이 높은 순으로 정렬해 상위 5개 구를 추출하도록 한다.


사실 나같은 초보자는 아래와 같은 코드를 짤 확률이 매우매우 높다...!

위에서 부여한 등급 코드에, 구별로 나누라고 했고 나름 완벽하다고 하는데 아래는 에러가 난다.


seoulair %>% 

  filter(!is.na(pm10)) %>% 

  mutate(pm_grade = ifelse(pm10<=30, 'good',

                           ifelse(pm10<=81, 'normal',

                                  ifelse(pm10<=150, 'bad','worse')))) %>% 

  group_by(district) %>% 

  summarise(cnt = n()) %>% 

  mutate(total = sum(cnt), 

         pct = round(cnt/total*100,1)) %>%

  filter(pm_grade == 'good') %>% 

  select(district, cnt, pct) %>% 

  arrange(desc(pct)) %>% 

  head(5)


위의 에러를 보면, 두 번째 filter에 기입한 pm_grade를 찾을 수 없다는 뜻!


해답은 아래와 같다. group_by에 구별로 미세먼지 등급(pm_grade)를 교차 분류하는 것이 이 문제의 핵심이라고 한다. 이렇게 복수의 구별은 정말 초보자에게 너무 생각하기 어렵다... 그저 이것저것 넣고 빼다가 출력되면 다행 ㅠㅠ


seoulair %>% 

  filter(!is.na(pm10)) %>% 

  mutate(pm_grade = ifelse(pm10<=30, 'good',

                           ifelse(pm10<=81, 'normal',

                                  ifelse(pm10<=150, 'bad','worse')))) %>% 

  group_by(district, pm_grade) %>% # district별로 pm_grade 등급 분류

  summarise(cnt = n()) %>%  

  mutate(total = sum(cnt), 

         pct = round(cnt/total*100,1)) %>% 

  filter(pm_grade == 'good') %>% 

  select(district, cnt, pct) %>% 

  arrange(desc(pct)) %>% 

  head(5)


앞서 서대문구의 미세먼지 평균이 가장 낮은 것과 연관 있듯, 서대문구의 '좋음'이 1년 중 222일 이였으며, 이는 61.8%로 가장 높았다. 서울시의 공기좋은 서대문구로 한번 이쁜 카페 탐방을 가고 싶다 :D


7. 1년간 미세먼지 추이를 그래프로 그리기

x축을 date, y축을 pm10으로 하는 선 그래프를 그린다.


ggplot(data = seoulair, aes(x = date, y = pm10)) + 

  geom_line() 


위와 같이 출력이 되었다. 앞서 미세먼지가 가장 안좋았던 강동구의 5월이 한눈에 보인다.

그런데 아무리 생각해도 x축 date의 값이 광범위해서 그런지 수치가 눈에 보이지 않아 불편해서 이리저리 구글링을 해봤지만 별다른 방법이 없었다.

해서 생각한 것이, 현재 date변수는 chr인데 이걸 date 유형으로 바꿔서 출력해보면 어떨까? 하는 호기심이 생겨 도전!


seoulair$date <- as.Date(seoulair$date)

ggplot(data = seoulair, aes(x = date, y = pm10)) + 

  geom_line() 


결과는 아래와 같다.

이쁘게 x축에 날짜가 분기별로 나타나는 것을 확인 할 수 있다! 

초보자로서는 이런 생각을 해냈다는게 장족의 발전... 장하다! ㅎㅎ

그런데 극단치 때문에 다른 데이터들의 차이가 좀 없어져서 구분히 힘들긴 하다, 시간 여유되시는 분들은 따로 극단치 제거 후 확인하시면 될 것 같다...  :3

2019년 통계에는 여름과 가을에 미세먼지량이 적고, 봄과 겨울에 많은 사실이 더 잘 나타났다.


8. 계절별 미세먼지 등급 비율 그래프 그리기


앞서 구한 미세먼지 등급 pm_grade 변수와 계절별로 나눈 season 변수를 활용하여 '계절별 미세먼지 등급비율'을 막대그래프로 그린다.

아직 이런 문장이 안 익숙한 나는 그저 교재를 보고 따라했다..


1) 그래프 그리기 전, season_grade 객체 생성

아래 코드와 같이, 앞서 구했던 구별 미세먼지 등급 처럼, 교차분류를 통해 계절별 미세먼지 등급을 구해 막대그래프에 할당할 객체를 만들어 준다.

계속해서 summarise와 mutate를 이용한 비율로 나타내는 것은 수치가 너무 커서 그렇다.


season_grade <- seoulair %>% 

  filter(!is.na(pm10)) %>% 

  mutate(pm_grade = ifelse(pm10<=30, 'good',

                           ifelse(pm10<=81, 'normal',

                                  ifelse(pm10<=150, 'bad','worse')))) %>% 

  group_by(season, pm_grade) %>% 

  summarise(cnt = n()) %>%  

  mutate(total = sum(cnt), 

         pct = round(cnt/total*100,1)) %>% 

  select(season, pm_grade, pct)


만들어진 season_grade의 값은 위와 같다. 계절별 미세먼지 등급을 분류한 후, 비율로 나타냈다. 교차 분류는 아직까지 좀 혼동하기 쉽다 나에겐...


2) 막대 그래프 그리기

계절별로 등급 비율을 구분해서 표시해야하기 때문에 fill에 pm_grade를 할당해준다.

복수의 세로 막대 그래프로 그리기 위해서 geom_col은 dodge로, x축의 season의 순서는  scale_x_discrete 함수를 이용해 지정해줬다.

범례의 등급 순서가 마음에 들지 않아 바꾸고 싶어 이리저리 알아봤지만, 범례 순서만 바뀌고 막대 그래프에는 영향이 안가 따로 움직이는 방법이 대부분이길래 그냥 냅뒀다. 나중에는 처리할 수 있는 실력으로 성장하길 희망한다.


ggplot(data = season_grade, aes(x = season, y = pct, fill = pm_grade)) + 

  geom_col(position = 'dodge') +

  scale_x_discrete(limits = c('spring', 'summer','autume','winter')) +

  ggtitle('2020~2021년 서울의 계절별 미세먼지 실태') +

  xlab('계절') +

  ylab('등급별 비율')


참고로,   geom_col()로 그래프를 그리면 아래와 같이 누적 막대 그래프로 나온다. 한눈에 보기에는 이것도 나쁘지 않다!



9. 25개 구의 연간 초미세먼지 평균을 구하고, 평균의 내림차순으로 막대그래프를 그려라. 막대 그래프는 90도로 회전시켜서 y축에 구 이름을 표시한다.


난 이 문제가 좀 어려웠다.

내가 처음 짠 코드는 아래와 같다.


1) 그래프 그리기 전, district_2 객체 생성


district_pm2.5<- seoulair %>% 

  filter(!is.na(pm2.5)) %>% 

  group_by(district) %>% 

  summarise(m = mean(pm2.5))


2) 막대 그래프 그리기

여기서 엄청 헤맸는데, 그 헤맴의 과정은 아래와 같다.


# 첫 번째 시도

: 우선 아래 까지만 출력해보면, 뭔가 그럴싸 하게 그림이 나온다.

ggplot(data = district_pm2.5, aes(x = reorder(m, district), y = district)) +

  geom_bar(stat='identity') 


그런데 문제대로라면, coord_flip()함수를 써서 이게 y축에 구이름이 있는 가로로 누운 그래프가 나와야 한다. 그래서 함수를 추가해 보니 엉뚱하게도 아래와 같은 그래프가 나온다...!


ggplot(data = district_pm2.5, aes(x = reorder(m, district), y = district)) +

  geom_bar(stat='identity') +

  coord_flip()

마치 세상이 뒤집어진 것 마냥, 그래프는 물론 x,y도 바껴서 여기서 멘붕이 왔다.


# 두 번째 시도

:geom_bar가 아닌가 싶어, geom_col를 넣었는데도 그래프는 여전히 세상이 뒤집힌 상태였다...ㅜ

ggplot(data = district_pm2.5, aes(x = reorder(m,district), y = district)) +

  geom_col() +

  coord_flip() 


# 세 번째 시도

: 그냥 해답을 보고 따라했다. 즉, 애초에 답으로 출력되어야 하는 형태가 y축 district + 누운 그래프 이기 때문에, 초장에 x축 district + 세로 그래프를 만들고 뒤집기!를 했었어야 했다는 것...


ggplot(data = district_pm2.5, aes(x = reorder(district,m), y = m)) +

  geom_bar(stat='identity') 


ggplot(data = district_pm2.5, aes(x = reorder(district,m), y = m)) +

  geom_bar(stat='identity') +

  coord_flip()

두둥,, 이렇게 이쁜 그림이 나온다.


10. 연간 초미세먼지 평균이 가장 높은 구와 가장 낮은 구의 평균이 통계적으로 차이가 있는지 검정하라.


t.test를 사용하는 문제이다. 그러니까 먼저 가설을 세우기 위해 단순한 평균 차이가 나는 두 항목을 구해 추출한 후 t.test를 하면 된다. 


seoulair %>% 

  filter(!is.na(pm2.5)) %>% 

  group_by(district) %>% 

  summarise(m = mean(pm2.5)) %>% # step1

  filter(m == max(m)|

         m == min(m)) # step2

위와 같이 가장 높은 미세먼지 평균을 가진 곳은 금천구이고, 반대는 성북구 였다. 이 두항목을 추출한 테이블을 따로 제작하여 테스트 하도록 한다.


아래와 같이 두 장소만 해당하는 값은 아래와 같다.

district_2 <- seoulair %>% 

  filter(district %in% c('금천구','성북구'))


이때, t.test는 두 집단의 차이 검정 시 쓰이며, 


" t.test(data = 데이터세트, 종속변수(비교값)~독립변수(비교대상)) "


비교하고자하는 pm2.5가 앞에 나오고, 명목척도인 district는 뒤에 넣는다. 

그러니까 두 지역에 따라 미세먼지 평균이 통계적으로 유의미한 차이가 있는지 알아 보는 것!


귀무가설(H0): 두 지역의 미세먼지 평균은 차이가 없다.

대립가설(H1): 두 지역의 미세먼지 평균은 차이가 있다.


t.test(data = district_2, pm2.5~district)


결과는 아래와 같다.


"유의수준이 p<.001이기 때문에 금천구와 성북구의 평균 차이는 통계적으로 유의미 하다"


상관관계 분석

cor.test(seoulair$pm10, seoulair$pm2.5)


귀무가설(H0): 미세먼지와 초미세먼지는 상관관계가 없다.

대립가설(H1): 미세먼지와 초미세먼지는 상관관계가 있다.


유의수준은 p<.001이어서 통계적으로 의미가 있고, 상관계수(r)은 0.71이다. 미세먼지와 초미세먼지는 정적인 상관관계가 있다고 할 수 있다.


단순회귀분석

lm(data = seoulair, pm2.5~pm10) 

귀무가설(H0): 미세먼지는 초미세먼지에 영향을 주지 않는다.

대립가설(H1): 미세먼지는 초미세먼지에 영향을 준다.

결과를 보면, pm10의 계수(Coefficients)는 0.2815이며, 절편은 9.7157이다. 

따라서 pm2.5 = 9.7157 + 0.2815*pm10 이라는 식을 세울 수 있다. 즉, 미세먼지가 1단위 올라갈 때마다 초미세먼지는 0.2815씩 증가한다(독립변수가 종속변수에 미치는 영향력).


RA <- lm(data = seoulair, pm2.5~pm10) 

summary(RA)


"미세먼지가 초미세먼지에 미치는 영향력(β)은 유의수준p>.001에서 0.281541이다. 회귀모형은 적합하여(p<.001) 회귀계수는 통계적으로 유의미하다. 따라서 미세먼지는 초미세먼지에 정적인 영향을 준다. 또한, 수정된 결정계수는 0.5181로 회귀식의 설명력에 신뢰성을 가지고 있다고 할 수 있다. "


lm(data = seoulair, pm10~pm2.5)

귀무가설(H0): 초미세먼지는 미세먼지에 영향을 주지 않는다.

대립가설(H1): 초미세먼지는 미세먼지에 영향을 준다.

동일하게, pm2.5의 계수(Coefficients)는 1.841이며, 절편은 0.772이다. 

따라서 pm10 = 0.772 + 1.841*pm2.5 이라는 식을 세울 수 있다. 즉, 초미세먼지가 1단위 올라갈 때마다 미세먼지는 1.841씩 증가한다(독립변수가 종속변수에 미치는 영향력).


RA2 <- lm(data = seoulair, pm10~pm2.5) 

summary(RA2)


"초미세먼지가 미세먼지에 미치는 영향력(β)은 유의수준p>.001에서 1.84057이다. 회귀모형은 적합하여(p<.001) 회귀계수는 통계적으로 유의미하다. 따라서 초미세먼지는 미세먼지에 정적인 영향을 준다. 또한, 수정된 결정계수는 0.5181로 회귀식의 설명력에 신뢰성을 가지고 있다고 할 수 있다. "

작가의 이전글 R 프로그래밍 - 2021 서울대기오염 분석 3
브런치는 최신 브라우저에 최적화 되어있습니다. IE chrome safari