한의대 생활/└ R studio 통계적 분석2

R markdown

데이터 가져오기

setwd("C:/Users/is2js/python_da/deep analysis(논문용 설문지 분석)")
table = read.csv('sens_spec_for_r_table.csv', header = T)
head(table)
##   Kind.of.medicinal.herbs Sensitivity Specitivity              Group
## 1                      SA        0.95       1.000 Ph.D. of Herbology
## 2                      SA        0.95       1.000 Ph.D. of Herbology
## 3                      SA        1.00       0.875 Ph.D. of Herbology
## 4                      SA        0.95       0.900 Ph.D. of Herbology
## 5                      SA        0.90       1.000 Ph.D. of Herbology
## 6                      SA        0.95       0.800 Ph.D. of Herbology

dplyr - filter + 사용자 정의 함수로 약재별로 나누기

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
herb <- function(table, herbname){
  a = filter(table, Kind.of.medicinal.herbs == herbname)
  return(a)
  
}

ac = herb(table,"AC")
amc = herb(table, "AMC")
sa = herb(table, "SA")

Kruskal-Wallis Test

  • 각 약재별로 kruskal.test(민감도칼럼 ~ 그룹칼럼)
  • Kruskal-Wallis test(비모수적 검정)을 사용해야한다. =>
  • H0 : 3그룹의 차이가 없다.
# n수가 적어 정규분포 안따르는 비모수검정 -> 3 그룹간의 차이를 검정
# (정규분포를 따르지 않는)AC에 대한 3그룹 민감도 차이가 없다(H0)를 검정
kruskal.test(ac$Sensitivity ~ ac$Group)  #p-value = 0.02285
## 
##  Kruskal-Wallis rank sum test
## 
## data:  ac$Sensitivity by ac$Group
## Kruskal-Wallis chi-squared = 7.5576, df = 2, p-value = 0.02285
kruskal.test(amc$Sensitivity ~ amc$Group)  #p-value = 0.3188
## 
##  Kruskal-Wallis rank sum test
## 
## data:  amc$Sensitivity by amc$Group
## Kruskal-Wallis chi-squared = 2.2864, df = 2, p-value = 0.3188
kruskal.test(sa$Sensitivity ~ sa$Group)  #p-value = 0.3099
## 
##  Kruskal-Wallis rank sum test
## 
## data:  sa$Sensitivity by sa$Group
## Kruskal-Wallis chi-squared = 2.3433, df = 2, p-value = 0.3099

사후검정 - 각 2군씩로 Mann-Whitney test by bonferroni’s method

  • “pgirmess” 패키지 - kruskalmc( 민감도칼럼, 그룹칼럼)
  • bonferroni’s method란 기존 유의수준(5%)에 nC2 = 2군씩 검정횟수를 나누어 유의수준을 낮추는 것
#install.packages("pgirmess")
library(pgirmess)

kruskalmc(ac$Sensitivity,ac$Group)
## Multiple comparison test after Kruskal-Wallis 
## p.value: 0.05 
## Comparisons
##                                    obs.dif critical.dif difference
## KMD-Ph.D. of Herbology            44.75000     44.74725       TRUE
## KMD-Undergraduates                 2.22973     27.19549      FALSE
## Ph.D. of Herbology-Undergraduates 42.52027     37.51044       TRUE
kruskalmc(amc$Sensitivity,amc$Group)
## Multiple comparison test after Kruskal-Wallis 
## p.value: 0.05 
## Comparisons
##                                      obs.dif critical.dif difference
## KMD-Ph.D. of Herbology            23.7916667     44.74725      FALSE
## KMD-Undergraduates                 0.4087838     27.19549      FALSE
## Ph.D. of Herbology-Undergraduates 23.3828829     37.51044      FALSE
kruskalmc(sa$Sensitivity,sa$Group)
## Multiple comparison test after Kruskal-Wallis 
## p.value: 0.05 
## Comparisons
##                                     obs.dif critical.dif difference
## KMD-Ph.D. of Herbology            26.875000     44.74725      FALSE
## KMD-Undergraduates                 4.600225     27.19549      FALSE
## Ph.D. of Herbology-Undergraduates 22.274775     37.51044      FALSE

#### 실업률 그래프 분석 ####
# http://kosis.kr/index/index.do 
# 국내통계 > 주제별통계 > 고용임금 > 경제활동 인구조사 > 실업률 > 성/연령별 실업률
# csv로 저장후 unemployment_rate.csv

# 1. 데이터 준비
raw_df <- read.csv("data/unemployment_rate.csv", header = T, fileEncoding = 'CP949')
head(raw_df, 20)
summary(raw_df)
image
image

# 2. reshape2(melt)을 통해 좌우로 넓은 데이터 --> long한 통합 데이터 형태로 수정하기
http://nittaku.tistory.com/349

library(tidyr)
library(reshape2)
library(ggplot2)

#    id.var = c("")옵션을 통해 남아있을 칼럼들을 지정해주자.
#    녹는 칼럼은 칼럼명집합칼럼 : variable /  값 집합칼럼 : value로 생성된다.
df_m <- melt(raw_df, id.var = c("성별", "연령계층별"))
head(df_m)
image


# 3. 데이터 클리닝 - substr() -> paste()
# - 지저분한(R에서 자동생성한)하게 녹아있는 칼럼명수정해줘야한다.
# - 1) 년/월이 모여있는 variable칼럼에서 년/월의 숫자만 따로 추출하여 year칼럼 /month칼럼을 생성
# - 2) 따로 추출한 year와 month를 원하는 형태(yyyy-mm)로 합해주는 year_month 칼럼을 만들어 Date타입으로 형변환시킬 준비를 하자.
df_m$year = substr(df_m$variable, 2, 5)
df_m$month = substr(df_m$variable, 8, 9)
head(df_m)

df_m$year_month =paste(df_m$year, df_m$month, sep="-")
head(df_m)
image

# - 3) year_month 칼럼의 값을 as.Date()타입으로 바꿔주자.
# -    하지만, Data()타입으로 변환하기 위해선 -dd(일)의 형태가 반드시 필요하다.
# -    그래서 임의로 -01 을 만들어서 yyyy-mm-dd 형태를 만들어주자.
df_m$year_month = as.Date(paste(df_m$year_month, "01", sep="-"))
head(df_m)
str(df_m)
image
image

# 4) 필요한 칼럼만 챙기자.
df_select = df_m[, c("성별", "연령계층별", "year_month", "value")]
head(df_select)
summary(df_select)

# 5) NA가 있다면,, 제거하자.*** 행인덱싱 자리에 complete.cases( )***
df_cleaned = df_select[ complete.cases(df_select), ]
summary(df_cleaned)
image


# 칼럼명을 영어로.. 임시로 바꾸자.. 반드시x
colnames(df_cleaned) = c("sex", "age_group", "year_month", "value")
head(df_cleaned)


# factor범주형 데이터( sex, age_group)는 summary()에 나타난 순서가 levels이다.(수치형이라면 통계값이 있을 것임)
# levels순으로 범주형 데이터 값들도 levels(칼럼인덱싱) [인덱싱]으로  영어로 바꿔주자.

str(df_cleaned)
levels(df_cleaned$sex)
levels(df_cleaned$sex)[1] = "Total"
levels(df_cleaned$sex)[2] = "Male"
levels(df_cleaned$sex)[3] = "Female"

# age_group은 범주가 너무 많다. levels()로 대입하지말고, gsub()를 이용해 문자열 값을 대체하자.
# 1) 범주가 너무 많을 때는, gsub("발견할문자", "바뀔문자", 칼럼명)을 이용해 일단 간단하게 만들자
#    즉, 문자열 안에 "세" 나 "이상" 등을 지운 새로운칼럼을 하나 생성하자.

levels(df_cleaned$age_group)
df_cleaned$new_age_group <- gsub("세|이상", "", df_cleaned$age_group)
df_cleaned$new_age_group
df_cleaned$new_age_group <- gsub("계", "Total", df_cleaned$new_age_group)


head(df_cleaned)

df_final <- df_cleaned[, c("sex", "new_age_group", "year_month", "value")]
head(df_final)


# ggplot으로 그래프 그리기
- aes()에는 x축의 칼럼, y축의 칼럼 이외에 ~별로 나눌 그룹과, 색을 표시할 그룹을 지정할 수 있다.

ggplot(df_final, aes(x=year_month, y=value, group = sex, colour= sex)) + geom_point() + geom_line()# x축 칼럼, y축칼럼, ~별 그룹, 색으로 지정할 그룹
image

image

       
# 그래프를 못알아 먹겠으니 인덱싱을 통해 df를 연령대 잘라서 각각 그리기
df_20s <- df_final[which(df_final$new_age_group == "20 - 29"), ]
df_30s <- df_final[which(df_final$new_age_group == "30 - 39"), ]
df_40s <- df_final[which(df_final$new_age_group == "40 - 49"), ]
df_50s <- df_final[which(df_final$new_age_group == "50 - 59"), ]

ggplot(df_20s, aes(x=year_month, y=value, group = sex, colour= sex)) + geom_point() + geom_line()# x축 칼럼, y축칼럼, ~별 그룹, 색으로 지정할 그룹
image


# function으로 만들어서 연령대별로 넣기
# data.frame과 title을 받아서 그리기
plotF <- function(df, title){
   ggplot(df, aes(x=year_month, y=value, group = sex, colour= sex)) + geom_point() + geom_line() + ggtitle(title)
}

plotF(df_20s, "20대 실업률")
plotF(df_30s, "30대 실업률")
plotF(df_40s, "40대 실업률")
plotF(df_50s, "50대 실업률")

image

#### F(피셔) test ####
# rbook p354 : 2개의 데이터의 분산-퍼진정도를 비교한다
# 세미나에서는 two sample t-test전에 f - test로 등분산/이분산을 구분해서 계산해야한다고 하였다.


f.test.data <- read.table('data/f.test.data.txt', header = T)
f.test.data

# 1. 각 데이터(칼럼별)의 분산 알아보기
var(f.test.data$gardenB) # 1.3
var(f.test.data$gardenC) # 14.2

#### 2. f-test를 통해 통계학적으로 유믜미적으로 분산의 차이가 있는지 보기***####

# 1) 분산이 큰 값을 분자***로 하여 ***분산비율.ratio 구하기
attach(f.test.data)
F.ratio = var(gardenC) / var(gardenB) # 10.6

# 2) F.ratio의 그래프(F 분포의 확률밀도 함수)에서 하위 2.5% 이하거나 (~95%~) 상위 2.5%(97.5% 이상) = p-value(0.05이하) 에 속하면 유의미한 차이가있다!***
# 즉, pf(F.ratio, n, n) 의 값이  0.25 이하  or 0.975(0.25+0.95)이상일 경우, 분산차이는 유의미하다.
# H0 : 2개의 분산이 다르지 않다.
# 만약 하위2.5% 이하나 상위 2.5%이상에 속한다면( = p-value : 0.05이하), H0거절 H1(분산이 다르다)를 채택image

# 여기서 자유도에 따라 달라진다. 개수가 늘어나면 그래프 봉우리가 오른쪽으로 이동한다.
# 자유도 = 데이터개수 - 1


image
length(gardenB) # 10

# pf()로 F.ration를 이용하여 그래프상(F 분포의 확률밀도함수)의 넓이(density, %) 구하기
pf(F.ratio, 9, 9) # 0.9991879 -> 왼쪽부터 2.5 + 95 = 97.5보다 크니 상위25% 이상으로 분산차이가 유의미하다.


# 참고) qf()로 pf값을(넓이, %)을 이용한  F.ratio기(F-test그래프의 x값) 구하기
qf(0.9991879, 9, 9) #10.66666 -> F.ratio = 그래프상의 x 구하기


# pf값은 그래프상의 왼쪽부터 넓이다. 전체넓이인 1에서 빼주면, 오른쪽 나머지 넓이 상위 몇퍼센트인지 알 수 있다.
# 이것에 2를 곱하면, 양쪽 넓이가 나올 것인데, 이것이 5% 0.05보다 작아야 유의미한것이다.
# p-value에 대해 계산하는 것

2 * ( 1- pf(F.ratio, 9, 9)) # 0.0016 < 0.05  ==> H0를 거절하고 분산이 서로 유의미하게 차이가 난다.

# **** 한번에 F-test하는 방법****
# F.ratio 값만 다르다. 왜냐하면 막.. 순서대로.. 앞에 것을 분자로 넣기 때문에..(원래는 큰값이 분자로)
# 그러나 p-value값으로 비교하는 것이 최종목적이라 상관은 없다고 한다.

var.test(gardenB, gardenC)
# F test to compare two variances

# data:  gardenB and gardenC
# F = 0.09375, num df = 9, denom df = 9, p-value = 0.001624
# alternative hypothesis: true ratio of variances is not equal to 1
# 95 percent confidence interval: 0.02328617 0.37743695
# sample estimates: ratio of variances  0.09375

### R-squared의 종류 ####
#점들이 퍼져있는 정도, 분산의정도
# multiple ->변수와 관계없는 R-squared / adjusted-> x값이 늘어날수록 수치가 줄어드는 R-suared
year <- c(2000:2004)
value <- c(2.3, 3.2, 5.6, 5.4, 5.8)

par(mfrow=c(1,1))
plot(year, value)

fit <- lm(value ~ year)
abline(fit, col="red")

summary(fit)# R-squared 0.82/0.76
image


# 좀 더 선형에 가깝게 3번째 값 조절 -> R-squared 값이 올라간다.
year <- c(2000:2004)
value <- c(2.3, 3.2, 4.6, 5.4, 5.8)

par(mfrow=c(1,1))
plot(year, value)

fit <- lm(value ~ year)
abline(fit, col="red")

summary(fit)# R-squared 0.96/0.95
image


#### 2. plot(fit - lm()결과 ) ####
# *** 원래는 기본plot() 위 에다가  abline( lm결과 )으로 넣어주면 lm선으로 나왔음.
# 그래프 4개를 return한다. -> 콘솔에서 enter쳐야함
plot(fit)
# 4개를 par(mfrow)를 이용해 한번에 subplots처럼 그리기
par(mfrow = c(2,2))
plot(fit)

# 1번째 그래프 : 잔차그래프 -> 0에가깝에 수렴할 수록 좋은데 왔다갔다함
# 2번째 그래프 : QQplot( 정규분포의 정도를 보는 )
# 3번째 그래프 : 잔차에 루트씌운 그래프
# 4번째 그래프 : outlier처럼 영향을 크게 주는 데이터를 확인할 수 있다. -> Cook's distance에 가깝거나 포함되는 outlier 데이터
image


# 선형모델로 설명된다고 하지만,, 왔다갔다 패턴을 가진 예1
par(mfrow = c(1,1))
x=c(1, 2, 3, 4, 5, 7, 8, 9, 10)
y=c(2, 1, 4, 3, 6, 5, 8, 7, 9)
plot(x, y) #올라갔다가 내려갔다가 하는 패턴이 있음. # 어떤 변수를 놓쳤다는 것을 의미할 수도.. 날씨-계절변수
fit <- lm(y~x)
abline(fit, col="red")

summary(fit) # 별표있고, R-squared 높고, p-value가 낮으나 잔차가 왔다갔다하므로 < 좋은 모델이라고 할 수 없을 것이다.. >
image

# *** plot(fit)
# 잔차가 왔다갔다한다 -> 좋은그래프가 아니다.
# 이런식으로 패턴이 있는 데이터 -> 좋은 모델이 아니다.
# 그러나 확인할 수 있는 것에는 괜찮으니 --> 놓친 변수를 차자.

par(mfrow = c(2,2))
plot(fit)
image



# 선형모델로 설명된다고 하지만,, 왔다갔다 패턴을 가진 예2
par(mfrow = c(1,1))
x=c(1, 2, 3, 4, 5, 7, 8, 9, 10)
y=c(2, 3, 4, 5, 0, 2, 3, 4, 5)
plot(x, y) #올라갔다가 내려갔다가 하는 패턴이 있음. # 어떤 변수를 놓쳤다는 것을 의미할 수도.. 날씨-계절변수
image
fit <- lm(y~x)
abline(fit, col="red")
image

summary(fit)
# 별표 x
# 점이 퍼진정도가 심하여 R-squared가 매우 작다.
# p-value가 매우커서 선형모델로 설명x

par(mfrow = c(2,2))
plot(fit)
# 잔차가 왔다갔다 한다.
# QQ는 괜찮아보이지만,, 잔차때문에.. 이 선형모델은 버려야한다.


# iris데이터에서 1종의 Species만 뽑아  length에 따른 width의  선형모델 알아보기
data(iris)
head(iris)

length = iris[which(iris$Species=="setosa"), ]$Sepal.Length
witdh = iris[which(iris$Species=="setosa"), ]$Sepal.Width
par(mfrow = c(1,1))
plot(length, witdh, col="blue")
image

# ***overlab된 데이터를 살짝 흔들고 싶을 땐 x변수에 jitter()
plot( jitter(length), witdh, col="blue")

# 선형모델 만들기
fit <- lm(witdh ~ length)
abline(fit, col="red")

image

#확인하기
summary(fit)
par(mfrow=c(2,2))
plot(fit) # 완벽한 1자는 아니더라도 잔차가 0에 가까운 선을 보인다.
# ***4번째 그래프를 보니 42번째 데이터가 영향을 많이 준다 = outlier일 것이다.

image


# ***boxplot()로 y값 그려 outlier확인해보기
#참고 : outlier확인 시각화 / http://nittaku.tistory.com/353
par(mfrow=c(1,1))
boxplot(witdh) # boxplot에는 y값만 넣기***! 하나만 티어나와있다..
image

# *** 아웃라이어 제거해보기
length_new = length[-42]
witdh_new = witdh[-42]
boxplot(witdh_new)
image

plot(jitter(length), witdh, col="blue")
# 비교
plot( jitter(length_new), witdh_new, col="blue")
image
# lm라인 추가해보기
fit2 <-lm(witdh_new ~ length_new)
abline(fit2, col="red")

# 기존lm라인을 노란색으로 추가해보기
abline(fit, col="yellow")

image

# fit결과 전체 그리기 *** 원래는 abline으로 넣어주면 선으로 나왔음.
par(mfrow=c(2,2))
plot(fit2) # cook's distance에 다가가는 점없어짐.
image


#### 선형회귀모델 최종 확인 ####
# 변수에 별표 / R-Square / p-value  확인
# plot(fit)으로 잔차의 패턴 확인
# 왔다갔다한다면 , 선형모델로 나오더라도 다른 변수 집어넣기

#### *** 변수 추가해보기 *** ####
x=c(1, 2, 3, 4, 5, 7, 8, 9, 10)
y=c(2, 1, 4, 3, 12, 5, 8, 7, 9)
z=c(2, 1, 4, 3, 12, 5, 8, 7, 10)
# lm(z ~ x + y) 형태로 넣기 나머지 2개는 +로 연결해줌. 더하기는 아님 ****
fit3 <- lm(z~ x+y)
summary(fit3)  # y변수에만 별표 3개이다. x는 버려도 될 듯

par(mfrow=c(2,2))
plot(fit3)

# 별표없는 x를 지우고 해보자..
fit3 <- lm(z~ y)
plot(fit3) # 9번째 잔차만 많이 티어나와있다.

#### 부트스트래핑 ####
# ****쓰는 이유 : 많은 test들이 정규성을 가정하고 쓴다(student's t) 하지만,
부트스트레핑을 쓰면, 적은 데이터라도 정규분포를 형성시켜 --> 모집단의 평균을 추정한다.
# https://thebook.io/006723/ch10/08-2/
# 부트스트래핑은 주어진 데이터로부터 복원 표본을 구하는 작업을 여러 번 반복해 원하는 값을 추정한다. 한 가지 예로 평균의 신뢰 구간을 구하는 경우를 생각해보자. 평균의 신뢰 구간은 평균이 속할 범위를 95% 신뢰도로 찾는 것이다. 부트스트래핑은 데이터로부터 표본을 복원 추출로 구하고 이들의 평균을 기록하는 일을 반복한다. 이렇게 구한 평균들을 나열한 뒤 가장 작은 값으로부터 2.5%, 가장 큰 값으로부터 97.5% 지점의 평균 두 개를 구한다. 그러면 두 평균 사이의 구간이 평균의 95% 신뢰 구간이 된다.

# 지난시간 정규성이 없어보이던 데이터
light <- read.table("data/light.txt", header = T) # 하위폴더에 접근할때는 앞에 / 안붙히고 바로..
light

hist(light$speed, col = "green")
image

# QQplot으로 확인하는 정규분포
qqnorm(light$speed)
qqline(light$speed, col = "red")
shapiro.test(light$speed) # p-value = 0.09876 -> H0정규분포를 채택하긴 했으나 찜찜..
image


#### 부트스트래핑 ####
# 0) 데이터의 개수 알아놓기
length(light$speed) # 데이터 개수 : 20개


# 1) 몇번 쌤플링(추출)할 것인지를 --> 빈 숫자 데이터에 정해놓기
a <- numeric(10000) # 10000개의 빈 숫자데이터 생성 -> 샘플링할 횟수가 될 것임

# 2) 20개씩 복원추출(1개뺀뒤 넣고 다시뺌x20)한 sample의 분포를 히스토그램으로 한번 봐보기
hist( sample(light$speed, size=20, replace = T), col="green" )
hist( sample(light$speed, size=20, replace = T), col="red" )
hist( sample(light$speed, size=20, replace = T), col="blue" )
imageimageimage

# 3) 10000개의 빈공간에 < 각 sample들( - *20개를 - 복원추출(replace)한 것) 의 평균mean*** >을 집어넣음
  my) 데이터 20개짜리에 대해,  10000번을,  20개씩 복원추출하여 나온 각 평균을  for문을 통해서 집어넣는다.

# 20개를 복원추출한 샘플의 평균구해보기
mean(  sample(light$speed, size=20, replace = T))

for( i in 1:10000){
   a[i] <- mean(  sample(light$speed, size=20, replace = T)  )
}

# 4) 이제 10000개의 20개복원추출의평균  을  히스토그램으로 분포를 보자.
hist(a, col='yellow') # 정규분포를 띈다.

# 5) 10000개짜리 정규분포의 평균으로 -> 모집단의 평균(H0는 990일 것이다)를 추정하여 H0를 reject하거나 채택가능

# 6) 각 샘플들에서는 1000짜리도 등장하지만, 그것(20개-복원추출)의 평균한 것들은 최대값도 990(H0)에는 못간다.
#    결과적으로 빛의 속도 990이라는 H0는 reject할 수 있다.
#    자료가 충분하지 않을 때는, 부트스트래핑을 통해 데이터를 잔뜩 늘일고 + 정규분포를 만들어서-> 모집단의 평균을 추정한다.
#    중심값 극한 정리에 쓰인다?
#    1개의 샘플의 평균은 별 의미가 없으나 1000개 이상한 부트스트레핑의 평균은--> 모집단의 평균을 추정가능하다.
#  
max(a) # 978.5
mean(a) # 909.2256 --> 1. 윌콕슨검정 시간에 909은 H0를 채택했었는데, 그 값과 비슷

# 데이터 다운로드 : http://www.bio.ic.ac.uk/research/mjcraw/therbook/
# 책 pdf 주소 : https://www.cs.upc.edu/~robert/teaching/estadistica/TheRBook.pdf
# chapter 8 classic test 부터 다룬다.

# 데이터 준비
getwd()
# header = T를 안붙힐 경우, txt에서온 칼럼명이 row로 들어가고 V1이라는 칼럼명이 생성된다.
light <- read.table("data/light.txt", header = T) #폴더에 접근할때는앞에 / 안붙히고 바로..
light


#### normal QQ plot ####
# 2개의 그래프를 그리고 접선과 일치할수록, 정규분포를 따르는 데이터라 할 수 있다.

# 1) 눈 대중..으로 데이터의 분포를 보기 위한 hist() 히스토그램확인
hist(light$speed, col = "green")
image


# 2) rnorm()으로 생성한 정규분포 데이터로 그래프 비교해보기
x <- rnorm(1000)
hist(x)

# 3) 데이터 통계를 보기 위한 summary() + boxplot()****
summary(light) # 통계정보
boxplot(light) # 통계정보를 그래프로
image

# 4) hist()그램보다 더 정확하게 정규분포인지를 확인할 수 있는 QQplot *** ####
# qqnorm() : 데이터들을 점으로 표시하여, qqline()과 비교할 준비를 한다.
# qqline() : quantile을 기준으로 줄을 세워, 정규분포를 가정하는 선을 생성
# 생성된 qqline과 점들의 분포가 비슷하면, normal distribution을 "부정하기 어렵다"라고 말함.

# 1단계 : H0 : 정규분포 되어있다. 를 가정
# 2단계 : qqnorm()+qqline()-> 그려보고 선에 안맞으면 reject할 준비
# 3단계 : 실제적으로 shapiro.t()를 넣어 p-value를 확인하여 0.05(5%)보다 큰 것을 확인한다.
 

# H0 : 정규분포를 가진다

qqnorm(light$speed)
qqline(light$speed, col = "red")
# 선에 벗어나 있으므로 H0 거절할 준비하기
image

# p-value확인
shapiro.test(light$speed) # (대립가설이)0.09를 가져 0.05보다 크다. => H0를 채택해야하는 상황.
# 정규분포라 봐야한다...


# 5) 완전 정규분포인 데이터를 가지고 qqplot을 그려보자.
# 1단계
# H0 : 정규분포 되어있다. 를 부정할 수 없을 것이다..예상
x <- rnorm(1000)

# 2단계
qqnorm(x)
qqline(x, col = "red")
# 정규분포 선과 점들이 거의 일치
image

# 3단계
shapiro.test(x) # p-value가 0.05보다 커서 H0채택 H1거절 -> 정규분포   아니라고말할수x
# cf) 카이제곱검정 H0 : chisq.test()로 검정 : 남/녀가 영향안줄 것이다. 결과들 독립이다. 를 검정
#      - p가 0.05보다 작으면 H0거절 H1(남여영향줌)을 채택
# cf) 정규분포검정 H0 : shapiro.test()로 검정정규분포를 따를 것이다.
#      - p가 0.05보다 작으면 H0거절 H1(정규분포라고 말할 수 없다.)을 채택
      

#  다시 한번 정규분포x인 데이터로 연습해보자.
x <- exp(  rnorm(1000)  )

# 1단계
# H0 : 정규분포 되어있다. 를 부정(reject)


# 2단계
qqnorm(x)
qqline(x, col = "red")
# 정규분포 선과 점들이 일치x
image

# 3단계
shapiro.test(x) # p-value가 엄청 작다 -> H0 거절 H1채택 -> 정규분포 라고 보기 힘듬.



# 7) wilcox.test()를 통해, 정규분포가 아닌 측정된 샘플데이터에서 [H0:빛의 속도(모집단)의 평균(mu)는 299 990 k/ms]이다.를 검정해보자
# ***정규분포를 따르지 않으므로 Student's t test(데이터가 정규분포일 때, 모집단의 평균 검정)는 하지 못하여 wilcox.test()를 사용한다.
# 데이터 설명
# 빛의 속도라고 알려진 299990에서 299000을 뺀 수치들이다.
# 그 결과 정규분포를 따르지 않게 되고
# 정규분포를 가정하고 실험하는 Student's t test를 하지 못한다.***

#*** 윌콕슨 검정 H0 : 샘플로 본 빛의 속도(모집단)의 평균은 990 -> 299 990일 것이다.
wilcox.test(light$speed, mu = 990 )
# 그 결과 p-value가 0.05보다 작다 -> H0:거절 / H1 채택-> 빛의 속도평균 299990 아님

# summary()로 확인한 샘플들의 Mean값을 가지고 윌콕슨 검정해보자.
summary(light$speed) # Mean = 909
#*** 윌콕슨 검정 H0 : 샘플로 본 빛의 속도(모집단)의 평균은 909 -> 299 909일 것이다.
wilcox.test(light$speed, mu = 909 )
# 그 결과 p-value가 0.69 -> H0: 채택-> 빛의 속도평균 299 909 이 아니라고할수x

+ Recent posts