CI CV

Confidence Interval of Cramer’s V

Cramer’s V の信頼区間を求める方法は３つ知られている。

Fisher’s Z 変換を用いる　
近傍data を作成して bootstrap
Non central chisq 分布を用いる

この３通りのいずれもRで実行できる。各々次の通り

userfriendlyscience::confIntV(x,y=NULL,conf.level=0.95)
RVAideMemoire::cramer.test(x,y,nrep=1000,conf.level=0.95)
DescTools::CramerV(x,y=NULL,conf.level=NA)

なお、第一の confIntV は、bootstrap にもつかえるはずえあるが、現状ではうまく作動しない。また、あるサイトに 3 の方法とされるスクリプトがあるが、以下の検算とは一致しない。

以下、いくつかの検算を行う。

Fisher’s Z transformation

Fisher’s Z 変換とは、２次元正規分布に従う分布からのn点標本の相関係数 r を逆双曲線正接関数で変換したz=atanh(r) を取り扱う。それが分散 \(1/\sqrt{n-3}\) の正規分布にしたがうとして信頼区間等を求め、tanh(z) により元に戻す。

これを、行列 mat の Cramer’s V に対して \(n=sum(mat)\) として適用する。実際のスクリプトはたとえば次の通り。

fzV <- function(x,y,biasadj=FALSE){
  #"mat" being a r x c matrix/table
  if (is.matrix(x) || is.data.frame(x)){
    mat <- as.matrix(x)
  } else {
    mat <- table(x,y)
  }
  chicalc <- chisq.test(mat,correct = FALSE)
  # calculate Cramer's v -
  K <- min(nrow(mat),ncol(mat))
  crv <- sqrt(chicalc$statistic / (sum(mat)*(K-1)))
  names(crv) <- "CramersV"
  # convert the Cramer's V to a Fisher's Z
  fz <- atanh(crv)
  
  # calculate 95% conf.int around Fisher Z
  conf.level <- 0.95
  bias=0
  if (biasadj) 
    bias <- 0.5*crv/(sum(mat)-1)
  se <- 1/sqrt(sum(mat)-3) * qnorm(1-(1-conf.level)/2)
  cifz <- fz + c(-se,se) + bias
  # convert it back to conf.int around Cramer's V
  cicrv <- tanh(cifz)
  if (cicrv[1]<0) cicrv[1] <- 0
  list(crv=crv, cicrv=cicrv, p.value=chicalc$p.value)
}

実際に適用してみるため、いくつかのデータを設定する。

library(userfriendlyscience)

## Loading required package: ggplot2

# library(RVAideMemoire)
library(DescTools)

## Loading required package: manipulate

HEC12 <- HairEyeColor[,,1] + HairEyeColor[,,2]
yield <- read.table(
text="fertilizerBlend Wheet Corn Soy Rice
X1               X   123  128 166  151
X2               X   156  150 178  125
X3               X   112  174 187  117
X4               X   100  116 153  155
X5               X   168  109 195  158
Y1               Y   135  175 140  167
Y2               Y   130  132 145  183
Y3               Y   176  120 159  142
Y4               Y   120  187 131  167
Y5               Y   155  184 126  168
Z1               Z   156  186 185  175
Z2               Z   180  138 206  173
Z3               Z   147  178 188  154
Z4               Z   146  176 165  191
Z5               Z   193  190 188  169")
yield2 <- as.matrix(yield[,-1])
require(MASS)

## Loading required package: MASS

data(caith)
dat2 <- matrix(c(805, 414, 226, 99, 38, 12), nrow = 3)
rownames(dat2) <- c("Level 1", "Level 2", "Level 3")
colnames(dat2) <- c("Correct", "Error")

適用すれば、次の通りであり、confIntV の結果と一致する。また、バイアス補正はされていないことがわかる。

fzV(HEC12)

## $crv
##  CramersV 
## 0.2790446 
## 
## $cicrv
## [1] 0.2030262 0.3517195
## 
## $p.value
## [1] 2.325287e-25

confIntV(HEC12)

## Cramer's V 95% confidence interval (point estimate = 0.279):
## Using Fisher's z: [0.203, 0.3517]

fzV(yield2)

## $crv
##   CramersV 
## 0.07558305 
## 
## $cicrv
## [1] 0.05550204 0.09560294
## 
## $p.value
## [1] 5.521135e-16

confIntV(yield2)

## Cramer's V 95% confidence interval (point estimate = 0.07558):
## Using Fisher's z: [0.0555, 0.0956]

Noncentral chi-sq distribution

非心カイ２乗分布を用い、Cramer’s V の元となっているカイ２乗値の信頼区間を求め、それに自由度を加えて適用する。実際のスクリプトの例としてはたとえば次の通り。ここで lochi, hichi は既存のものを用いる。出典 Chiscript.SSC http://psychology3.anu.edu.au/people/smithson/details/CIstuff/CI.html

lochi <-
  function(chival,df,conf)
  {
    ulim <- 1 - (1-conf)/2
    #   This first part finds a lower value from which to start.
    lc <- c(.001,chival/2,chival)
    while(pchisq(chival,df,lc[1])<ulim) {
      if(pchisq(chival,df)<ulim)
        return(c(0,pchisq(chival,df)))
      lc <- c(lc[1]/4,lc[1],lc[3])
    }
    #   This next part finds the lower limit for the ncp.
    diff <- 1
    while(diff > .00001) {
      if(pchisq(chival,df,lc[2])<ulim) 
        lc <- c(lc[1],(lc[1]+lc[2])/2,lc[2]) 
      else lc <- c(lc[2],(lc[2]+lc[3])/2,lc[3])
      diff <- abs(pchisq(chival,df,lc[2]) - ulim)
      ucdf <- pchisq(chival,df,lc[2])
    }
    c(lc[2],ucdf)
  }
#
#
hichi <-
  function(chival,df,conf)
  {
    #   This first part finds upper and lower startinig values.
    uc <- c(chival,2*chival,3*chival)
    llim <- (1-conf)/2
    while(pchisq(chival,df,uc[1])<llim) {
      uc <- c(uc[1]/4,uc[1],uc[3])
    }
    while(pchisq(chival,df,uc[3])>llim) {
      uc <- c(uc[1],uc[3],uc[3]+chival)
    }
    #   This next part finds the upper limit for the ncp.
    diff <- 1
    while(diff > .00001) {
      if(pchisq(chival,df,uc[2])<llim) 
        uc <- c(uc[1],(uc[1]+uc[2])/2,uc[2]) 
      else uc <- c(uc[2],(uc[2]+uc[3])/2,uc[3])
      diff <- abs(pchisq(chival,df,uc[2]) - llim)
      lcdf <- pchisq(chival,df,uc[2])
    }
    c(uc[2],lcdf)
  }

################

noncentVCI <- function(x, y, conf.level=0.95, dfshift=TRUE){
  if (is.matrix(x) || is.data.frame(x)){
    dat <- as.matrix(x)
  } else {
    dat <- table(x,y)
  }
  chians <- chisq.test(dat, correct=FALSE)
  chival <- chians$statistic
  df = (nrow(dat)-1)*(ncol(dat)-1)
  K = min(nrow(dat), ncol(dat))
  N = sum(dat)
  lo <- lochi(chival, df, 0.95)
  hi <- hichi(chival, df, 0.95)
  if (dfshift){
    ans <- sqrt(c(chival,lo[1]+df,hi[1]+df)/(N*(K-1)))
  } else {
    ans <- sqrt(c(chival,lo[1],hi[1])/(N*(K-1)))
  }
  names(ans) <- c("V", "low", "high")
  ans
}

実際に適用すると次のとおりであり、DescToolsとは一致しない。実は、DescTools はdf を加えずに平方根を取っており、誤りであってかなり偏っている。

noncentVCI(HEC12)

##         V       low      high 
## 0.2790446 0.2345873 0.3258573

CramerV(HEC12, conf.level = 0.95)

##  Cramer V    lwr.ci    upr.ci 
## 0.2790446 0.2235255 0.3179865

noncentVCI(yield2)

##          V        low       high 
## 0.07558305 0.06514759 0.08679269

CramerV(yield2, conf.level = 0.95)

##   Cramer V     lwr.ci     upr.ci 
## 0.07558305 0.05255717 0.07778829

なお、K.Y.氏の blog にでていた計算で、一部のみ変更すれば、正しく同じ計算結果になります。 http://langstat.hatenablog.com/entry/20150103/1420256759

vci <- function(data, alpha=0.05)
{
  require(vcd)
  require(MBESS)
  data <- as.matrix(data)
  V <- assocstats(data)
  cramer <- V$cramer
  K <- min(nrow(data), ncol(data))
#  cat("V =", cramer,"\n")
  nc.chisq  <-   # non-centrality parameter
    conf.limits.nc.chisq(Chi.Square  =  V$chisq_tests[2, 1],
                         df = V$chisq_tests[2, 2],
                         conf.level = 1 - alpha)
  vci.lower <- sqrt((V$chisq_tests[2,2] +
                       as.numeric(nc.chisq[1])) /
                      ((K - 1) * sum(data)))　　　　　　　####
  vci.upper <- sqrt((V$chisq_tests[2,2] + 
                       as.numeric(nc.chisq[3])) /
                      ((K - 1) * sum(data)))              ####
  Vci <- c(cramer, vci.lower,vci.upper)
  names(Vci) <- c("V","ciV L", "ciV U")
  print(Vci)
  res <- list(Vci=Vci, V=V, nc=nc.chisq)
  if (ncol(data)==2 && nrow(data)==2){
    odds <- oddsratio(data, log=FALSE)
    oci <- confint(oddsratio(dat1, log = F))
    res <- list(res, odds=odds, oci=oci)
  }
  return(invisible(res))
}

vci(HEC12)

##         V     ciV L     ciV U 
## 0.2790446 0.2345888 0.3258589

vci(yield2)

##          V      ciV L      ciV U 
## 0.07558305 0.06514742 0.08679186

Bootstrap

RVAideMemoire::cramer.test　のスクリプトを利用し、bootstrap サンプルや値を取得できるように cramer2.test を設計した。ただし、knitr でエラーするので、その結果のみ提示する。黒破線がbootstrap , 緑戦が DescTools, 青線が noncentrality, 赤線が Fisher’s Z. なお、グラフをここに提示しようとしたが、不都合にて中止した。

CI CV

T.Yano

2015年9月28日

Confidence Interval of Cramer’s V

Fisher’s Z transformation

Noncentral chi-sq distribution

Bootstrap