Skip to content

Instantly share code, notes, and snippets.

@ysaito8015
Created May 11, 2014 01:19
Show Gist options
  • Save ysaito8015/1279a91859d3afd40a0b to your computer and use it in GitHub Desktop.
Save ysaito8015/1279a91859d3afd40a0b to your computer and use it in GitHub Desktop.
オッズレシオ計算スクリプト
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)
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment