https://www.youtube.com/watch?v=AWtFFSoaEAM&list=PLY0OaF78qqGAxKX91WuRigHpwBU0C2SB_&index=14

# 공분산분석

# 공분산분석의 특징
- 분산분석은 집단간의 차이를 검정을 함
- 예)일원분산분석을 이용해 세 종류의 비만치료제 간의 효과가 차이가 있는지 검증을 할 수가 있음
- 그런데 비만은 스트레스 와 밀접한 관련이 있다고 알려져 있음
- 그리고 3가지 치료제 집단의 포함된 참여자들의 스트레스 수준은 서로 다를 수 있음
- 따라서 "일원분산분석을 통해 비만치료제 간 치료효과의 차이가 있다"고 결론을 내렸으면 그 이유는 치료제 때문만이 아니라 어쩌면 참여자들의 스트레스 수준이 달라서 그랬을 수도 있음
- 참여자의 스트레스 수준을 사전에 측정해서 그 데이터를 갖고 있다면 참여자의 스트레스 수준을 통제해 치료제로 인한 순수한 차이를 검증 가능. 스트레스의 영향을 배제한 상태에서, 즉, 스트레스 수준이 일정하다는 가정하에서 치료제 간의 차이가 존재하는지 검증
- 이때 이 스트레스 변수를 공변량이라고 하고 이렇게 공변량을 모델에 포함시켜서 분산분석모델을 확장한 것을 공분산분석이라고 함
- 공변량은 일반적으로 연속형 변수를 가정 함
# 공분산분석 수행절차 (R 사용)
https://youtu.be/AWtFFSoaEAM?list=PLY0OaF78qqGAxKX91WuRigHpwBU0C2SB_&t=100
# 13.통계데이터분석 - 분산분석 - 공분산분석 🔑 ANCOVA | analysis of covariance | 공변량(covariate) | 독립변수와 종속변수 간의 순수한 영향관계
# https://www.youtube.com/watch?v=AWtFFSoaEAM&list=PLY0OaF78qqGAxKX91WuRigHpwBU0C2SB_&index=14
### 공분산분석 수행
# https://youtu.be/AWtFFSoaEAM?list=PLY0OaF78qqGAxKX91WuRigHpwBU0C2SB_&t=100
## faraway 데이터셋의 구조
install.packages("faraway")
library(faraway)
str(sexab)
# Result
# 'data.frame': 76 obs. of 3 variables:
# $ cpa : num 2.048 0.839 -0.241 -1.115 2.015 ...
# $ ptsd: num 9.71 6.17 15.16 11.31 9.95 ...
# $ csa : Factor w/ 2 levels "Abused","NotAbused": 1 1 1 1 1 1 1 1 1 1 ...
?sexab
# sexab 데이터셋은 아동기 성폭력 경험이 성인의 정신건강에 미치는 영향을 분석한 연구에서 사용됨
# 76개의 관측 값과 3개의 변수(cpa, ptsd, csa)를 포함
# csa 변수: 독립변수인 집단변수로 아동기 성폭력 경험을 나타냄. 아동기 성폭력 경험 여부에 따라 성인여성을 2개의 집단으로 구분. 성폭력을 경험했으면 abused, 경험하지 않았으면 not abused로 코딩
# ptsd 변수: 종속변수로써 외상후스트레스장애를 나타냄
# cpa 변수: 공변량으로 사용되는 변수로써 아동기 신체적 학대를 의미
# 여기서는 csa 집단변수로 해서 아동기 성폭력 경험 유무가 성인의 외상후스트레스장애에 미치는 영향을 분석하는데 공변량으로서 아동기 신체적 학대를 포함시켜 아동기 성폭력 경험 유무가 외상후스트레스장애에 미치는 순수한 영향관계를 파악하고자 함
## tapply()함수-집단별 요약 통계량 계산
# 여기서는 아동기 성폭력 경험 유무에 따라서 두 개 집단으로 구분되고 그 집단별로 표본크기, 평균, 표준편차 계산
# 1) 집단별 표본크기
tapply(sexab$ptsd, sexab$csa, length)
# 1st 인수: 대상벡터 지정
# 2nd 인수: 집단변수 지정
# 3rd 인수: 집단별로 적용할 함수 지정. length() 함수 지정하면 집단별로 관측 값 개수를 구할 수 있음
# Result
# Abused NotAbused
# 45 31
# 2) 집단별 평균
tapply(sexab$ptsd, sexab$csa, mean)
# 아동기 성폭력 경험유무에 따라서 외상후스트레스장애의 평균값 계산
# Result
# Abused NotAbused
# 11.941093 4.695874
# 3) 표준편차
tapply(sexab$ptsd, sexab$csa, sd)
# Result
# Abused NotAbused
# 3.440152 3.519743
# 결과해석1) 총 76명의 참여자 중 아동기 성폭력 경험자는 45명, 그렇지 않은 자는 30명
# 결과해석2) 평균값에서 외상후스트레스장애 정도는 아동기 성폭력 경험 유무에 따라서 큰 차이 보임
# 아동기 성폭력 경험한 경우가 11.9, 그렇지 않은 경우가 4.7
# 그래서 아동기 성폭력 경험 여성이 그렇지 않은 여성에 비해 2배 이상 높은 외상후스트레스장애를 겪었음
# 결과해석3) 반면 표준편차는 두 집단간의 큰 차이가 없음
# 이 표본 데이터를 바탕으로 아동기성폭력 경험과 외상후스트레스장애 간 관계를 통계적으로 검정
## (공분산분석 전) 기본적인 일원분산분석 수행
# 먼저 기본적인 일원분산분석 수행
sexab.aov <- aov(ptsd ~ csa, data=sexab)
# 1st 인수: 포뮬러 형식
# 종속변수: 외상후스트레스장애(ptsd)
# 독립변수: 아동기 성폭력 경험 유무(csa)
# data 인수: 사용하는 데이터셋 sexab 지정
sexab.aov
summary(sexab.aov)
# Result
# Df Sum Sq Mean Sq F value Pr(>F)
# csa 1 963.5 963.5 79.9 2.17e-13 ***
# Residuals 74 892.4 12.1
# ---
# Signif. codes:
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 결과해석1) 분산분석표에서 F값이 79.9, 그에 대응되는 유의확률이 2.17 곱하기 10의 -13임. 유의확률이 0.05나 0.01에 비해 매우 작기 때문에 귀무가설 기각하고 대립가설 채택
# 귀무가설은 "아동기 성폭력 경험에 따라 외상후스트레스장애는 평균값의 차이가 없다"이기 때문에
# 이 귀무가설 기각한다는 것은 "아동기 성폭력 경험 유무에 따라 외상후스트레스장애 정도는 차이가 있다" 라고 결론 내릴 수 있음
# 독립변수인 아동기 성폭력 경험은 중속변수인 외상후스트레스장애에 영향을 미침
## (일원분산분석 수행 후) 공분산분석 수행
# 그런데 어떤 사람이 얘기를 하기를 "외상후스트레스장애는 아동기 성폭력 경험 뿐만 아니라 아동기 신체적 학대에서도 영향 받을 것이라고 얘기함
# 아동기 성폭력 경험이 외상후스트레스장애 미치는 순수한 영향을 파악하기 위해 아동기 신체적 학대를 고려해야 된다고 주장
# 아동기 신체적 학대를 추가로 고려하기 위해 공분상분석을 수행
# 아동의 신체적 학대를 공변량으로 모델에 투입해 통제하여 아동기 성폭력 경험과 외상후스트레스장애 간의 순수한 영향 관계 분석가능
sexab.aov <- aov(ptsd ~ cpa + csa, data=sexab)
summary(sexab.aov)
# 틸드 앞에는 검증하고자 하는 종속변수 외상후스트레스장애 지정
# 틸드 뒤에는 공변량 추가로 더 지정. 그리고 + 로 연결해 집단 변수인 아동기 성폭력 경험 유무를 나타내는 변수 지정
# 그러면 이 모델에서는 기존의 독립변수인 아동기 성폭력 경험과 함께 새롭게 공변량이 추가돼, 아동기 신체적 학대가 일정하다는 가정 하에서, 아동기 성폭력 경험과 외상후스트레스장애 간의 순수한 관계를 파악을 해 볼 수 있음
# Result
# Df Sum Sq Mean Sq F value Pr(>F)
# cpa 1 449.8 449.8 41.98 9.46e-09 ***
# csa 1 624.0 624.0 58.25 6.91e-11 ***
# Residuals 73 782.1 10.7
# ---
# Signif. codes:
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 결과해석1) 분산분석표를 보면 cpa와 csa 모두 p값이 유의수준 0.05에 비해 매우 작음. 따라서
# "아동기 신체적 학대는 외상후스트레스장애와 관련이 있고",
# 또 아동기 신체적 학대를 통제한 상태에서 "아동기 성폭력 경험은 외상후스트레스장해에 영향을 미침"
# 즉, 외상후스트레스장애에 영향을 미치는 신체적 학대를 통제한 후에도, 아동기 성폭력을 경험한 집단과 그렇지 않은 집단은 외상후 스트레 스장애에 있어서 통계적으로 유의한 차이를 나타냄
## effect() 함수로 계산-공변량 영향 제거 후, 조정된 집단별 종속변수(외상후스트레스장애의) 평균값 계산
# 공변량인 "아동기 신체적 학대" 영향 제거 후 "조정된 외상후스트레스장애의 집단별 평균이 어떻게 달라지는지" 확인가능
install.packages("effects")
library(effects)
effect("csa", sexab.aov)
# 1st 인수: 집단변수 지정
# 2nd 인수: 공분산 분석 결과 생성된 모델 객체 지정
# Result
# csa effect
# csa
# Abused NotAbused
# 11.544429 5.271677
# 결과해석1) 앞서 tapply() 함수로 구한 조정 전 집단평균에 비해 성폭력 경험 집단은 약간 낮아졌고 반대로 성폭력 비경험 집단은 약간 높아 짐
# 그래서 공변량 아동기 신체적 학대 영향을 제거하면 아동기 성폭력 경험집단과 그렇지 않은 집단 간의 외상후스트레스장애의 평균값 차이가 종전의 비해 작아짐
## ancova()-공분산 분석 결과를 그래프로 표현
# HH패키지의 ancove() 함수로 공분산분석 결과를 그래프로 표현
# aov() 함수로 공분산분석 수행 시 지정했던 것과 동일한 방식으로 종속변수-공변량-독립변수 간의 관계를 포뮬러 형식으로 인수 지정
library(HH)
windows(width=12, height=8)
ancova(ptsd ~ cpa + csa, data=sexab)
# Result
# Analysis of Variance Table
#
# Response: ptsd
# Df Sum Sq Mean Sq F value Pr(>F)
# cpa 1 449.80 449.80 41.984 9.462e-09 ***
# csa 1 624.03 624.03 58.247 6.907e-11 ***
# Residuals 73 782.08 10.71
# ---
# Signif. codes:
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 그래프: https://i.ibb.co/jZNTGSd/13-ancova.png
# 결과해석1) ancova() 함수 공분산분석에서 봤던 것과 동일한 분산분석표와 그래프 함께 출력
# 도표해석1)
# x축: 공변량으로 사용하고 있는 "아동기 신체적 학대"(cpa) 배치
# y축: 종속변수인 "외상후스트레스장애" 배치
# 패널: 독립변수 "아동기 성폭력 경험 유무"는 두 개 패널로 표현. "아동기 성폭력 경험 유무"에 따라 각각 별도 패널의 집단 별로 그래프 생성
# 도표해석2)
# 왼쪽은 아동기 성폭력 경험 있는 집단, 오른쪽은 없는 집단임
# 아동기 신체적 학대로부터 외상후스트레스장애를 예측하는 이 두 개의 회귀선은 두 집단 모두에서 동일 기울기, 다른 절편으로 표현
# 두 회기선 기울기가 같은 이유는 아동기 신체적 학대 수준이 외상후스트레스장애에 미치는 영향이 두 집단에서 일정하도록 공변량으로 통제 했기 때문임
# 두 그래프에서 아동기 신체적 학대 경험이 증가할수록 외상후스트레스장애 역시 증가함
# 아동기 성폭력 경험 집단이 그렇지 않은 집단보다 더 큰 외상후 스트레스 장애를 겪는 것을 확인됨