scratch/test-sampleSizeSignificance.R

In florafauna/RStest: Design and Analysis of Replication Studies

rm(list=ls())
library(ReplicationSuccess)
library(testthat)
library(dplyr)
## sapply(list.files("../R", pattern='\\.R$', full.names = TRUE), source)


## from new-formulas.R
## slighlty adapted, e.g.,
## - shrinkage >= 1 - u*sqrt(H)/zo changes to shrinkage * zo >= zo - u*sqrt(H)
##   to avoid division by 0
sampleSizeSignificance2 <- function(zo,
                                        power = NA,
                                        d = NA,
                                        level = 0.025,
                                        alternative = "one.sided",
                                        designPrior = "conditional",
                                        h = 0,
                                        shrinkage = 0) {

    ## check that the length of all aguments is <= 1
    if(any(length(zo) > 1, length(power) > 1, length(d) > 1, length(level) > 1,
           length(alternative) > 1, length(designPrior) > 1, length(h) > 1, length(shrinkage) >1))
        stop("all arguments have to be of length one")
    
    ## checks that only one of 'power' and 'd' has to be specified
    if (is.na(d) && is.na(power)) stop("either 'power' or 'd' has to be specified")
    if (!is.na(d) && !is.na(power)) stop("only one of 'power' or 'd' has to be specified")

    ## sample size calculation based on power
    if (is.na(d)) {
        ## input checks
        if (!(designPrior %in% c("conditional", "predictive", "EB")))
            stop('designPrior must be either "conditional", "predictive", or "EB"')
        if (!is.numeric(power) || (power <= 0 || power >= 1))
            stop("power must be numeric and in (0, 1)")
        if (!is.numeric(level) || (level <= 0 || level >= 1))
            stop("level must be numeric and in (0, 1)")
        if (!is.numeric(h) || h < 0)
            stop("h must be numeric and cannot be negative")
        if (!is.numeric(shrinkage) || (shrinkage < 0 || shrinkage > 1))
            stop("shrinkage must be numeric and in [0, 1]")
        u <- qnorm(p = power)
        v <- p2z(p = level, alternative = alternative)
        ## for conditional designPrior use analytical solution
        if (designPrior == "conditional") {
            c <- (u + v)^2*(1/((1 - shrinkage)*zo))^2
        }
        if (designPrior %in% c("predictive", "EB")) {
            if (designPrior == "EB") {
                shrinkage <- pmin((1 + h)/zo^2, 1)
                H <- 1 - shrinkage + 2*h - shrinkage*h
            } else {
                H <- 1 + 2*h
            }
            if (shrinkage * zo >= zo - u*sqrt(H)) { ## changed to avoid division by 0 !!!!
                c <- NaN
            } else {
                zos <- (1 - shrinkage)*zo
                num <- zos*v + u*sqrt(zos^2 + H*(v^2 - u^2))
                denom <- zos^2 - u^2*H
                sqrtc <- num/denom
                c <- sqrtc^2
            }
        }
    } else { # sample size calculation based on relative effect size
        ## input checks
        if (!is.numeric(d))
            stop("d must be numeric")
        if (!is.numeric(level) || (level <= 0 || level >= 1))
            stop("level must be numeric and in (0, 1)!")
        zalpha <- qnorm(1 - level)
        zalpha <- p2z(p = level, alternative = alternative)
        c <- zalpha^2/(d^2*zo^2)
    }
    return(c)
}


## test if sampleSizeSignificance2() yields same results as 


vec01 <- c(0.001, 0.2532, 0.99)
vec01bound <- c(0, 0.0386, 0.5031, 1)
vec55 <- c(-5, -2.6288, 0, 0.0427, 4)
alternative <- c("two.sided", "one.sided", "less", "greater")
designPrior <- c("conditional", "predictive", "EB")
pars_grid_power <- expand.grid(zo=vec55,
                               power=vec01,
                               d=NA,
                               level=vec01,
                               alternative=alternative,
                               designPrior=designPrior,
                               h=abs(vec55),
                               shrinkage=vec01bound)
pars_grid_d <- expand.grid(zo=vec55,
                           power=NA,
                           d=vec01,
                           level=.025,
                           alternative="one.sided",
                           designPrior="conditional",
                           h=0,
                           shrinkage=0)
pars_grid <- rbind(pars_grid_power, pars_grid_d)

## test all configurations separately
pars_grid <- cbind(pars_grid, new=NA, legacy=NA)
sampleSizeSignificance <- ReplicationSuccess::sampleSizeSignificance
for(i in seq_len(nrow(pars_grid))){
    cat(".")
    pars_grid[i,9] <- do.call("sampleSizeSignificance2", args = pars_grid[i,1:8])
    pars_grid[i,10] <- do.call("sampleSizeSignificance", args = pars_grid[i,1:8])
}


pars_grid %>% filter(abs(new - legacy) > 0.001 |
                     (is.finite(new)  & !is.finite(new)) |
                     (!is.finite(new) & is.finite(new)))     -> problems
problems %>% nrow()
problems %>% select(-d)