Skip to content

Instantly share code, notes, and snippets.

@aurora-mareviv

aurora-mareviv/server.R

Last active August 29, 2015 14:14
Show Gist options
  • Save aurora-mareviv/84461fe8abe81ab7f78e to your computer and use it in GitHub Desktop.
Save aurora-mareviv/84461fe8abe81ab7f78e to your computer and use it in GitHub Desktop.
Download ZIP
Chi-square power / sample estimations with R!
Raw
server.R
#under GNU General Public License
#server.R
library(shiny)
library(pwr)
shinyServer(function(input, output) {
## My pwr function:
f <- function(eff.size, n.observs, degfree, signif) {
require(pwr)
signif <- as.numeric(signif)
dat <- pwr.chisq.test(w = eff.size, N = n.observs, df = degfree, sig.level = signif)
dati <- data.frame(dat$w, dat$N, dat$df, dat$sig.level, dat$power)
colnames(dati) <- c("effect.size", "cases", "degrees.of.freedom", "significance.level", "power")
dati
}
## My function to calc w from two proportions:
g <- function(diff.percent.control, diff.percent.case, n.observs) {
ntot <- n.observs
nrow <- ntot/2
# Observed ROW proportions:
p1r.control1 <- diff.percent.control
p1r.case1 <- diff.percent.case
p1r.control2 <- 1- p1r.control1
p1r.case2 <- 1- p1r.case1
.TableOPr <- matrix(c(p1r.control1,p1r.control2,p1r.case1,p1r.case2), 2, 2, byrow=TRUE)
rownames(.TableOPr) <- c('ctrl', 'case')
colnames(.TableOPr) <- c('1', '2')
.TableOPr
# Observed counts:
c1.control1 <- nrow*p1r.control1
c1.case1 <- nrow*p1r.case1
c1.control2 <- nrow-c1.control1
c1.case2 <- nrow-c1.case1
.TableOC <- matrix(c(c1.control1,c1.control2,c1.case1,c1.case2), 2, 2, byrow=TRUE)
rownames(.TableOC) <- c('ctrl', 'case')
colnames(.TableOC) <- c('1', '2')
.TableOC
# Observed table proportions:
p1.control1 <- c1.control1*1/ntot
p1.case1 <- c1.case1*1/ntot
p1.control2 <- c1.control2*1/ntot
p1.case2 <- c1.case2*1/ntot
.TableOP <- matrix(c(p1.control1,p1.control2,p1.case1,p1.case2), 2, 2, byrow=TRUE)
rownames(.TableOP) <- c('ctrl', 'case')
colnames(.TableOP) <- c('1', '2')
.TableOP
# Expected counts:
c0.control1 <- (c1.control1+c1.case1)*(c1.control1+c1.control2)/ntot
c0.control2 <- (c1.control2+c1.case2)*(c1.control1+c1.control2)/ntot
c0.case1 <- (c1.control1+c1.case1)*(c1.case1+c1.case2)/ntot
c0.case2 <- (c1.control2+c1.case2)*(c1.case1+c1.case2)/ntot
.TableEC <- matrix(c(c0.control1,c0.control2,c0.case1,c0.case2), 2, 2, byrow=TRUE)
rownames(.TableEC) <- c('ctrl', 'case')
colnames(.TableEC) <- c('1', '2')
.TableEC
# Expected proportions:
p0.control1 <- c0.control1*1/ntot
p0.control2 <- c0.control2*1/ntot
p0.case1 <- c0.case1*1/ntot
p0.case2 <- c0.case2*1/ntot
.TableEP <- matrix(c(p0.control1,p0.control2,p0.case1,p0.case2), 2, 2, byrow=TRUE)
rownames(.TableEP) <- c('ctrl', 'case')
colnames(.TableEP) <- c('1', '2')
.TableEP
# calculate recommended Cohen's w:
w <- (((p0.control1 - p1.control1)^2/p0.control1) + ((p0.control2 - p1.control2)^2/p0.control2) + ((p0.case1 - p1.case1)^2/p0.case1) + ((p0.case2 - p1.case2)^2/p0.case2))^(1/2)
# Chi <- w^2*ntot
# X <- ((c1.control1 - c0.control1)^2/c0.control1) + ((c1.control2 - c0.control2)^2/c0.control2) + ((c1.case1 - c0.case1)^2/c0.case1) + ((c1.case2 - c0.case2)^2/c0.case2)
w
}
## My sample size McNemar calculator
# in my case, I am assuming that the 0 to 1 change will be 0.20 and that no one will change from 1 to 0
samsize.mcnemar <- function(pi.0.to.1, pi.1.to.0, alpha, beta, sided)
{
pi.01 <- pi.0.to.1
pi.10 <- pi.1.to.0
alpha <- as.numeric(alpha)
sided <- as.numeric(sided)
pi.d <- (pi.01 + pi.10)
N <- (qnorm(1 - alpha/sided) * sqrt(pi.d) + qnorm(1 - beta) *
sqrt(pi.d - (pi.01 - pi.10)^2))^2/(pi.01 - pi.10)^2
sample.size <- ceiling(N)
size.dat <- data.frame(pi.01, pi.10, alpha, beta, sided, sample.size)
size.dat
}
mydata <- reactive(f(input$eff.size,
input$n.observs,
input$degfree,
input$signif))
mydata2 <- reactive(g(input$diff.percent.control,
input$diff.percent.case,
input$n.observs))
mydatanemar <- reactive(samsize.mcnemar(input$pi.0.to.1,
input$pi.1.to.0,
input$alpha,
input$beta,
input$sided))
# Show the final calculated values from Pearson
output$powTable <- renderTable(
{mydata <- f(input$eff.size,
input$n.observs,
input$degfree,
input$signif)
mydata}
)
# Show the final calculated values from McNemar
output$powTableNemar <- renderTable(
{mydatanemar <- samsize.mcnemar(input$pi.0.to.1,
input$pi.1.to.0,
input$alpha,
input$beta,
input$sided)
mydatanemar}
)
output$text1 <- renderText({
paste("We need ", input$n.observs, " patients", " to detect an Odds-Ratio of approximately ", round(input$eff.size/0.065, 2),
" (effect size w = ", input$eff.size, ") [1]",
", using a Chi-square test of ", input$degfree, " degrees of freedom",
" with a ", input$signif, " significance and a ", round(mydata()$power*100, 1), "% minimum power.",
sep="")
})
output$text2 <- renderText({
paste("For a difference in proportions between ", input$diff.percent.control, " in controls", " and ", input$diff.percent.case,
" in cases (in a 2x2 table), the calculated Cohen's w = ", round(mydata2(), 2), ".",
sep="")
})
output$textNemar <- renderText({
paste("We expect that: if the ", round(input$pi.0.to.1*100, 1), "% of the individuals change from not having the condition to having it, and ",
round(input$pi.1.to.0*100, 1), "% of the sample changes to not having the condition, the estimated sample size is ", mydatanemar()$sample.size,
" to detect a significant effect with an alpha error of ", input$alpha, " and a power of ", round((1-input$beta)*100, 1), "%.",
sep="")
})
# Show the Download handler:
output$downloadData <- downloadHandler(
filename = function() { paste(input$mydata, 'chisq_power_table.csv', sep='') },
content = function(file) {
write.csv(mydata(), file, na="")
}
)
})
Raw
ui.R
#under GNU General Public License
# ui.R
library(markdown)
shinyUI(navbarPage("Chi-square power / sample estimations with R!",
tabPanel("Pearson",
sidebarLayout(
sidebarPanel(
# Simple integer interval
sliderInput("degfree", "Degrees of freedom:",
min=1, max=10, value=1),
sliderInput("eff.size", "Effect size to detect:",
min=0.1, max=1.4, value = 0.28, step=0.01),
sliderInput("n.observs", "Total number of observations:",
min=25, max=1000, value = 100, step=5),
selectInput("signif", "Significance level:",
choices = c("0.05", "0.01", "0.001")),
conditionalPanel(
condition = "input.degfree == '1'",
sliderInput("diff.percent.control", "Observed proportion of controls without the condition (for df = 1):",
min=0, max=1, value=0.95),
sliderInput("diff.percent.case", "Observed proportion of cases without the condition (for df = 1):",
min=0, max=1, value=0.75)
),
helpText("The Pearson's Chi-squared test is a statistical test used on independent samples, applied to contingency tables."),
helpText("Effect size (w): the overall difference among groups with respect to a factor: 0.1 is a low diference; 0.25 is a mild effect or difference; 0.5 is a high difference among groups."),
helpText("Degrees of freedom (df): if the number of rows = r and the number of columns = c, then df = (r-1)*(c-1). Typically, df = 1 in a 2x2 table, and df = 2 in a 2-row, 3-column table."),
helpText("Significance level: the p-value threshold of significance."),
helpText("Bibliography:"),
helpText("1. Jacob Cohen. Statistical Power Analysis for the Behavioral Sciences (2nd Edition). Psychology Press, 1988."),
helpText("2. pwr: Basic functions for power analysis. Stephane Champely, 2012. R package version 1.1.1."),
helpText("Written in R/Shiny by A. Baluja."),
downloadButton('downloadData', 'Download dataset')
),
mainPanel(
uiOutput("powTable"),
textOutput("text1"),
conditionalPanel(
condition = "input.degfree == '1'",
textOutput("text2")
)
)
)
),
tabPanel("Mc Nemar",
sidebarLayout(
sidebarPanel(
# Simple integer interval
sliderInput("pi.0.to.1", "Proportion of individuals that change from factor=no to factor=yes:",
min=0, max=1, value=0.25, step=0.05),
sliderInput("pi.1.to.0", "Proportion of individuals that change from factor=yes to factor=no:",
min=0, max=1, value=0, step=0.05),
sliderInput("beta", "Error beta (power=1-beta):",
min=0, max=1, value=0.2, step=0.05),
selectInput("alpha", "Significance level:",
choices = c("0.05", "0.01", "0.001")),
selectInput("sided", "One/ two-sided test:",
choices = c("1", "2")),
helpText("McNemar's test is a statistical test used on paired nominal data, applied to 2 × 2 contingency tables."),
helpText("Bibliography:"),
helpText("1. Siegel S and Castellan Jr. N J (1988) Nonparametric Statistics for the Behavioral Sciences 2nd. Ed. McGraw Hill, Inc. Boston Massachusetts. ISBN 0-07-057357-3 p 75."),
helpText("2. Machin D, Campbell M, Fayers, P, Pinol A (1997) Sample Size Tables for Clinical Studies. Second Ed. Blackwell Science IBSN 0-86542-870-0 p. 70-71."),
helpText("3. Code: http://stats.stackexchange.com/questions/87829/sample-size-function-for-mcnemar-test-of-repeated-proportions-in-r."),
helpText("Written in R/Shiny by A. Baluja.")
),
mainPanel(
mainPanel(
uiOutput("powTableNemar"),
textOutput("textNemar")
)
)
)
)))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment