Permalink
Browse files

added covariate adjustment to 2propD

  • Loading branch information...
jm3594 committed Jan 5, 2018
1 parent 715fe41 commit 553ada423656de3c7adc8cf5d8dcddb42499b6af
Showing with 93 additions and 17 deletions.
  1. +40 −8 R/crtpwr.2propD.R
  2. +31 −4 inst/shiny-examples/app/app.R
  3. +2 −2 inst/shiny-examples/app/helpers.R
  4. +7 −1 inst/shiny-examples/app/labels.R
  5. +13 −2 man/crtpwr.2propD.Rd
View
@@ -28,6 +28,14 @@
#' subject level. This should be used for a cohort design or a mixture of cohort
#' and cross-sectional designs. In a purely cross-sectional design (baseline subjects
#' are completely different from post-test subjects), this value should be 0.
#' @param covdf The degrees of freedom used by group-level covariates. A value of 0 means no regression
#' adjustment using covariates.
#' @param pvar_c The expected cluster-level proportion of variance to be explained by regression
#' adjustment for covariates. If covariate adjustment isn't used or its effects not known, this value should
#' be 0.
#' @param pvar_s The expected subject-level proportion of variance to be explained by regression
#' adjustment for covariates. If covariate adjustment isn't used or its effects not known, this value should
#' be 0.
#' @param tol Numerical tolerance used in root finding. The default provides
#' at least four significant digits.
#' @return The computed argument.
@@ -49,6 +57,7 @@ crtpwr.2propD <- function(alpha = 0.05, power = 0.80,
m = NA, n = NA,
p = NA, d = NA, icc = NA,
rho_c = NA, rho_s = NA,
covdf = 0, pvar_c = 0, pvar_s = 0,
tol = .Machine$double.eps^0.25){
if(!is.na(m) && m <= 1) {
@@ -69,12 +78,12 @@ crtpwr.2propD <- function(alpha = 0.05, power = 0.80,
}
}
needlist <- list(alpha, power, m, n, p, d, icc, rho_c, rho_s)
neednames <- c("alpha", "power", "m", "n", "p", "d", "icc", "rho_c", "rho_s")
needlist <- list(alpha, power, m, n, p, d, icc, rho_c, rho_s, covdf, pvar_s, pvar_c)
neednames <- c("alpha", "power", "m", "n", "p", "d", "icc", "rho_c", "rho_s", "covdf", "pvar_s", "pvar_c")
needind <- which(unlist(lapply(needlist, is.na))) # find NA index
if (length(needind) != 1) {
stop("Exactly one of 'alpha', 'power', 'm', 'n', 'p', 'd', 'icc', 'rho_c', or 'rho_s' must be NA.")
stop("Exactly one of 'alpha', 'power', 'm', 'n', 'p', 'd', 'icc', 'rho_c', 'rho_s', 'covdf', 'pvar_s', or 'pvar_c' must be NA.")
}
target <- neednames[needind]
@@ -87,9 +96,11 @@ crtpwr.2propD <- function(alpha = 0.05, power = 0.80,
# between cluster: varb = p*(1-p)*icc
# within cluster: varw = p*(1-p)*(1 - icc)
# 2*2*(p*(1-p)*(1 - icc)*(1 - rho_s) + n*p*(1-p)*icc*(1 - rho_c))/(n*m)
ncp <- d/sqrt(2*2*(p*(1-p)*(1 - icc)*(1 - rho_s) + n*p*(1-p)*icc*(1 - rho_c))/(n*m))
varb <- p*(1-p)*icc
varw <- p*(1-p)*(1 - icc)
ncp <- d/sqrt(2*2*(varw*(1 - rho_s)*(1 - pvar_s) + n*varb*(1 - rho_c)*(1 - pvar_c))/(n*m))
pt(tcrit, 2*(m - 1), ncp, lower.tail = FALSE)
pt(tcrit, 2*(m - 1) - covdf, ncp, lower.tail = FALSE)
})
# calculate alpha
@@ -107,8 +118,8 @@ crtpwr.2propD <- function(alpha = 0.05, power = 0.80,
# calculate m
if (is.na(m)) {
m <- stats::uniroot(function(m) eval(pwr) - power,
interval = c(2 + 1e-10, 1e+07),
tol = tol)$root
interval = c(2 + covdf + 1e-10, 1e+07),
tol = tol, extendInt = "upX")$root
}
# calculate n
@@ -142,17 +153,38 @@ crtpwr.2propD <- function(alpha = 0.05, power = 0.80,
# calculate rho_c
if (is.na(rho_c)){
rho_c <- stats::uniroot(function(rho_c) eval(pwr) - power,
interval = c(1e-07, 1e+07),
interval = c(1e-07, 1 - 1e-7),
tol = tol, extendInt = "upX")$root
}
# calculate rho_s
if (is.na(rho_s)){
rho_s <- stats::uniroot(function(rho_s) eval(pwr) - power,
interval = c(1e-07, 1 - 1e-7),
tol = tol, extendInt = "upX")$root
}
# calculate covdf
if (is.na(covdf)){
covdf <- stats::uniroot(function(covdf) eval(pwr) - power,
interval = c(1e-07, 1e+07),
tol = tol, extendInt = "upX")$root
}
# calculate pvar_c
if (is.na(pvar_c)){
pvar_c <- stats::uniroot(function(pvar_c) eval(pwr) - power,
interval = c(1e-07, 1 - 1e-7),
tol = tol, extendInt = "upX")$root
}
# calculate pvar_s
if (is.na(pvar_s)){
pvar_s <- stats::uniroot(function(pvar_s) eval(pwr) - power,
interval = c(1e-07, 1 - 1e-7),
tol = tol, extendInt = "upX")$root
}
structure(get(target), names = target)
}
@@ -14,7 +14,7 @@ names2mean <- c("alpha","power","m","n","cv","d","icc","varw","method")
names2meanD <- c("alpha","power","m","n","d","icc","rho_c","rho_s","varw")
names2meanM <- c("alpha","power","m","n","d","icc","varw","rho_m")
names2prop <- c("alpha","power","m","n","cv","p1","p2","icc","pooled","p1inc")
names2propD <- c("alpha","power","m","n","p","d","icc","rho_c","rho_s")
names2propD <- c("alpha","power","m","n","p","d","icc","rho_c","rho_s","covdf","pvar_c","pvar_s")
names2propM <- c("alpha","power","m","n","p1","p2","cvm","p1inc")
names2rate <- c("alpha","power","m","py","r1","r2","cvb","r1inc")
@@ -329,6 +329,21 @@ ui <- fluidPage(
bsTooltip("rho_s2propD", rho_stooltip,
'right', options = list(container = "body")),
#----------------------------------------------------------
fluidRow(textInput("covdf2propD", covdftext,
value = "0", width = "100%")),
bsTooltip("covdf2propD", covdftooltip,
'right', options = list(container = "body")),
#----------------------------------------------------------
fluidRow(textInput("pvar_c2propD", pvar_ctext,
value = "0", width = "100%")),
bsTooltip("pvar_c2propD", pvar_ctooltip,
'right', options = list(container = "body")),
#----------------------------------------------------------
fluidRow(textInput("pvar_s2propD", pvar_stext,
value = "0", width = "100%")),
bsTooltip("pvar_s2propD", pvar_stooltip,
'right', options = list(container = "body")),
#----------------------------------------------------------
fluidRow(
column(6, style='padding:0px;', actionButton("default2propD", defaulttext, width = "100%")),
column(6, style='padding:0px;', actionButton("clear2propD", clearalltext, width = "100%"))
@@ -1028,6 +1043,9 @@ server <- function(input, output, session){
updateTextInput(session, inputId = "p2propD", value = "")
updateTextInput(session, inputId = "d2propD", value = "")
updateTextInput(session, inputId = "icc2propD", value = "")
updateTextInput(session, inputId = "covdf2propD", value = "0")
updateTextInput(session, inputId = "pvar_c2propD", value = "0")
updateTextInput(session, inputId = "pvar_s2propD", value = "0")
}
) # end observeEvent(input$default2propD ...
@@ -1045,6 +1063,9 @@ server <- function(input, output, session){
updateTextInput(session, inputId = "p2propD", value = "")
updateTextInput(session, inputId = "d2propD", value = "")
updateTextInput(session, inputId = "icc2propD", value = "")
updateTextInput(session, inputId = "covdf2propD", value = "")
updateTextInput(session, inputId = "pvar_c2propD", value = "")
updateTextInput(session, inputId = "pvar_s2propD", value = "")
}
) # end observeEvent(input$clear2propD ...
@@ -1062,6 +1083,9 @@ server <- function(input, output, session){
p <- make_sequence(isolate(input$p2propD))
d <- make_sequence(isolate(input$d2propD))
icc <- make_sequence(isolate(input$icc2propD))
covdf <- make_sequence(isolate(input$covdf2propD))
pvar_c <- make_sequence(isolate(input$pvar_c2propD))
pvar_s <- make_sequence(isolate(input$pvar_s2propD))
if(!is.na(power)){
validate(
@@ -1088,14 +1112,17 @@ server <- function(input, output, session){
icc,
rho_c,
rho_s,
covdf,
pvar_c,
pvar_s,
stringsAsFactors = FALSE)
# record the column index of the target parameter
needind <- which(is.na(tab[1,]))
# validate that only one input is blank
validate(
need(length(needind) == 1,
"Exactly one of 'alpha', 'power', 'p', 'd', 'm', 'n', 'icc', 'rho_c', and 'rho_s' must be left blank."
"Exactly one of 'alpha', 'power', 'p', 'd', 'm', 'n', 'icc', 'rho_c', 'rho_s', 'covdf', 'pvar_c', and 'pvar_s' must be left blank."
)
)
names(tab) <- names2propD
@@ -1120,7 +1147,7 @@ server <- function(input, output, session){
# create 2propD output table
output$table2propD <- DT::renderDataTable(
res2propD()[, c(1:10)],
res2propD()[, c(1:13)],
server = FALSE,
extensions = 'Buttons',
filter = 'top',
@@ -1130,7 +1157,7 @@ server <- function(input, output, session){
buttons = list(list(extend = 'csv', filename = paste('data-2propD-', Sys.time(), sep=''), text = 'Download')),
autoWidth = TRUE,
columnDefs = list(list(className = 'dt-center', targets = '_all'),
list(width = '500px', targets = 10)),
list(width = '500px', targets = 13)),
pageLength = 10,
lengthMenu = list(c(10, 25, 100, -1), list('10', '25', '100', 'All'))
)
@@ -58,11 +58,11 @@ crtpwr.2prop.safe <- function(alpha,power,m,n,cv,p1,p2,icc,pooled,p1inc){
res
}
crtpwr.2propD.safe <- function(alpha,power,m,n,p,d,icc,rho_c,rho_s){
crtpwr.2propD.safe <- function(alpha,power,m,n,p,d,icc,rho_c,rho_s,covdf,pvar_c,pvar_s){
# make safe version
fun <- safely(crtpwr.2propD, otherwise = NA)
# store result
res <- fun(alpha,power,m,n,p,d,icc,rho_c,rho_s)
res <- fun(alpha,power,m,n,p,d,icc,rho_c,rho_s,covdf,pvar_c,pvar_s)
# if res$error NULL, set to NA, otherwise set to message
if(is.null(res$error)){
res$error = 'None'
@@ -68,14 +68,20 @@ pooledtooltip <- 'Select to indicate if pooled variance is desired.'
ptext <- "Expected mean proportion (p)"
ptooltip <- "The expected mean proportion at the post-test, averaged across treatment and control arms."
covdftext <- "Covariate degrees of freedom (covdf)"
covdftooltip <- "The degrees of freedom used by group-level covariates. A value of 0 means no regression adjustment using covariates."
pvar_ctext <- "Cluster-level proportion of variance explained (pvar_c)"
pvar_ctooltip <- "The expected cluster-level proportion of variance to be explained by regression adjustment for covariates."
pvar_stext <- "Subeject-level proportion of variance explained (pvar_s)"
pvar_stooltip <- "The expected subject-level proportion of variance to be explained by regression adjustment for covariates."
# -----------------------------------------------------------------------------
# specific to 2propM:
cvmtext <- "Within-pair Outcome CV (cvm)"
cvmtooltip <- "The coefficient of variation in the outcome within matched clusters."
# specific to 2rate:
r1text <- 'Rate 1 (r1)'
r1tooltip <- 'The expected rate in the treatment group.'
View

Some generated files are not rendered by default. Learn more.

Oops, something went wrong.

0 comments on commit 553ada4

Please sign in to comment.