내용정리

• 1. 상관관계 : 숫자형-숫자형 변수의 관계, 직선성에 대한 척도, 그래프시 산점도 / 분석시 상관관계
• 2. t-test : 범주(요인) 1개 속 2개의 집단(수준, class, 범주의 종류)별 숫자형의 평균차이 검정
• 3. anova : 범주(요인) 1개 or 2개 속 3개의 집단(수준, class, 범주의 종류)별 숫자형의 평균차이를 분산분석으로 검정

# 0. 데이터 준비

library(moonBook)

data(acs)

head(acs)

##   age    sex cardiogenicShock   entry              Dx   EF height weight

## 1  62   Male               No Femoral           STEMI 18.0    168     72

## 2  78 Female               No Femoral           STEMI 18.4    148     48

## 3  76 Female              Yes Femoral           STEMI 20.0     NA     NA

## 4  89 Female               No Femoral           STEMI 21.8    165     50

## 5  56   Male               No  Radial          NSTEMI 21.8    162     64

## 6  73 Female               No  Radial Unstable Angina 22.0    153     59

##        BMI obesity  TC LDLC HDLC  TG  DM HBP smoking

## 1 25.51020     Yes 215  154   35 155 Yes  No  Smoker

## 2 21.91381      No  NA   NA   NA 166  No Yes   Never

## 3       NA      No  NA   NA   NA  NA  No Yes   Never

## 4 18.36547      No 121   73   20  89  No  No   Never

## 5 24.38653      No 195  151   36  63 Yes Yes  Smoker

## 6 25.20398     Yes 184  112   38 137 Yes Yes   Never

# 1. 숫자형 필요칼럼만 뽑기

acs2 <- acs[, c("height", "weight", "BMI")]



# 2. 상관계수 행렬

cor( acs2, use = "na.or.complete")

##              height    weight          BMI

## height  1.000000000 0.6315767 -0.009049596

## weight  0.631576661 1.0000000  0.762441135

## BMI    -0.009049596 0.7624411  1.000000000

# 3. 산점도+상관계수 행렬 pairplot 1

library(psych)

pairs.panels(acs2)

# 4. 산점도+상관계수 행렬 pairplot 2

library(PerformanceAnalytics)

## Loading required package: xts

## Loading required package: zoo

## 

## Attaching package: 'zoo'

## The following objects are masked from 'package:base':

## 

##     as.Date, as.Date.numeric

## 

## Attaching package: 'PerformanceAnalytics'

## The following object is masked from 'package:graphics':

## 

##     legend

chart.Correlation(acs2, histogram = TRUE, pch = 19)

# 5. 상관계수 heatmap

library(corrplot)

## corrplot 0.84 loaded

corrplot(cor(acs2,use="na.or.complete"))

#### 1. 독립표본 t-test

# 데이터 준비

# 1. 독립성 생략

# 2. 정규성 생략

# 3. 등분산성 검정 :  집단별로 BMI를 인덱싱한 뒤, 등분산성 체크

var.test(acs[ acs$sex == "Male", ]$BMI, acs[acs$sex == "Female",]$BMI)

## 

##  F test to compare two variances

## 

## data:  acs[acs$sex == "Male", ]$BMI and acs[acs$sex == "Female", ]$BMI

## F = 0.82798, num df = 508, denom df = 254, p-value = 0.07756

## alternative hypothesis: true ratio of variances is not equal to 1

## 95 percent confidence interval:

##  0.6663069 1.0210005

## sample estimates:

## ratio of variances 

##          0.8279795

# 4. 독립표본 t-test

t.test(BMI ~ sex, data = acs, var.equal = T)

## 

##  Two Sample t-test

## 

## data:  BMI by sex

## t = -0.50823, df = 762, p-value = 0.6114

## alternative hypothesis: true difference in means is not equal to 0

## 95 percent confidence interval:

##  -0.6348532  0.3737344

## sample estimates:

## mean in group Female   mean in group Male 

##             24.19492             24.32548

#### 2. 보정된 t-test 

summary(aov(BMI ~ age + sex, data = acs)) # 기존 sex가 +뒤에 있으며, 이 것의 p-value 보는 것이 목적

##              Df Sum Sq Mean Sq F value   Pr(>F)    

## age           1    386   385.8  36.097 2.91e-09 ***

## sex           1     26    25.9   2.427     0.12    

## Residuals   761   8134    10.7                     

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 93 observations deleted due to missingness

#### 3. 대응표본 t-test

# 데이터 준비 : 같은 길이를 가진 데이터야한다.

Before_score = c(80,30,50,40,40,10,20,30,60,70)

After_score = c(95,40,60,45,40,15,20,25,67,90)



t.test(Before_score, After_score, paired=T)

## 

##  Paired t-test

## 

## data:  Before_score and After_score

## t = -2.8423, df = 9, p-value = 0.01933

## alternative hypothesis: true difference in means is not equal to 0

## 95 percent confidence interval:

##  -12.032489  -1.367511

## sample estimates:

## mean of the differences 

##                    -6.7

#### 1. 1way anova

# 1.숫자형으로 되어있다면 범주형으로 타입 변경 생략

# 2. 범주별 숫자의 boxplot, mean, sd 생략

# 3. 등분산성 검정 생략

# 4. anova

aov <- aov(TG ~ smoking, data= acs)

summary(aov)

##              Df  Sum Sq Mean Sq F value Pr(>F)  

## smoking       2   65695   32848   4.008 0.0185 *

## Residuals   839 6876082    8196                 

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 15 observations deleted due to missingness

# 5. 사후분석 Duncan

library(laercio)

LDuncan(aov, "group")

## 

##  DUNCAN TEST TO COMPARE MEANS 

##  

##  Confidence Level:  0.95 

##  Dependent Variable:  TG

##  Variation Coefficient:  72.28546 % 

##  

## 

##  Independent Variable:  smoking 

##   Factors   Means              

##   Smoker    136.490566037736 a 

##   Never     119.496913580247  b

##   Ex-smoker 116.65            b

#### 2. 2way anova

# 1. 2way anova

aov2 <- aov( TG ~ smoking + Dx, data = acs)

summary(aov2) # 두 요인 모두 유의

##              Df  Sum Sq Mean Sq F value   Pr(>F)    

## smoking       2   65695   32848   4.119   0.0166 *  

## Dx            2  201756  100878  12.651 3.87e-06 ***

## Residuals   837 6674326    7974                     

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 15 observations deleted due to missingness

aov2 <- aov( TG ~ Dx + smoking , data = acs)

summary(aov2) # 두 요인 모두 유의 but 수준차이가 난다...

##              Df  Sum Sq Mean Sq F value   Pr(>F)    

## Dx            2  166365   83182  10.432 3.35e-05 ***

## smoking       2  101087   50543   6.338  0.00185 ** 

## Residuals   837 6674326    7974                     

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 15 observations deleted due to missingness

# 2. interactioplot

par(mfrow=c(1,2))

mean.rm.na <- function(x) mean(x,na.rm=T) # interaction.plot에 들어가는 집계함수를 직접 정의***

interaction.plot(acs$smoking,acs$Dx, acs$TG,fun=mean.rm.na)

interaction.plot(acs$Dx,acs$smoking, acs$TG,fun=mean.rm.na)

칼럼명 확인과 변경

# 데이터 준비
iris2 <- iris
head(iris2)

# 칼럼명확인
names(iris2)

## [1] "Sepal.Length" "Sepal.Width"  "Petal.Length" "Petal.Width" 
## [5] "Species"

# 칼럼명 변경
names(iris2) <- c( "first", "second", "third", "fourth", "fifth")
names(iris2)

## [1] "first"  "second" "third"  "fourth" "fifth"

사용자 정의함수로 ( 반복적인 작업인 ) scaling 해보기

• 정규분포에서 random한 숫자 10를 뽑아 4개의 칼럼 만들기
• 그 값들을 scaling을 통해 0~1사이 값으로 만들기
• x - min(x) / max(x) - min(x) => 분모가 가장 크므로 항상 1보다 작다.

# tibble을 활용한 데이터프레임 만들기
# library(tibble)
# 라이브러리 적용과 동시에 해당 함수 사용하기 package :: function
df <- tibble::tibble(
  a = rnorm(10),
  b = rnorm(10),
  c = rnorm(10),
  d = rnorm(10)
)
head(df)

summary(df)

##        a                 b                  c                  d          
##  Min.   :-1.5702   Min.   :-1.55135   Min.   :-1.97816   Min.   :-1.0510  
##  1st Qu.:-0.6281   1st Qu.:-0.52394   1st Qu.:-1.05749   1st Qu.:-0.2120  
##  Median : 0.4213   Median : 0.21181   Median :-0.23076   Median : 0.3042  
##  Mean   : 0.1132   Mean   : 0.06085   Mean   :-0.03817   Mean   : 0.3449  
##  3rd Qu.: 0.7566   3rd Qu.: 0.83831   3rd Qu.: 1.02052   3rd Qu.: 0.7530  
##  Max.   : 1.8824   Max.   : 1.32477   Max.   : 2.16831   Max.   : 1.7938

# scaling은 각 칼럼마다 해줘야하므로 반복적이다.
(df$a - min(df$a, na.rm = T))  /  (max(df$a, na.rm = T) -  min(df$a, na.rm = T))

##  [1] 0.447926340 0.002240828 0.680114530 0.583405776 0.214523807
##  [6] 1.000000000 0.570234402 0.000000000 0.721868086 0.655403845

df$a <- (df$a - min(df$a, na.rm = T))  /  (max(df$a, na.rm = T) -  min(df$a, na.rm = T))
df$b <- (df$b - min(df$b, na.rm = T))  /  (max(df$b, na.rm = T) -  min(df$b, na.rm = T))
df$c <- (df$c - min(df$c, na.rm = T))  /  (max(df$c, na.rm = T) -  min(df$c, na.rm = T))
df$d <- (df$d - min(df$d, na.rm = T))  /  (max(df$d, na.rm = T) -  min(df$d, na.rm = T))

head(df)

# range()함수는 [1]에는 최소값 [2]에는 최대값을 내포한다.
rng <- range( rnorm(10) , na.rm = T)
rng[1] # 최소값

## [1] -0.7202099

rng[2] # 최대값

## [1] 1.028419

# 사용자 정의함수로 scaling 정의하기

rescale <- function (x){
  rng <- range(x, na.rm = T)
  ( x- rng[1] ) / ( rng[2] - rng[1])
  
}

#테스트해보기
rescale(df$a)

##  [1] 0.447926340 0.002240828 0.680114530 0.583405776 0.214523807
##  [6] 1.000000000 0.570234402 0.000000000 0.721868086 0.655403845

rescale( c(0, 5, 10) )

## [1] 0.0 0.5 1.0

# 각 칼럼마다 적용해보기
df$a <- rescale(df$a)
df$b <- rescale(df$b)
df$c <- rescale(df$c)
df$d <- rescale(df$d)

head(df)

one-way ANOVA(요인 = 범주1개) 분석해보기

# 데이터 준비
anova_data<-data.frame(
  group=c(1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3),
  score=c(50.5, 52.1, 51.9, 52.4, 50.6, 51.4, 51.2, 52.2, 51.5, 50.8,47.5, 47.7, 46.6, 47.1, 47.2, 47.8, 45.2, 47.4, 45.0, 47.9,46.0, 47.1, 45.6, 47.1, 47.2, 46.4, 45.9, 47.1, 44.9, 46.2))

# 데이터 확인
head(anova_data)

str(anova_data)

## 'data.frame':    30 obs. of  2 variables:
##  $ group: num  1 1 1 1 1 1 1 1 1 1 ...
##  $ score: num  50.5 52.1 51.9 52.4 50.6 51.4 51.2 52.2 51.5 50.8 ...

# 1. 범주형 데이터를 (숫자로 나뉘어도) 문자형 변수로 변경***
# - as.character ( 칼럼인덱싱 ) 
anova_data$group <- as.character(anova_data$group)
str(anova_data)

## 'data.frame':    30 obs. of  2 variables:
##  $ group: chr  "1" "1" "1" "1" ...
##  $ score: num  50.5 52.1 51.9 52.4 50.6 51.4 51.2 52.2 51.5 50.8 ...

# 2. 집단별 종속변수(숫자형) 데이터의 boxplot 확인 및 tapply(종속변수인덱싱, 범주칼럼인덱싱, 집계함수)로 집단별 종속변수의 평균과 표준편차.
boxplot(score~group,data=anova_data)
tapply(anova_data$score,anova_data$group, mean)

##     1     2     3 
## 51.46 46.94 46.35

tapply(anova_data$score,anova_data$group, sd)

##         1         2         3 
## 0.6834553 1.0415800 0.7763876

# 3. 아노바분석 전, 3집단 등분산성 
# 3_1. bartlett.test ***
bartlett.test( score ~ as.factor(group), data= anova_data) #  0.4368

## 
##  Bartlett test of homogeneity of variances
## 
## data:  score by as.factor(group)
## Bartlett's K-squared = 1.6565, df = 2, p-value = 0.4368

# 3_2. userfriendlyscience패키지의 oneway()함수로 등분산성 테스트 및 바로 oneway 아노바 적용
# install.packages("userfriendlyscience")
library(userfriendlyscience)

#      등분산인지 확인
oneway(factor(anova_data$group),y=anova_data$score,fullDescribe=T,correction=T)

## ### Oneway Anova for y=score and x=group (groups: 1, 2, 3)
## 
## Omega squared: 95% CI = [.78; .92], point estimate = .88
## Eta Squared: 95% CI = [.8; .92], point estimate = .89
## 
##                                    SS Df    MS      F     p
## Between groups (error + effect) 156.3  2 78.15 108.81 <.001
## Within groups (error only)      19.39 27  0.72             
## 
## 
## ### Welch correction for nonhomogeneous variances:
## 
## F[2, 17.55] = 136.29, p < .001.
## 
## ### Brown-Forsythe correction for nonhomogeneous variances:
## 
## F[2, 23.76] = 108.81, p < .001.

#      one-way anova 시행
oneway.test(score~group,data=anova_data,var.equal = T)

## 
##  One-way analysis of means
## 
## data:  score and group
## F = 108.81, num df = 2, denom df = 27, p-value = 1.199e-13

# 4. oneway 아노바 분석 aov( 종속변수 ~ 요인, data= ) : 그룹간 변동, 그룹내변동(Residual) 만 계산
# - H0 : Ma = Mb = Mc ( 모두 평균이 동일 )
# - Sum of Squares = 변동 = 편차(X-m)의 제곱 -> 제곱을 해야 합이 0이 안되므로.
# - group = 요인1에 대한 그룹간 변동 = ( 전체평균 - 그룹의 평균의  ) 의 제곱
# - Residuals = error = 그룹내 변동 =( 각 그룹별 x - 각 그룹별 평균) 의 제곱
# - aov()결과는 p-value가 안나온다.
aov( score ~ group, data = anova_data) # Sum of Square = 변동, Deg. of Freedom -> 변동의 평균 구해야하므로 필요.

## Call:
##    aov(formula = score ~ group, data = anova_data)
## 
## Terms:
##                   group Residuals
## Sum of Squares  156.302    19.393
## Deg. of Freedom       2        27
## 
## Residual standard error: 0.8475018
## Estimated effects may be unbalanced

# 5. summaray(  aov())를 통한  ANOVA 테이블 도출 : 여기서 p-value가 나온다.
# - Df : 자유도 -> F분포를 결정 + 변동들의 평균구할때 분모로 나눠짐
# - Sum Sq : 변동들
# - Mean sq : 변동들의 평균
# - F value : F검정 통계량 = effect(그룹간 변동의 평균) / error(그룹내 변동의 평균)
summary( aov( score ~ group, data = anova_data)  ) # 1.2e-13

##             Df Sum Sq Mean Sq F value  Pr(>F)    
## group        2 156.30   78.15   108.8 1.2e-13 ***
## Residuals   27  19.39    0.72                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# 6. 사후분석1 : laercio 패키지의 LDuncan ( 3집단 비교시 어느것이 차이가 있는지 알아야함)
# - aov()의 결과를 넣는다. 
# install.packages("laercio")
library(laercio)

a1 <- aov( score ~ group, data = anova_data)
LDuncan( a1, "group" ) # 1번 그룹만 차이가 있는 그룹으로 나왔다.

## 
##  DUNCAN TEST TO COMPARE MEANS 
##  
##  Confidence Level:  0.95 
##  Dependent Variable:  score
##  Variation Coefficient:  1.75648 % 
##  
## 
##  Independent Variable:  group 
##   Factors Means   
##   1       51.46 a 
##   2       46.94  b
##   3       46.35  b

# 7. 사후분석2 : 튜키 정직유의검정(TukeyHSD) 및 튜키검정에 의한 평균 차이에 대한 95% 신뢰구간 plot
#    마찬가지로, aov()결과가 들어간다.
#    반드시 칼럼이 문자형이어야한다!
#    자동으로 Holm-Bonferroni correction을 통해 조절된 p-value ( p * nC2 )
TukeyHSD(aov(score~as.character(group),data=anova_data)) # 1-2 ,1-3 만 차이.. 1만 다른 집단!

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = score ~ as.character(group), data = anova_data)
## 
## $`as.character(group)`
##      diff       lwr        upr     p adj
## 2-1 -4.52 -5.459735 -3.5802652 0.0000000
## 3-1 -5.11 -6.049735 -4.1702652 0.0000000
## 3-2 -0.59 -1.529735  0.3497348 0.2813795

plot(TukeyHSD(aov(score~as.character(group),data=anova_data)))

COPD데이터로 one-way 연습 (no quiz)

• 3집단을 가진 범주 : SMK_STAT_TYPE
• 종속변수 : HMG

# 데이터 준비
data <- read.csv("week4/4주차_코드/COPD_Lung_Cancer.csv")
str(data)

## 'data.frame':    944 obs. of  20 variables:
##  $ PERSON_ID            : int  10000065 10000183 10000269 10000471 10000788 10001096 10001122 10001260 10001546 10001919 ...
##  $ SEX                  : int  1 1 1 1 2 2 1 1 1 1 ...
##  $ DTH                  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ CTRB_PT_TYPE_CD      : int  8 9 8 1 7 3 10 10 8 8 ...
##  $ BMI                  : num  25.7 22.8 24.5 21.6 25.6 ...
##  $ BP_LWST              : int  90 100 70 100 54 60 100 73 90 70 ...
##  $ BLDS                 : int  112 105 76 70 113 137 89 98 73 233 ...
##  $ TOT_CHOLE            : int  177 177 156 228 199 182 243 284 221 133 ...
##  $ HMG                  : num  12.2 15.7 13.1 14 12.1 12 17 12.1 15 8.7 ...
##  $ FMLY_CANCER_PATIEN_YN: int  1 1 1 NA 1 1 1 1 1 1 ...
##  $ SMK_STAT_TYPE        : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ DRNK_HABIT           : int  1 5 1 1 1 1 1 1 1 1 ...
##  $ EXERCI_FREQ          : int  1 1 1 1 1 1 1 2 1 1 ...
##  $ lungC                : int  0 0 1 0 0 0 0 0 0 0 ...
##  $ copd                 : int  0 1 1 1 0 0 0 0 0 1 ...
##  $ Asthma               : int  0 1 0 0 1 0 1 1 0 1 ...
##  $ DM                   : int  1 1 1 1 1 0 1 1 1 1 ...
##  $ Tub                  : int  0 0 1 0 0 0 0 0 0 0 ...
##  $ before_op_score      : int  44 44 37 50 31 34 31 45 35 43 ...
##  $ after_op_score       : int  52 53 60 44 43 49 41 51 59 42 ...

# 데이터 확인
tapply(data$HMG, data$SMK_STAT_TYPE, mean, na.rm=T)

##        1        2        3 
## 13.00319 13.78158 13.34242

tapply(data$HMG, data$SMK_STAT_TYPE, sd, na.rm=T)

##        1        2        3 
## 1.349230 1.600059 1.393871

# 1. 범주형인데 int/num 타입으로 된 것의 인덱스만 *** 가져와서, for문을 통해 
#    반복작업인 as.character로 변형해서 덮어씌우기*************8
index <- c(2, 3, 4, 10:18)
for (i in index) {
  data[,i] <- as.character( data[,i] )
}
str(data)

## 'data.frame':    944 obs. of  20 variables:
##  $ PERSON_ID            : int  10000065 10000183 10000269 10000471 10000788 10001096 10001122 10001260 10001546 10001919 ...
##  $ SEX                  : chr  "1" "1" "1" "1" ...
##  $ DTH                  : chr  "0" "0" "0" "0" ...
##  $ CTRB_PT_TYPE_CD      : chr  "8" "9" "8" "1" ...
##  $ BMI                  : num  25.7 22.8 24.5 21.6 25.6 ...
##  $ BP_LWST              : int  90 100 70 100 54 60 100 73 90 70 ...
##  $ BLDS                 : int  112 105 76 70 113 137 89 98 73 233 ...
##  $ TOT_CHOLE            : int  177 177 156 228 199 182 243 284 221 133 ...
##  $ HMG                  : num  12.2 15.7 13.1 14 12.1 12 17 12.1 15 8.7 ...
##  $ FMLY_CANCER_PATIEN_YN: chr  "1" "1" "1" NA ...
##  $ SMK_STAT_TYPE        : chr  "1" "1" "1" "1" ...
##  $ DRNK_HABIT           : chr  "1" "5" "1" "1" ...
##  $ EXERCI_FREQ          : chr  "1" "1" "1" "1" ...
##  $ lungC                : chr  "0" "0" "1" "0" ...
##  $ copd                 : chr  "0" "1" "1" "1" ...
##  $ Asthma               : chr  "0" "1" "0" "0" ...
##  $ DM                   : chr  "1" "1" "1" "1" ...
##  $ Tub                  : chr  "0" "0" "1" "0" ...
##  $ before_op_score      : int  44 44 37 50 31 34 31 45 35 43 ...
##  $ after_op_score       : int  52 53 60 44 43 49 41 51 59 42 ...

# 2. 3집단 등분산성 검정
bartlett.test( HMG ~ as.factor(SMK_STAT_TYPE), data = data) # 0.3899

## 
##  Bartlett test of homogeneity of variances
## 
## data:  HMG by as.factor(SMK_STAT_TYPE)
## Bartlett's K-squared = 4.3473, df = 2, p-value = 0.1138

# 3. 3집단 정규성 검정
unique(data$SMK_STAT_TYPE) # "1" "3" "2" NA 

## [1] "1" "3" "2" NA

data1 <- data[ data$SMK_STAT_TYPE == "1", ]
data2 <- data[ data$SMK_STAT_TYPE == "2", ]
data3 <- data[ data$SMK_STAT_TYPE == "3", ]

shapiro.test(data1$HMG)

## 
##  Shapiro-Wilk normality test
## 
## data:  data1$HMG
## W = 0.99681, p-value = 0.1639

shapiro.test(data2$HMG) # 0.002898 정규성 만족안하지만, 한다고 가정

## 
##  Shapiro-Wilk normality test
## 
## data:  data2$HMG
## W = 0.94628, p-value = 0.002898

shapiro.test(data3$HMG)

## 
##  Shapiro-Wilk normality test
## 
## data:  data3$HMG
## W = 0.99057, p-value = 0.7168

# 3. 1way anova
one_way <- aov(HMG ~ SMK_STAT_TYPE, data=data)
summary(one_way) # 3.97e-06 -> 평균차이가 명확히 차이남

##                Df Sum Sq Mean Sq F value   Pr(>F)    
## SMK_STAT_TYPE   2   47.8  23.915   12.61 3.97e-06 ***
## Residuals     892 1691.3   1.896                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 49 observations deleted due to missingness

# 4. 사후검정
library(laercio)
LDuncan(one_way, "group") # a, b, c  3집단 모두 평균이 다른 것으로 나오

## 
##  DUNCAN TEST TO COMPARE MEANS 
##  
##  Confidence Level:  0.95 
##  Dependent Variable:  HMG
##  Variation Coefficient:  10.50583 % 
##  
## 
##  Independent Variable:  SMK_STAT_TYPE 
##   Factors Means               
##   2       13.7815789473684 a  
##   3       13.3424242424242  b 
##   1       13.0031944444444   c

2way ANOVA 와 interaction plot

• 종속변수 : BMI
• 요인 : DM 과 SMK_STAT_TYPE

# 각 요인별 boxplot
par(mfrow=c(1,2))
boxplot(BMI ~ DM,data=data, 
        boxwex = 0.7,
        main="BMI~DM", 
        col = c("yellow","orange", "green"))
boxplot(BMI ~ SMK_STAT_TYPE,data=data, 
        boxwex = 0.5,
        main="BMI~Smkoing", 
        col = c("blue", "coral"))

# 요인1(DM이력-3가지집단), 요인2(흡연이력-3가지집단)에 따른 BMI 평균의 차이
two_way <- aov( BMI ~ DM + SMK_STAT_TYPE ,data = data)
summary(two_way) # DM : 1.66e-05 *** , SMK_STAT_TYPE : 0.136

##                Df Sum Sq Mean Sq F value   Pr(>F)    
## DM              1    199  198.96  18.748 1.66e-05 ***
## SMK_STAT_TYPE   2     42   21.19   1.997    0.136    
## Residuals     896   9509   10.61                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 44 observations deleted due to missingness

# 요인2개를 순서바꿔서
two_way2 <- aov( BMI ~ SMK_STAT_TYPE + DM  ,data = data)
summary(two_way2) # DM :  3.1e-05 ***, SMK_STAT_TYPE : 0.0746

##                Df Sum Sq Mean Sq F value  Pr(>F)    
## SMK_STAT_TYPE   2     55   27.63   2.604  0.0746 .  
## DM              1    186  186.08  17.534 3.1e-05 ***
## Residuals     896   9509   10.61                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 44 observations deleted due to missingness

# interaction plot
# - interaction.plot( x축-요인1  -  범례 - 요인2  -  종속변수)
# 요인1-종속변수에  요인2가 영향을 끼치는지 보기
par(mfrow=c(1,2))
interaction.plot(data$SMK_STAT_TYPE, data$DM, data$BMI)
interaction.plot(data$DM, data$SMK_STAT_TYPE, data$BMI)

등분산성, 정규성 검정 후 t-test 종류별로 해보기

# 데이터 준비
t_data <- data.frame(
  group = c(1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2),
  score = c(100,100,80,50,40,90,20,50,50,70,30,40,30,70,30,40,30,60,30,60),
  age = c(70,70,70,60,50,60,60,60,50,50,20,30,20,20,25,10,13,10,10,10)
)

head(t_data)

# 2집단(2범주)별 따로 데이터 인덱싱 -> 그래야 등분산성 var.test 및 t.test :  
# - var.test, t.test ( 집단1의 데이터, 집단2의 데이터)
# 범주별 특정범주만 가져오기 -> 행인덱싱 자리에 넣기
t_data1 <- t_data[ t_data$group == 1, ]
t_data2 <- t_data[ t_data$group == 2, ]


# t-test 전에 확인사항 3가지 중 2가지 - 1. 등분산성 2. 정규성 => 둘다 p-value 0.05보다 커야 만족

#### 1. 등분산성 확인후 독립표본 T-검정
####    H0 : 두 집단간 평균의 차이가 없다.
# 1_1. 등분산성 by F검정( = effect(그룹간 변동의 평균) / error(그룹내 변동의 평균) )
var.test( t_data1$score , t_data2$score ) # 0.1093 -> 2그룹간 등분산 H0 기각x -> 등분산성 만족

## 
##  F test to compare two variances
## 
## data:  t_data1$score and t_data2$score
## F = 3.0787, num df = 9, denom df = 9, p-value = 0.1093
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##   0.7647065 12.3948431
## sample estimates:
## ratio of variances 
##           3.078704

# 1_2. 등분산성 확인후 t-test는 var.equal = T 옵션을 넣어야한다.
t.test( t_data1$score , t_data2$score , var.equal = TRUE) #  0.03199 : 평균차이 유의

## 
##  Two Sample t-test
## 
## data:  t_data1$score and t_data2$score
## t = 2.3247, df = 18, p-value = 0.03199
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   2.213727 43.786273
## sample estimates:
## mean of x mean of y 
##        65        42

# 1_3. 등분산성 옵션을 안주면 결과가 달라진다. t.test default가 등분산성 = FALSE => Welch t-test
t.test( t_data1$score , t_data2$score, var.equal = F) # 0.03531

## 
##  Welch Two Sample t-test
## 
## data:  t_data1$score and t_data2$score
## t = 2.3247, df = 14.289, p-value = 0.03531
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   1.819881 44.180119
## sample estimates:
## mean of x mean of y 
##        65        42

# 1_4. t.test()함수 다른 사용법
# : t.test( 종속변수(관측치) ~ 독립변수(범주칼럼명, 요인), dataframe, 유형(등분산성, 대응표본 등))
t.test( score ~ group, data = t_data, var.equal = T)

## 
##  Two Sample t-test
## 
## data:  score by group
## t = 2.3247, df = 18, p-value = 0.03199
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   2.213727 43.786273
## sample estimates:
## mean in group 1 mean in group 2 
##              65              42

#### 2. 정규성 검정의 3가지 방법 : 0.05보다 커서 H0 기각 안되어야 만족

# 2_1. shapiro-Wilk test 
shapiro.test(t_data1$score) # 0.4604 -> 정규성 만족

## 
##  Shapiro-Wilk normality test
## 
## data:  t_data1$score
## W = 0.93126, p-value = 0.4604

shapiro.test(t_data2$score) # 0.006914 -> H0기각 -> 정규성x -> Mann whitney or Wilcoxon rank-sum?

## 
##  Shapiro-Wilk normality test
## 
## data:  t_data2$score
## W = 0.77358, p-value = 0.006914

# 2_2. nortest 패키지의 ad.test
# install.packages("nortest")
library(nortest)
ad.test(t_data1$score) # 0.4458

## 
##  Anderson-Darling normality test
## 
## data:  t_data1$score
## A = 0.32911, p-value = 0.4458

ad.test(t_data2$score) # 0.006013

## 
##  Anderson-Darling normality test
## 
## data:  t_data2$score
## A = 1.0264, p-value = 0.006013

# 2_3. Kolmogrov-Smirnov test
# ks.test (test, pnorm, mean (test), sd (test))
# pnorm : 누적 정규 분포
?ks.test

## starting httpd help server ... done

ks.test(t_data1$score, "pnorm", mean(t_data1$score), sd(t_data1$score)) # 0.7726

## Warning in ks.test(t_data1$score, "pnorm", mean(t_data1$score),
## sd(t_data1$score)): ties should not be present for the Kolmogorov-Smirnov
## test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  t_data1$score
## D = 0.20947, p-value = 0.7726
## alternative hypothesis: two-sided

ks.test(t_data2$score, "pnorm", mean(t_data2$score), sd(t_data2$score)) # 0.41

## Warning in ks.test(t_data2$score, "pnorm", mean(t_data2$score),
## sd(t_data2$score)): ties should not be present for the Kolmogorov-Smirnov
## test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  t_data2$score
## D = 0.28071, p-value = 0.41
## alternative hypothesis: two-sided

ks.test(t_data1$score, pnorm)

## Warning in ks.test(t_data1$score, pnorm): ties should not be present for
## the Kolmogorov-Smirnov test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  t_data1$score
## D = 1, p-value = 4.122e-09
## alternative hypothesis: two-sided

ks.test(t_data2$score, pnorm)

## Warning in ks.test(t_data2$score, pnorm): ties should not be present for
## the Kolmogorov-Smirnov test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  t_data2$score
## D = 1, p-value = 4.122e-09
## alternative hypothesis: two-sided

ks.test(t_data2$score, pnorm)[2]$p.value

## Warning in ks.test(t_data2$score, pnorm): ties should not be present for
## the Kolmogorov-Smirnov test

## [1] 4.122307e-09

#### 3. 대응표본 T-검정 (실험전/후를 범주로 주어진 상태라 가정)
# 2_1. 대응표본은 정규성 검사가 필수이다.
# 2_1_1. 샤피로 테스트
shapiro.test(t_data1$score) # 0.4604 => H0 기각x => 정규성 만족

## 
##  Shapiro-Wilk normality test
## 
## data:  t_data1$score
## W = 0.93126, p-value = 0.4604

shapiro.test(t_data2$score) # 0.006914 => H0 기각 => 정규성 x => 비모수적 방법인 

## 
##  Shapiro-Wilk normality test
## 
## data:  t_data2$score
## W = 0.77358, p-value = 0.006914

# 2_2. 대응표본일 경우 옵션 paired = T
t.test( t_data1$score, t_data2$score, paired = TRUE) # 0.05106

## 
##  Paired t-test
## 
## data:  t_data1$score and t_data2$score
## t = 2.2493, df = 9, p-value = 0.05106
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1311024 46.1311024
## sample estimates:
## mean of the differences 
##                      23

# 2_3. t.test 다른 사용법
t.test( score ~ group, data = t_data, paired = TRUE) # 0.05106

## 
##  Paired t-test
## 
## data:  score by group
## t = 2.2493, df = 9, p-value = 0.05106
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1311024 46.1311024
## sample estimates:
## mean of the differences 
##                      23

#### 3. 보정된 t검정 : 요인1에 대한 t.test가 끝난 뒤 -> 해당 요인(범주)이외에 다른 변수까지 같이 고려한다.
#### - method로서 선형회귀 -> 2way-ANOVA를 2집단에 적용해서, 분산분석으로 -> 2집단의 평균차이의 유의성을 검증(p-value값이 t-test와 똑같이 나온다.
#### - 그러기 위해선 lm()의 선형회귀로 기울기+절편 결과값을 summary()에 넣어, ANOVA테이블 형태를 얻어야한다.

# 3_1. 1개 요인(t-test와 같은 모델)에 대해 선형회귀->ANOVA 테이블에서 독립표본 t-test와 같은 결과값 얻기
#  선형회귀한 결과값을 분산분석(2way-요인2개 가능)을 2집단에 적용해 -> 평균 차이를 알아본다
mod1 <- lm( score ~ group , data = t_data) 
anova(mod1) # 0.03199 : 아노바 테이블 상 독립표본 t-test와 같은 p 결과값이 나온다.

t.test( t_data1$score , t_data2$score , var.equal = TRUE) # 0.03199

## 
##  Two Sample t-test
## 
## data:  t_data1$score and t_data2$score
## t = 2.3247, df = 18, p-value = 0.03199
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   2.213727 43.786273
## sample estimates:
## mean of x mean of y 
##        65        42

# 3_2. t-test에서 할 수 없는, 새로운 요인(숫자형)을 추가해서 lm->anova 테이블로 p값을 보자.
#      숫자형 요인은 단독으로는 t-test를 할 수 없다(집단 못나눔). 보정시에만 들어감
#      보정변수를 + 앞에 주는 것과 뒤에 주는 것과 결과가 다르다! ㅠㅠ;
#      예제는 추가한 요인 age를 group 앞에 주었다. 그대로 해보자.
mod2 <- lm( score ~ age + group, data = t_data)
anova(mod2) # 0.58633 : age를 앞쪽 요인으로 추가해줬더니, group에 의한 차이는 유의미하지않았다 -> 아! 기존 t-test에서 유의미한 것은 group 때문이 아니라 age땜에 난 것임을 알 수 있다.

# 3_3. 다른 방식 : aov() -> summary()
#      기존 방식 : 기존 lm() -> anova()
summary( aov( score ~ group , data = t_data) ) # 0.032

##             Df Sum Sq Mean Sq F value Pr(>F)  
## group        1   2645  2645.0   5.404  0.032 *
## Residuals   18   8810   489.4                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary( aov( score ~ age + group , data = t_data) ) # 0.586

##             Df Sum Sq Mean Sq F value Pr(>F)  
## age          1   3405    3405   7.322  0.015 *
## group        1    143     143   0.308  0.586  
## Residuals   17   7906     465                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# 3_4. my) 만약 age + group 의 자리를 바꾸면 결과값이 달라진다.
summary( aov( score ~ age + group , data = t_data) ) # 0.586

##             Df Sum Sq Mean Sq F value Pr(>F)  
## age          1   3405    3405   7.322  0.015 *
## group        1    143     143   0.308  0.586  
## Residuals   17   7906     465                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary( aov( score ~ group + age , data = t_data) ) # 0.029 *

##             Df Sum Sq Mean Sq F value Pr(>F)  
## group        1   2645  2645.0   5.687  0.029 *
## age          1    904   903.6   1.943  0.181  
## Residuals   17   7906   465.1                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

연습해보기

• COPD 데이터를 가져와서, lungC(폐암)에 대한 모든 변수(요인, id칼럼만 제외! )들 종속변수로을 자동검정해보기
• mytable( 요인 ~ 종속변수 시리즈 )
• lungC(폐암) 여부(이력, 2집단으로 이루어진 범주)이 BMI(종속변수)에 주는 영향에 대한 t-test(2집단의 평균차이 유의성)확인해보기

# 데이터 준비
getwd()

## [1] "C:/Users/is2js/R_da/myR"

data <- read.csv("week4/4주차_코드/COPD_Lung_Cancer.csv")
head(data)

# 폐암칼럼에 대한 다른칼럼들을 모두 검정해보기
# - 칼럼명이 적히는 곳에 . : 모든 칼럼을 의미 ,  .-칼럼명 : 해당칼럼 제외 모든 칼럼을 의미!!***
library(moonBook)
mytable( lungC ~ .-PERSON_ID, data = data)

## 
##           Descriptive Statistics by 'lungC'          
## ______________________________________________________ 
##                             0            1         p  
##                          (N=870)       (N=74)   
## ------------------------------------------------------ 
##  SEX                                             0.002
##    - 1                 350 (40.2%)   44 (59.5%)       
##    - 2                 520 (59.8%)   30 (40.5%)       
##  DTH                                             1.000
##    - 0                 864 (99.3%)  74 (100.0%)       
##    - 1                  6 ( 0.7%)     0 ( 0.0%)       
##  CTRB_PT_TYPE_CD        5.7 ±  3.4   5.3 ±  3.5  0.307
##  BMI                   22.7 ±  3.3  22.8 ±  2.9  0.710
##  BP_LWST               81.0 ± 12.1  82.4 ± 11.5  0.326
##  BLDS                  102.4 ± 33.2 97.9 ± 26.6  0.178
##  TOT_CHOLE             200.3 ± 39.9 193.1 ± 35.6 0.132
##  HMG                   13.1 ±  1.4  13.2 ±  1.5  0.503
##  FMLY_CANCER_PATIEN_YN                           0.805
##    - 1                 752 (95.4%)   63 (96.9%)       
##    - 2                  36 ( 4.6%)   2 ( 3.1%)        
##  SMK_STAT_TYPE                                   0.654
##    - 1                 670 (80.7%)   54 (77.1%)       
##    - 2                  69 ( 8.3%)   8 (11.4%)        
##    - 3                  91 (11.0%)   8 (11.4%)        
##  DRNK_HABIT                                      0.495
##    - 1                 673 (80.0%)   52 (74.3%)       
##    - 2                  47 ( 5.6%)   6 ( 8.6%)        
##    - 3                  44 ( 5.2%)   6 ( 8.6%)        
##    - 4                  27 ( 3.2%)   1 ( 1.4%)        
##    - 5                  50 ( 5.9%)   5 ( 7.1%)        
##  EXERCI_FREQ                                     0.051
##    - 1                 674 (80.4%)   48 (68.6%)       
##    - 2                  68 ( 8.1%)   10 (14.3%)       
##    - 3                  27 ( 3.2%)   3 ( 4.3%)        
##    - 4                  9 ( 1.1%)    3 ( 4.3%)        
##    - 5                  60 ( 7.2%)   6 ( 8.6%)        
##  copd                                            0.000
##    - 0                 534 (61.4%)   29 (39.2%)       
##    - 1                 336 (38.6%)   45 (60.8%)       
##  Asthma                                          0.006
##    - 0                 422 (48.5%)   23 (31.1%)       
##    - 1                 448 (51.5%)   51 (68.9%)       
##  DM                                              0.012
##    - 0                 335 (38.5%)   40 (54.1%)       
##    - 1                 535 (61.5%)   34 (45.9%)       
##  Tub                                             0.012
##    - 0                 807 (92.8%)   62 (83.8%)       
##    - 1                  63 ( 7.2%)   12 (16.2%)       
##  before_op_score       40.1 ±  6.0  40.1 ±  5.8  0.914
##  after_op_score        50.1 ±  6.0  50.0 ±  6.2  0.907
## ------------------------------------------------------

# t-test 전 준비하기(범주(이력, 여부)별로 데이터 가르기)
data_lung0 <- data[ data$lungC == 0, ]
data_lung1 <- data[ data$lungC == 1, ]

# 1. lungC 여부에 따른 BMI의 등분산성 검사
var.test(data_lung0$BMI, data_lung1$BMI) # 0.1736 -> 등분산성 만족 -> var.equal = T 넣을 것

## 
##  F test to compare two variances
## 
## data:  data_lung0$BMI and data_lung1$BMI
## F = 1.2864, num df = 869, denom df = 73, p-value = 0.1736
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.8933904 1.7632502
## sample estimates:
## ratio of variances 
##           1.286364

# 2. 각각 정규성 검사
shapiro.test( data_lung0$BMI ) # 0.297

## 
##  Shapiro-Wilk normality test
## 
## data:  data_lung0$BMI
## W = 0.99777, p-value = 0.297

shapiro.test( data_lung1$BMI ) # 0.1953 -> 정규성 만족

## 
##  Shapiro-Wilk normality test
## 
## data:  data_lung1$BMI
## W = 0.97698, p-value = 0.1953

# 3. 독립표본 t-test
t.test( data_lung0$BMI, data_lung1$BMI, var.equal = T) # 0.7096 

## 
##  Two Sample t-test
## 
## data:  data_lung0$BMI and data_lung1$BMI
## t = -0.37253, df = 942, p-value = 0.7096
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.9266823  0.6309930
## sample estimates:
## mean of x mean of y 
##  22.70189  22.84973

t.test( BMI ~ lungC , data = data, var.equal = T ) # 0.7096 ( 종속변수 ~ 요인(범주, 이력, 여부))

## 
##  Two Sample t-test
## 
## data:  BMI by lungC
## t = -0.37253, df = 942, p-value = 0.7096
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.9266823  0.6309930
## sample estimates:
## mean in group 0 mean in group 1 
##        22.70189        22.84973

# 4. 대응표본 t-test 
# t.test ( data_lung0$BMI, data_lung1$BMI, paired = T ) # 데이터가 갯수가 달라서, 실험 전/후의 데이터 가 아님



# 5. 보정된 t-test
# 5_1. 요인1개
summary( aov( BMI ~ lungC , data = data)) # 0.71 

##              Df Sum Sq Mean Sq F value Pr(>F)
## lungC         1      1   1.491   0.139   0.71
## Residuals   942  10118  10.741

anova( lm( BMI ~ lungC, data = data)) # 0.7096

# 5_2. 요인 1개 더 추가(숫자형)
summary( aov( BMI ~ TOT_CHOLE + lungC , data = data)) # TOT_CHOLE 0.00121, lung_C 0.60168

##              Df Sum Sq Mean Sq F value  Pr(>F)   
## TOT_CHOLE     1    112  111.62  10.535 0.00121 **
## lungC         1      3    2.89   0.273 0.60168   
## Residuals   931   9864   10.59                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 10 observations deleted due to missingness

summary( aov( BMI ~  lungC + TOT_CHOLE , data = data)) #lung_C 0.71773, TOT_CHOLE 0.00112

##              Df Sum Sq Mean Sq F value  Pr(>F)   
## lungC         1      1    1.39   0.131 0.71773   
## TOT_CHOLE     1    113  113.12  10.677 0.00112 **
## Residuals   931   9864   10.59                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 10 observations deleted due to missingness