Skip to content

Instantly share code, notes, and snippets.

@ysaito8015
Created May 11, 2014 01:21
Show Gist options
  • Save ysaito8015/ee8cc2a8a7e09e53fd37 to your computer and use it in GitHub Desktop.
Save ysaito8015/ee8cc2a8a7e09e53fd37 to your computer and use it in GitHub Desktop.
相対危険度スクリプト
calcRelativeRisk <- function(mymatrix,alpha=0.05,referencerow=2)
{
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
# 相対危険度の計算
relativeRisk <- probDiseaseGivenExposed/probDiseaseGivenUnexposed
print(paste("category =", rowname, ", relative risk = ",relativeRisk))
# 信頼区間の計算
confidenceLevel <- (1 - alpha)*100
sigma <- sqrt((1/DiseaseExposed) - (1/totExposed) + (1/DiseaseUnexposed) - (1/totUnexposed)) # sigmaは相対危険度の対数の推定値の標準誤差
z <- qnorm(1-(alpha/2))
lowervalue <- relativeRisk * exp(-z * sigma)
uppervalue <- relativeRisk * exp( z * sigma)
print(paste("category =", rowname, ", ", confidenceLevel,"% 信頼区間 = [",lowervalue,",",uppervalue,"]"))
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment