pump_demo.R
In PUMP: Power Under Multiplicity Project

## ----initialize, include = FALSE----------------------------------------------
library( PUMP )

# This allows us to stash precalculated results so the vignette
# compiles in a timely manner for CRAN.
vignette_cache <- here::here("vignettes/output", "MDEScalc.RDS")
recompile <- !file.exists(vignette_cache)
# recompile <- TRUE

knitr::opts_chunk$set(
  cache = FALSE,
  collapse = TRUE,
  comment = "#>",
  warning = FALSE,
  message = FALSE,
  fig.width = 7,
  fig.height = 4,
  fig.align = "center"
)
options(knitr.kable.NA = '')
library( dplyr )
library( ggplot2 )
library( ggthemes )
library( here )
library( kableExtra )
library( knitr )
set.seed( 0440044 )

## ----echo = FALSE-------------------------------------------------------------
info <- pump_info(comment = FALSE)
kable( info$Context, format = 'latex', booktabs = TRUE ) %>%
  kable_styling(position = "center")

## ----echo = FALSE-------------------------------------------------------------
kable(info$Parameters, format = 'latex', booktabs = TRUE ) %>%
  kable_styling(position = "center")

## ----MDEScalc, eval = FALSE---------------------------------------------------
# m <- pump_mdes(
#   d_m = "d3.2_m3fc2rc",         # choice of design and analysis strategy
#   MTP = "HO",                   # multiple testing procedure
#   target.power = 0.80,          # desired power
#   power.definition = "D1indiv", # power type
#   M = 5,                        # number of outcomes
#   J = 3,                        # number of schools per block
#   K = 21,                       # number districts
#   nbar = 258,                   # average number of students per school
#   Tbar = 0.50,                  # prop treated
#   alpha = 0.05,                 # significance level
#   numCovar.1 = 5,               # number of covariates at level 1
#   numCovar.2 = 3,               # number of covariates at level 2
#   R2.1 = 0.1, R2.2 = 0.7,       # explanatory power of covariates for each level
#   ICC.2 = 0.05, ICC.3 = 0.4,    # intraclass correlation coefficients
#   rho = 0.4 )                   # how correlated outcomes are

## ----MDEScalc-compile, echo = FALSE-------------------------------------------
if (recompile)
{
    m <- pump_mdes(
      d_m = "d3.2_m3fc2rc",         # choice of design and analysis strategy
      MTP = "HO",                   # multiple testing procedure
      target.power = 0.80,          # desired power
      power.definition = "D1indiv", # power type
      M = 5,                        # number of outcomes
      J = 3,                        # number of schools per block
      K = 21,                       # number districts
      nbar = 258,                   # average number of students per school
      Tbar = 0.50,                  # prop treated
      alpha = 0.05,                 # significance level
      numCovar.1 = 5,               # number of covariates at level 1
      numCovar.2 = 3,               # number of covariates at level 2
      R2.1 = 0.1, R2.2 = 0.7,       # explanatory power of covariates for each level
      ICC.2 = 0.05, ICC.3 = 0.4,    # intraclass correlation coefficients
      rho = 0.4 ) 
    dir.create( here::here("vignettes/output") )
    saveRDS(m, here::here("vignettes/output", "MDEScalc.RDS"))
} else {
    m <- readRDS( here::here("vignettes/output", "MDEScalc.RDS") )
}

## ----echo = TRUE--------------------------------------------------------------
knitr::kable( m, digits = 3 ) %>%
  kableExtra::kable_styling( position = "center" )

## ----MDEScalcmin1, echo = TRUE, eval = FALSE----------------------------------
# m2 <- update( m, power.definition = "min1" )

## ----MDEScalcmin1-compile, echo = FALSE---------------------------------------
if (recompile)
{
    m2 <- update( m, power.definition = "min1" )
    saveRDS(m2, here::here("vignettes/output", "MDEScalcmin1.RDS"))
} else
{
    m2 <- readRDS(here::here("vignettes/output", "MDEScalcmin1.RDS"))
}
print( m2 )

## ----MDESwithNumZero, echo = TRUE, eval = FALSE-------------------------------
# m3 <- update( m2, numZero = 2 )

## ----MDESwithNumZero-compile, echo = FALSE------------------------------------
if (recompile)
{
    m3 <- update( m2, numZero = 2 )
    saveRDS(m3, here::here("vignettes/output", "MDESwithNumZero.RDS"))
} else
{
    m3 <- readRDS(here::here("vignettes/output", "MDESwithNumZero.RDS"))
}
print( m3 )

## ----samplesizecalc, eval = FALSE---------------------------------------------
# smp <- pump_sample(
#   d_m = "d3.2_m3fc2rc",
#   MTP = "HO",
#   typesample = "K",
#   target.power = 0.80, power.definition = "min1", tol = 0.01,
#   MDES = 0.10, M = 5, nbar = 258, J = 3,
#   Tbar = 0.50, alpha = 0.05, numCovar.1 = 5, numCovar.2 = 3,
#   R2.1 = 0.1, R2.2 = 0.7, ICC.2 = 0.05, ICC.3 = 0.40, rho = 0.4 )
# 
# print( smp )

## ----samplesizecalc-compile, echo = FALSE-------------------------------------
if (recompile)
{
    smp <- pump_sample(
      d_m = "d3.2_m3fc2rc",
      MTP = "HO",
      typesample = "K",
      target.power = 0.80, power.definition = "min1", tol = 0.01,
      MDES = 0.10, M = 5, nbar = 258, J = 3,
      Tbar = 0.50, alpha = 0.05, numCovar.1 = 5, numCovar.2 = 3,
      R2.1 = 0.1, R2.2 = 0.7, ICC.2 = 0.05, ICC.3 = 0.40, rho = 0.4 )
    saveRDS(smp, here::here("vignettes/output", "samplesizecalc.RDS"))
} else
{
    smp <- readRDS(here::here("vignettes/output", "samplesizecalc.RDS"))
}
print( smp )

## ----samplesizeverify, eval = FALSE-------------------------------------------
# p_check <- update( smp, type = "power", tnum = 20000,
#                    long.table = TRUE )

## ----samplesizeverify-compile, echo = FALSE-----------------------------------
if (recompile)
{
    p_check <- update( smp, type = "power", tnum = 20000,
                   long.table = TRUE )
    saveRDS(p_check, here::here("vignettes/output", "samplesizeverify.RDS"))
} else
{
    p_check <- readRDS(here::here("vignettes/output", "samplesizeverify.RDS"))
}
knitr::kable( p_check, digits = 2 ) %>%
  kableExtra::kable_styling(position = "center")

## ----plotsamplepowercurve, fig.height = 3.5, fig.width=5, echo=TRUE-----------
plot( smp )

## ----powbase, eval = FALSE----------------------------------------------------
# pow <- update( p_check, tnum = 10000 )

## ----powbase-compile, echo = FALSE--------------------------------------------
if (recompile)
{
    pow <- update( p_check, tnum = 10000 )
    saveRDS(pow, here::here("vignettes/output", "powbase.RDS"))
} else
{
    pow <- readRDS(here::here("vignettes/output", "powbase.RDS"))
}

## ----othercorrections, eval = FALSE-------------------------------------------
# p2 <- update( pow,
#               MTP = c( "BF", "HO", "WY-SD" ) )
# plot( p2 )

## ----othercorrections-compile, echo = FALSE-----------------------------------
if (recompile)
{
    p2 <- update( pow,
      MTP = c( "BF", "HO", "WY-SD" ) )
    saveRDS(p2, here::here("vignettes/output", "othercorrections.RDS"))
} else
{
    p2 <- readRDS(here::here("vignettes/output", "othercorrections.RDS"))
}
plot( p2 )

## ----powICC, eval = FALSE-----------------------------------------------------
# p_b <- update( pow, ICC.2 = 0.20, ICC.3 = 0.25 )
# print( p_b )

## ----powICC-compile, echo = FALSE---------------------------------------------
if (recompile)
{
    p_b <- update( pow, ICC.2 = 0.20, ICC.3 = 0.25 )
    saveRDS(p_b, here::here("vignettes/output", "powICC.RDS"))
} else
{
    p_b <- readRDS(here::here("vignettes/output", "powICC.RDS"))
}
print( p_b )

## ----powR2, eval = FALSE------------------------------------------------------
# p_d <- update( pow,
#           	   R2.1 = c( 0.1, 0.3, 0.1, 0.2, 0.2 ),
#           	   R2.2 = c( 0.4, 0.8, 0.3, 0.2, 0.2 ) )
# print( p_d )

## ----powR2-compile, echo = FALSE----------------------------------------------
if (recompile)
{
    p_d <- update( pow,
          	   R2.1 = c( 0.1, 0.3, 0.1, 0.2, 0.2 ),
          	   R2.2 = c( 0.4, 0.8, 0.3, 0.2, 0.2 ) )
    saveRDS(p_d, here::here("vignettes/output", "powR2.RDS"))
} else
{
    p_d <- readRDS(here::here("vignettes/output", "powR2.RDS"))
}
print( p_d )

## -----------------------------------------------------------------------------
summary( p_d )

## ----ICCgrid, eval = FALSE----------------------------------------------------
# grid <- update_grid( pow,
#         	ICC.2 = seq( 0, 0.3, 0.05 ),
#         	ICC.3 = seq( 0, 0.60, 0.20 ) )
# 
# plot( grid, power.definition = "min1" )

## ----ICCgrid-compile, echo = FALSE--------------------------------------------
if (recompile)
{
    grid <- update_grid( pow,
        	ICC.2 = seq( 0, 0.3, 0.05 ),
        	ICC.3 = seq( 0, 0.60, 0.20 ) )
    saveRDS(grid, here::here("vignettes/output", "ICCgrid.RDS"))
} else
{
    grid <- readRDS(here::here("vignettes/output", "ICCgrid.RDS"))
}
plot( grid, power.definition = "min1" )

## ----rhogrid, eval = FALSE----------------------------------------------------
# gridRho <- update_grid( pow,
#         	  MTP = c( "BF", "WY-SD" ),
#         	  rho = seq( 0, 0.9, by = 0.15 ),
#         	  tnum = 1000,
#         	  B = 3000 )

## ----rhogrid-compile, echo = FALSE--------------------------------------------
if (recompile)
{
    gridRho <- update_grid( pow,
        	  MTP = c( "BF", "WY-SD" ),
        	  rho = seq( 0, 0.9, by = 0.15 ),
        	  tnum = 1000,
        	  B = 3000 )
    saveRDS(gridRho, here::here("vignettes/output", "rhogrid.RDS"))
} else
{
    gridRho <- readRDS(here::here("vignettes/output", "rhogrid.RDS"))
}

## -----------------------------------------------------------------------------
plot( gridRho )

## ----numzerogrid, fig.height = 2.5, eval = FALSE------------------------------
# gridZero <- update_grid( pow,
#         	             numZero = 0:4,
#                          M = 5 )
# plot( gridZero, nrow = 1 )

## ----numzerogrid-compile, echo = FALSE, fig.height = 2.5----------------------
if (recompile)
{
    gridZero <- update_grid( pow,
                             numZero = 0:4,
                             M = 5 )
    saveRDS(gridZero, here::here("vignettes/output", "numzerogrid.RDS"))
} else
{
   gridZero <- readRDS(here::here("vignettes/output", "numzerogrid.RDS"))
}
plot( gridZero, nrow = 1 )