R, Nonparametric test program example (15' LSG & LSM)

# 1. SIGN Test - 1Group (correct)

한 그룹에서의 싸인 검정을 할 때의, 중요한 점은 다른 강의나 과제에서 공부했듯이, 'Handling Zeros in test' 를 잘 숙지해야 한다는점

따라서, 각 방향 (side 를 말함) 과 검정할 때 본인이 사용할 'Q star (test statistic 을 말함)의 부호' 에 따라서 equal data 의 카운트가 되는지 안되는지를 잘 인지해야 한다는 점이다 

library(DescTools)

# Data
DATA <- c(35,65,48,40,70,50,58,36,47,41,49,39,34,33,31
             ,55,55,55,55,55)

# f.t
sm.st <- function(x, mzero, side) {
  equal=0 ; low=0 ; high=0

  for (i in 1: length(x)) {
    if(x[i]==mzero) {equal=equal+1} 
    else if (x[i]<mzero) {low=low+1} 
    else {high= high+1} 
  }
  qstar.a <- low           ; qstar.c <- high
  qstar.b <- low+equal     ; qstar.d <- high+equal

  if (side=="two.sided") {
    if(low<high) {pvalue <- min((1- pbinom(qstar.b-1, length(x),0.5))*2, pbinom(qstar.b, length(x),0.5)*2)}
  else {pvalue <- min((1- pbinom(qstar.a-1, length(x),0.5))*2, pbinom(qstar.a, length(x),0.5)*2)}
  }

  if (side=="less") 
    {pvalue <- 1-pbinom(qstar.a - 1, length(x), 0.5)}

  else if (side=="greater") 
  {pvalue <- pbinom(qstar.b , length(x), 0.5)}

  if (pvalue>1) {pvalue <- 1} 
  cat("\n",side ,"of Sign-Test","\np-value =",pvalue)
}

sm.st(DATA, 55, "two.sided")
sm.st(DATA, 55, "less")
sm.st(DATA, 55, "greater")

---

# 2. SIGNED Rank Test - 1Group (correct)

한 그룹에서의 부호 검정을 실시할 때는, 1에서 다룬 'SIGN Test' 와 비슷하게 equal data를 다루되, 'W star & mean' 이 두 값이 어떻게 결정되는지만 잘 알고있다면, 까다로울 것이 없다

# Data
x <- c(115.1, 120.3, 117.8, 119, 116.5, 119.8, 121, 118.5)
wilcox.test(x, alternative = "g", mu=120, conf.level = 0.95)

# f.t
sm.srt1 <- function (x,M0,side) {
  n <- length(x) ; d <- x-M0 
  rs <- rank(abs(x-M0),ties.method = "a")*sign(x-M0)
  wp <- sum(rs[rs>0]) ; wm <- sum(abs(rs[rs<0]))
  
  we <- 0
  for(i in 1:n) {if(d[i]==0) {we <- we+rank(abs(d))[i]}}
  
  wstar.n <- wm ; wstar.e <- wm+we
  mean=n*(n+1)/4 ; sd=sqrt(n*(n+1)*(2*n+1)/24)
    
    if(side=="two.sided") {
      if (wm<wp) {pv <- min(pnorm(wstar.e,mean,sd,lower.tail = T)*2, pnorm(wstar.e,mean,sd,lower.tail = F)*2)}
      else {pv <- min(pnorm(wstar.n,mean,sd,lower.tail = T)*2, pnorm(wstar.n,mean,sd,lower.tail = F)*2)}}
    if(side=="less") {pv <- 1-pnorm(wstar.n,mean,sd)}
    if(side=="greater") {pv <- pnorm(wstar.e,mean,sd)}
  if (pv>1) pv <- 1 ; cat("1Group - case of",side,"\np-value =",pv)
}

sm.srt1(x,M0=120,side="greater")

---

# 3-1. SIGNED Rank Test Paired Data - 2Group To way.1 (주의! 아직 필자도 정확히 완성시키지를 못함...)

단순하게 두 개의 그룹에 대해서 그 차이를 d로 칭한 후, 1번 (SIGN Test) 대로 하면 된다고 생각하는데, 그 차이에 대한 랭크(순위)와 mu 기준에 대한 부호 카운트가 조금 까다로웠던 것으로 생각난다

# DATA
x=c(512,650,890,410,1050,1500,600,750,850)
y=c(500,600,890,400,1025,1400,625,710,850)
wilcox.test(x,y,alternative = "l", paired = T, conf.level = 0.90, correct = T)
d=x-y ; SignTest(d, alternative = "l" ,mu=0)

# f.t
sm.srtp <- function (x,y,side) {
 n=length(x) ; rs=rank(abs(x-y))*sign(x-y)
  wm <- length(which(rs<0)) ; wp <- length(which(rs>0)) ; we <- length(which(rs==0))
  wstar.a <- wm        ; wstar.c <- wp
  wstar.b <- wm+we     ; wstar.d <- wp+we
    
     if (side=="two.sided") {
    if(wm<wp) {pvalue <- min((1-pbinom(wstar.c-1, n,0.5))*2, pbinom(wstar.c, n,0.5)*2)}
  else {pvalue <- min((1-pbinom(wstar.d-1, n,0.5))*2, pbinom(wstar.d, n,0.5)*2)}
  }
  
  if (side=="less") {pvalue <- pbinom(wstar.d-1, n, 0.5)}
  
  else if (side=="greater") {pvalue <- 1-pbinom(wstar.c , n, 0.5)}
    
    if(pvalue>1) pvalue=1 ; cat("2Group - case of",side,"\np-value =",pvalue)
}

sm.srtp(x,y,side="less")

---

# 3-2. SIGNED Rank Test Paired Data - 2Group To way.2

mean, sd(=root(var)) 만을 이용하면 되기 때문에, 3-1보다는 까다로울 점이 없다고 생각됩니다

# DATA
x=c(512,650,890,410,1050,1500,600,750,850)
y=c(500,600,890,400,1025,1400,625,710,850)
wilcox.test(x,y,alternative = "g", paired = T, conf.level = 0.9, correct = T)

# f.t
sm.srtp <- function (x,y,side) {
  d <- x-y ; z <- abs(d) 
  n=length(x) ; rs=rank(abs(x-y))*sign(x-y)
  wp <- sum(rs[rs>0]) ; wm <- sum(abs(rs[rs<0]))
  
  we <- 0
  for (i in 1:n) {if (d[i]==0) {we <- we+rank(z)[i]}}
  
  wstar.n <- wm ; wstar.e <- wm+we
  
  mean=n*(n+1)/4 ; sd=sqrt(n*(n+1)*(2*n+1)/24)
  
  if(side=="two.sided") {
    if (wm<wp) {pv <- min(pnorm(wstar.e ,mean,sd,lower.tail = T)*2, pnorm(wstar.e,mean,sd,lower.tail = F)*2)}
    else {pv <- min(pnorm(wstar.n,mean,sd,lower.tail = T)*2, pnorm(wstar.n,mean,sd,lower.tail = F)*2)}}
  
  if(side=="less") {pv <- 1-pnorm(wstar.n,mean,sd)}
  if(side=="greater") {pv <- pnorm(wstar.e,mean,sd)}
  if (pv>1) pv <- 1 ; cat("2Group - case of",side,"\np-value =",pv)
}

sm.srtp(x,y,side="greater")

---

# 4. SIGNED Rank Sum Test - 2Group (correct)
3번보다는 식을 구성하는데에 훨씬 쉬웠을 거에요, 하지만 주의해야 할 점이 '두 그룹중에 원소개수가 상대적으로 적은 것을 t.s 로 이용'

해야함을 인지하고 있어야 합니다

# DATA
x <- 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)
y <- 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(x,y, alternative = "l", paired = F, correct = F, conf.level = 0.95)

# f.t
sm.srst <-function(x,y,side){
  all <- c(x,y) ; nx <- length(x) ; ny <- length(y)
 if (nx<=ny) {m <- nx ; n <- ny}
  else {m <- ny ; n <- nx}
  r <- rank(abs(all))
  ifelse (nx<=ny, wm <- sum(r[1:nx]), wm <- sum(r[(nx+1):(nx+ny)]))

  mean=m*(m+n+1)/2 ; sd=sqrt(m*n*(m+n+1)/12)
    
    if(side=="two.sided") {pv <- min(pnorm(wm,mean,sd,lower.tail = F)*2,pnorm(wm,mean,sd,lower.tail = T))*2}
    if(side=="less") {pv <- pnorm(wm,mean,sd,lower.tail = F)}
    if(side=="greater") {pv <- pnorm(wm,mean,sd,lower.tail = T)}
    
  cat("2Group - case of",side,"\np-value =",pv)
}
sm.srst(x,y,side="less")

---

# 5. Kruskal Test (R.C.D - nonparametric)
크루스칼에서 까다로웠던 점은 그나마 '데이터 리스트 형태에서 어떻게 랭크를 매기며, 같은 리스트에 속하는 원소의 랭크끼리 어떻게

합하는가?' 정도입니다

# DATA
z=list()
z[[1]] <- c(207,150,197,173,147,144,192)
z[[2]] <- c(194,146,175,186,223,143,170)
z[[3]] <- c(288,267,288,358,229,249,346,217,203,214)

ori <- unlist(z)
mat <- factor(rep(1:3,c(7,7,10)))
kruskal.test(ori,mat)

# f.t
kru <- function (x) {
  k<-length(x)
  n<-c()
  t<-c()
  
  y <- rank(unlist(x), ties.method = "average")
  M <- c()
  for (m in 1:k) { M[m] <- length(x[[m]]) }
  x <- list()
  for (m in 1:k) {
    ifelse (m==1, x[[m]] <- c( y[1:M[m]] ),x[[m]] <- c( y[(sum(M[0:(m-1)])+1):sum(M[1:m])] ))
  }
  
  for (i in 1:k) {
    n[i] <- length (x[[i]]) 
    t[i] <- sum(x[[i]])
    Ti <- sum (t^2/n)
  }
  N<-sum(n) 
  
  # d.f
  hstar <- 12/(N* (N+1)) *Ti -3*(N+1)
  
  pv <- pchisq(hstar, k-1, lower.tail = F)
  
  cat("test statistic =",hstar, "p-value =",pv)
}
kru(z)

---

# 6. Freidman Test (R.C.B.D - nonparametric)

사실 'tapply' 쓰시면 5번처럼 번거롭게 하실필요가 없어요.... ㅎㅎ

(왜냐하면, tapply 자체는 필자가 정의한 섹터에 대해서 지정한 연산을 진행하므로, 5번에서 걱정했던 점을 한번에 해결합니다)

# DATA
data=list()
data[[1]]=c(1,2,3,4,5,6,7,8)
data[[2]]=c(4,1,3,2,5,6,8,7)
data[[3]]=c(3,4,2,1,7,5,6,8)
data[[4]]=c(3,1,6,2,5,4,7,8)
trt <- factor(rep(1:8, len=32))
block <- gl(4,8)

# f.t
fri <- function(x,treatment,block){
  y <- unlist(x)
  k <- nlevels(treatment) ; b <- nlevels(block)
  
  for (i in 1:k) {Ti <- tapply(y,treatment,sum)}
  sstar <- 12/(b*k*(k+1))*sum(Ti^2) -(3*b*(k+1))
  
  pv <- pchisq(sstar,k-1,lower.tail = F)
  
  cat("test statistic =",sstar, "p-value =",pv)
}

fri(data,trt,block)

# Original
library(agricolae)
z <- 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)
judge <- factor(rep(1:4, len=32))
ter <- factor(rep(1:8,each=4))
friedman.test(z~ter|judge)