オッズレシオ計算スクリプト

calcOddsRatio.R

calcOddsRatio <- function(mymatrix,alpha=0.05,referencerow=2,quiet=FALSE)
{
  numrow <- nrow(mymatrix)
    myrownames <- rownames(mymatrix)
    for (i in 1:numrow)
    {
      rowname <- myrownames[i]
        DiseaseUnexposed <- mymatrix[referencerow,1]
        ControlUnexposed <- mymatrix[referencerow,2]
        if (i != referencerow)
        {
          DiseaseExposed <- mymatrix[i,1]
            ControlExposed <- mymatrix[i,2]
            totExposed <- DiseaseExposed + ControlExposed
            totUnexposed <- DiseaseUnexposed + ControlUnexposed
            probDiseaseGivenExposed <- DiseaseExposed/totExposed
            probDiseaseGivenUnexposed <- DiseaseUnexposed/totUnexposed
            probControlGivenExposed <- ControlExposed/totExposed
            probControlGivenUnexposed <- ControlUnexposed/totUnexposed
# オッズ比の計算
            oddsRatio <- (probDiseaseGivenExposed*probControlGivenUnexposed)/(probControlGivenExposed*probDiseaseGivenUnexposed)
            if (quiet == FALSE)
            {
              print(paste("category =", rowname, ", odds ratio = ",oddsRatio))
            }
# 信頼区間の計算
          confidenceLevel <- (1 - alpha)*100
            sigma <- sqrt((1/DiseaseExposed)+(1/ControlExposed)+(1/DiseaseUnexposed)+(1/ControlUnexposed)) # sigmaは,オッズ比の対数の推定値の標準誤差
            z <- qnorm(1-(alpha/2))
            lowervalue <- oddsRatio * exp(-z * sigma)
            uppervalue <- oddsRatio * exp( z * sigma)
            if (quiet == FALSE)
            {
              print(paste("category =", rowname, ", ", confidenceLevel, "% confidence interval = [",lowervalue,",",uppervalue,"]"))
            }
        }
    }
  if (quiet == TRUE && numrow == 2) # 処理水準が 2 つのみの場合 (曝露/非曝露)
  {
    return(oddsRatio)
  }
}