Power estimations in Mitochondrial DNA association studies: a Samuels et al. impementation

server.R

#under GNU General Public License
#server.R
library(shiny)
library(ggplot2)
 
shinyServer(function(input, output) {
## My function:
f <- function(n.cases, p0, OR.cas.ctrl, Nh, sig.level) {
ncases <- n.cases
p0 <- p0
Nh <- Nh
OR.cas.ctrl <- OR.cas.ctrl
sig.level <- sig.level
# Parameters related to sig.level, from [Table 2] of Samuels et al.
# For 90% power and alpha = .05, Nscaled = 8.5
if (sig.level == 0.05){
A <- -28 # Parameter A for alpha=.05
x0 <- 2.6 # Parameter x0 for alpha=.05
d <- 2.4 # Parameter d for alpha=.05
}
if (sig.level == 0.01){
A <- -13 # Parameter A for alpha=.01
x0 <- 5 # Parameter x0 for alpha=.01
d <- 2.5 # Parameter d for alpha=.01
}
if (sig.level == 0.001){
A <- -7 # Parameter A for alpha=.001
x0 <- 7.4 # Parameter x0 for alpha=.001
d <- 2.8 # Parameter d for alpha=.001
}
out.pow <- NULL # initialize vector
for(n.cases in n.cases){
OR.ctrl.cas <- 1 / OR.cas.ctrl # 1. CALCULATE P1 FROM A PREDEFINED P0, AND A DESIRED OR
OR <- OR.ctrl.cas
bracket.pw <- p0 / (OR - OR*p0) # obtained after isolating p1 in OR equation [3].
p1 <- bracket.pw / (1 + bracket.pw)
Nh037 <- Nh^0.37 # 2. CALCULATE NSCALED
num.n <- ncases*((p1-p0)^2)
den.n <- (p1*(1-p1) + p0*(1-p0))*Nh037
Nscaled <- num.n/den.n
num.power <- A - 100 # 3. CALCULATE POWER
den.power <- 1 + exp((Nscaled - x0)/d)
power <- 100 + (num.power/den.power) # The power I have to detect a given OR with my data, at a given alpha
}
OR <- OR.cas.ctrl
out.pow <- data.frame(ncases, Nh, Nscaled, p0, OR, sig.level, power)
out.pow
}
 

# Show the final calculated value
output$powTable <- renderTable(
{mydata <- f(input$n.cases,
input$p0,
input$OR.cas.ctrl,
input$Nh,
input$sig.level)
mydata}
)

})

ui.R

#under GNU General Public License
#ui.R
library(shiny)
library(ggplot2)
 
# Define UI for slider demo application
shinyUI(pageWithSidebar(
# Application title
headerPanel("mtDNA power estimation"),
# Sidebar with sliders that demonstrate various available options
sidebarPanel(
# Simple integer interval
sliderInput("n.cases", "number of cases (ncases):",
min=50, max=1000, value = c(200,400)),
sliderInput("p0", "Frequency of the haplogroup in controls (p0):",
min=0, max=1, value=0.5),
sliderInput("OR.cas.ctrl", "Odds ratio to detect (OR):",
min=1.25, max=4, value=2.3, step=0.05),
sliderInput("Nh", "Number of haplogroup categories (Nh):",
min=5, max=16, value=11),
numericInput("sig.level", "Signification level (only 0.05, 0.01, 0.001):", 0.05),
helpText("This calculator gives power and sample size estimates. It is intended for
genetic association studies on mitochondrial DNA haplogroups, involving Chi-square tests.
The number of controls must be equal to ncases."),
 
helpText("Made from formulae in Samuels et al., AJHG 2006, PMCID: PMC1424681.
Shiny App written by Mareviv.")
),
# Show a table summarizing the values entered
mainPanel(
uiOutput("powTable"))
)
)
 
# To execute this script from inside R, you need to have the ui.R and server.R files into the session's working directory. Then, type:
# runApp()
# To execute directly this Gist, from the Internet to your browser, type:
# shiny:: runGist('5900303')