Rで学ぶ統計講座(初級)

fisher.test.sample <- function(crossTable,alternative=c("two.sided","greater","less")){
    alternative <- match.arg(alternative)
    N <- sum(crossTable)
    if((nrow(crossTable) - 1)*(ncol(crossTable) - 1)==1){
        #観測値の生起確率
        prob <- choose(apply(crossTable,1,sum)[[1]],crossTable[1])*choose(apply(crossTable,1,sum)[[2]],crossTable[2])/choose(N,apply(crossTable,2,sum)[[1]])
        if(is.nan(prob))return(NA) #度数が大きすぎて計算できない時は終了

        prob.sum <- 0 #この変数に合計確率を計算していく
        prob.tmp <- 1 #試行ごとの生起確率、これを観測値の生起確率と比較する
        x <- 0
        if(alternative != "less"){
            #1行1列目の観測値を一つずつ減らしていき、生起確率を計算する
            while(TRUE){
                x <- x + 1
                prob.tmp <- choose(apply(crossTable,1,sum)[[1]],crossTable[1] - x)*choose(apply(crossTable,1,sum)[[2]],crossTable[2] + x)/choose(N,apply(crossTable,2,sum)[[1]])
                if(crossTable[1] - x < 0 || prob.tmp == 0) break
                if(alternative != "greater"){
                    #実測値の確率を下回った時点で終了
                    if(prob.tmp < prob) break
                }
                #実測値を上回る生起確率の時は合計確率に加えていく
                prob.sum <- prob.sum + prob.tmp
            }
        }
        x <- 0
        if(alternative != "greater"){
            #1行1列目の観測値を一つずつ増やしていき、生起確率を計算する
            while(TRUE){
                x <- x + 1
                prob.tmp <- choose(apply(crossTable,1,sum)[[1]],crossTable[1] + x)*choose(apply(crossTable,1,sum)[[2]],crossTable[2] - x)/choose(N,apply(crossTable,2,sum)[[1]])
                if(crossTable[2] - x < 0 || prob.tmp == 0) break
                if(alternative != "less"){
                    if(prob.tmp < prob) break
                }
                prob.sum <- prob.sum + prob.tmp
            }
        }
        p.value <- 1 - prob.sum #1から合計確率を引いた値がp値
    }else return(NA)
    return (p.value)
}

#fisher.test.sample関数の読み込み
source("http://kyoto-edu.sakura.ne.jp/weblesson/statistics/fisher.test.R", encoding="UTF-8")
#クロス集計をする元データの読み込み
data <- read.csv("http://kyoto-edu.sakura.ne.jp/weblesson/statistics/data/gender-opinion.csv", fileEncoding = "utf-8")
gender <- data$性別
opinion <- data$Aという意見に
#クロス集計表作成
crossTable <- table(gender,opinion)
#作成したクロス表のFisherの正確確率を計算する
fisher.test.sample(crossTable) #両側
fisher.test.sample(crossTable,alternative="less") #片側(<)
fisher.test.sample(crossTable,alternative="greater") #片側(>)