전산통계학 14주차 강의 - chisqtest, homogeneity, Independent, nonparametric, ranktest, signtest, Wilcox
코딩 공부/R-전산 통계학 2020. 1. 21. 21:57Non-parametric
우리가 확률론, 수리통계학을 하면서 데이터를 보고 테스팅을 했었는데 그 테스트들은 아마도 parametric (모수 방법) 일 것입니다.
즉, 어떤 일련의 가정이 맞다고 생각하고 진행된 것이므로 엄밀하게 말하자면 앞의 과정이 생략 되있는 셈이에요.
이에 대해서 우리들은 nonparametric (비모수 방법) 을 알아야 할 필요가 있고, 프로그래밍도 할줄 알아야겠죠?
# wilcox (286p.)
x <- c(115.1, 117.8, 116.5, 121, 120.3, 119, 119.8, 118.5)
wilcox.test(x, alternative = "l", mu=120, conf.level = 0.95)$p.value
<첨언>
에러가 나오는 이유는 x에서 120을 뺀 것 중에 랭크가 같은 것이 존재하기 때문입니다,
다시말해서 직접 프로그래밍을 해봄으로써 이러한 변수가 안나오도록 함수를 짜는 것이 꼭 필요할 것 입니다
wil <- function (x, mu) {
n <- length(x)
meanw <- n*(n+1)/4
varw <- n*(n+1)*(2*n+1)/24
y <- x-120
z <- abs(y)
wp <- sum (rank(z)[y>0])
wm <- sum (rank(z)[y<0])
if (wp>wm) {pvalue <- pnorm(wp, mean =meanw, sd=sqrt(varw))}
else {pvalue <- pnorm(wm, mean =meanw, sd=sqrt(varw), lower.tail = F) }
pvalue
}
wil (x, mu=120)
# two group wilcox (353p.)
A <- c(28.6, 25.1, 26.4, 34.9, 29.8, 28.4, 38.5, 30.2, 30.6,
31.8, 41.6, 21.1, 36.0, 37.9, 13.9)
B <- c(69.3, 56.0, 22.1, 47.6, 53.2, 48.1, 23.2, 13.8, 52.6, 34.4, 60.2, 43.8)
wilcox.test(A,B, alternative = "l", paired = F, correct = F)
# wilcox signed rank test (354p.)
A <- c(512, 650, 890, 410, 1050, 1500, 600, 750)
B <- c(500, 600, 890, 400, 1025, 1400, 625, 710)
wilcox.test(A,B,alternative = "two.sided", paired = T, conf.level = 0.90)
# Kruskal.test (553p. ; Nonparametric R.C.D) -> "H0 : M1 = ... = Mk"
x <- c(207,150,197,173,147,144,192,
194,146,175,186,223,143,170,
288,267,288,358,229,249,346,217,203,214)
mat <- factor(rep(1:3,c(7,7,10)))
kruskal.test(x,mat)$statistic
kruskal.test(x~mat)$p.value
<문제>
저희가 가장 많이 보았던 예시인 "514p. 석탄문제" 를 kruskal test 를 통해서 분석해보고, aov 를 통해서도 분석해서 두 값의 차이 구하세요
y<-c(1.51, 1.92, 1.08, 2.04, 2.14, 1.76, 1.17
,1.69, 0.64, 0.90, 1.41, 1.01, 0.84, 1.28, 1.59
,1.56, 1.22, 1.32, 1.39, 1.33, 1.54, 1.04, 2.25, 1.49
,1.30, 0.75, 1.26, 0.69, 0.62, 0.90, 1.20, 0.32
,0.73, 0.80, 0.90, 1.24, 0.82, 0.72, 0.57, 1.18, 0.54, 1.30)
metan <- factor(rep(1:5, c(7,8,9,8,10)))
kruskal.test(y~metan)$p.value ->a
model <- aov(y~metan)
summary(model)[[1]][[5]][[1]] ->b
a-b
<첨언>
지난 강의 때, 다루었던 패키지인 'agricolae' 내에서도 kruskal 명령어가 존재합니다, 그 방법 또한 알아보도록 합시다
하지만, 원소를 뽑아내서 분석하는 것이 번거롭기 때문에 자주 사용되지는 않는 방법입니다...
library(agricolae)
kruskal(y,metan)$statistic[[3]]
kruskal(y,metan)$means[[1]][3]
kruskal(y,metan,console = T)$comparison
하지만, 자료의 자세한 정도는 이쪽이 훨씬 좋음을 볼 수 있습니다.
# Friedman (555p. ; Nonparametric R.C.B.D)
data <- matrix(c(1,4,3,3,
2,1,4,1,
3,3,2,6,
4,2,1,2,
5,5,7,5,
6,6,5,4,
7,8,6,7,
8,7,8,8),nr=4)
block <- factor(rep(c("j1","j2","j3","j4"),len=32))
trt <- factor(rep(1:8,each=4,len=32))
friedman.test(data)$p.value
friedman.test(data,trt,blcok)$statistic
friedman.test(terminal~trt|block)$p.value
<첨언>
역시 friedman test도 앞서 언급한 패키지인 'agricolae' 내에서도 friedman 명령어가 존재합니다, 그 방법 또한 알아보도록 합시다
library(agricolae)
friedman(block,trt,terminal)$statistic[1,3] ; "p.value"
friedman(block,trt,terminal,console = T)
역시, 자세한 정보를 알아낼 수 있음을 확인 가능합니다
# chi-squared goodness (626p.)
x <- c(13,10,42,65,20)
p <- c(0.1, 0.05, 0.25, 0.4, 0.2)
chisq.test(x, p = p)
chisq.test(x, p = p)$p.value
chisq.test(x, p = p)$expect ; "기댓값"
pchisq(5.39,df=4,lower.tail = F)
# independent test (633p.) - example 15.3.4
M <- as.table(rbind(c(75, 60, 65), c(160, 115, 175), c(100, 65, 135), c(15, 10, 25)))
dimnames(M) <- list(Injury = c("None", "Minor", "Major", "Death"), restraint = c("seat belt","both", "None"))
chisq.test(M)$p.value
# homogeneity test (653p.) - column 데이터가 다 똑같다 = 칼럼 확률이 같을 것이다 - example 15.4.4
V <- matrix(c(3,4,13,
5,10,5,
7,11,2),byrow = T,nrow = 3,
dimnames = list("SO_2 level" = c("High", "Normal", "Low"),
"Chloroplast level" = c("High", "Normal", "Low"))
)
chisq.test(V)$p.value
chisq.test(V)$observed
<첨언>
지금까지 나온 비모수 방법들은 명령어를 이용해서 분석을 한 것 뿐입니다, 그 말인 즉슨... 직접 프로그래밍 할 수 있어야겠죠???
6개를 명령어로 만들어보는 연습을 하시면 본인에게 아주 좋을 겁니다~
TAG
