tests/testthat/test-ss_best_normal.R

In ssutil: Sample Size Calculation Tools

t2_1 <-
  read.csv(
    # table 2.1 from Bechhofer et al
    textConnection(
      "groups, power , delta , n
2,0.75,0.1,91
2,0.75,0.2,23
2,0.75,0.3,11
2,0.75,0.4,6
2,0.75,0.5,4
2,0.75,0.6,3
2,0.75,0.7,2
2,0.75,0.8,2
2,0.75,0.9,2
2,0.75,1.0,1
2,0.90,0.1,329
2,0.90,0.2,83
2,0.90,0.3,37
2,0.90,0.4,21
2,0.90,0.5,14
2,0.90,0.6,10
2,0.90,0.7,7
2,0.90,0.8,6
2,0.90,0.9,5
2,0.90,1.0,4
2,0.95,0.1,542
2,0.95,0.2,136
2,0.95,0.3,61
2,0.95,0.4,34
2,0.95,0.5,22
2,0.95,0.6,16
2,0.95,0.7,12
2,0.95,0.8,9
2,0.95,0.9,7
2,0.95,1.0,6
2,0.99,0.1,1083
2,0.99,0.2,271
2,0.99,0.3,121
2,0.99,0.4,68
2,0.99,0.5,44
2,0.99,0.6,31
2,0.99,0.7,23
2,0.99,0.8,17
2,0.99,0.9,14
2,0.99,1.0,11
3,0.75,0.1,206
3,0.75,0.2,52
3,0.75,0.3,23
3,0.75,0.4,13
3,0.75,0.5,9
3,0.75,0.6,6
3,0.75,0.7,5
3,0.75,0.8,4
3,0.75,0.9,3
3,0.75,1.0,3
3,0.90,0.1,498
3,0.90,0.2,125
3,0.90,0.3,56
3,0.90,0.4,32
3,0.90,0.5,20
3,0.90,0.6,14
3,0.90,0.7,11
3,0.90,0.8,8
3,0.90,0.9,7
3,0.90,1.0,5
3,0.95,0.1,735
3,0.95,0.2,184
3,0.95,0.3,82
3,0.95,0.4,46
3,0.95,0.5,30
3,0.95,0.6,21
3,0.95,0.7,15
3,0.95,0.8,12
3,0.95,0.9,10
3,0.95,1.0,8
3,0.99,0.1,1309
3,0.99,0.2,328
3,0.99,0.3,146
3,0.99,0.4,82
3,0.99,0.5,53
3,0.99,0.6,37
3,0.99,0.7,27
3,0.99,0.8,21
3,0.99,0.9,17
3,0.99,1.0,14
4,0.75,0.1,283
4,0.75,0.2,71
4,0.75,0.3,32
4,0.75,0.4,18
4,0.75,0.5,12
4,0.75,0.6,8
4,0.75,0.7,6
4,0.75,0.8,5
4,0.75,0.9,4
4,0.75,1.0,3
4,0.90,0.1,602
4,0.90,0.2,151
4,0.90,0.3,67
4,0.90,0.4,38
4,0.90,0.5,25
4,0.90,0.6,17
4,0.90,0.7,13
4,0.90,0.8,10
4,0.90,0.9,8
4,0.90,1.0,7
4,0.95,0.1,851
4,0.95,0.2,213
4,0.95,0.3,95
4,0.95,0.4,54
4,0.95,0.5,35
4,0.95,0.6,24
4,0.95,0.7,18
4,0.95,0.8,14
4,0.95,0.9,11
4,0.95,1.0,9
4,0.99,0.1,1442
4,0.99,0.2,361
4,0.99,0.3,161
4,0.99,0.4,91
4,0.99,0.5,58
4,0.99,0.6,41
4,0.99,0.7,30
4,0.99,0.8,23
4,0.99,0.9,18
4,0.99,1.0,15
5,0.75,0.1,341
5,0.75,0.2,86
5,0.75,0.3,38
5,0.75,0.4,22
5,0.75,0.5,14
5,0.75,0.6,10
5,0.75,0.7,7
5,0.75,0.8,6
5,0.75,0.9,5
5,0.75,1.0,4
5,0.90,0.1,676
5,0.90,0.2,169
5,0.90,0.3,76
5,0.90,0.4,43
5,0.90,0.5,28
5,0.90,0.6,19
5,0.90,0.7,14
5,0.90,0.8,11
5,0.90,0.9,9
5,0.90,1.0,7
5,0.95,0.1,934
5,0.95,0.2,234
5,0.95,0.3,104
5,0.95,0.4,59
5,0.95,0.5,38
5,0.95,0.6,26
5,0.95,0.7,20
5,0.95,0.8,15
5,0.95,0.9,12
5,0.95,1.0,10
5,0.99,0.1,1537
5,0.99,0.2,385
5,0.99,0.3,171
5,0.99,0.4,97
5,0.99,0.5,62
5,0.99,0.6,43
5,0.99,0.7,32
5,0.99,0.8,25
5,0.99,0.9,19
5,0.99,1.0,16
6,0.75,0.1,388
6,0.75,0.2,97
6,0.75,0.3,44
6,0.75,0.4,25
6,0.75,0.5,16
6,0.75,0.6,11
6,0.75,0.7,8
6,0.75,0.8,7
6,0.75,0.9,5
6,0.75,1.0,4
6,0.90,0.1,735
6,0.90,0.2,184
6,0.90,0.3,82
6,0.90,0.4,46
6,0.90,0.5,30
6,0.90,0.6,21
6,0.90,0.7,15
6,0.90,0.8,12
6,0.90,0.9,10
6,0.90,1.0,8
6,0.95,0.1,998
6,0.95,0.2,250
6,0.95,0.3,111
6,0.95,0.4,63
6,0.95,0.5,40
6,0.95,0.6,28
6,0.95,0.7,21
6,0.95,0.8,16
6,0.95,0.9,13
6,0.95,1.0,10
6,0.99,0.1,1610
6,0.99,0.2,403
6,0.99,0.3,179
6,0.99,0.4,101
6,0.99,0.5,65
6,0.99,0.6,45
6,0.99,0.7,33
6,0.99,0.8,26
6,0.99,0.9,20
6,0.99,1.0,17
7,0.75,0.1,426
7,0.75,0.2,107
7,0.75,0.3,48
7,0.75,0.4,27
7,0.75,0.5,18
7,0.75,0.6,12
7,0.75,0.7,9
7,0.75,0.8,7
7,0.75,0.9,6
7,0.75,1.0,5
7,0.90,0.1,783
7,0.90,0.2,196
7,0.90,0.3,87
7,0.90,0.4,49
7,0.90,0.5,32
7,0.90,0.6,22
7,0.90,0.7,16
7,0.90,0.8,13
7,0.90,0.9,10
7,0.90,1.0,8
7,0.95,0.1,1051
7,0.95,0.2,263
7,0.95,0.3,117
7,0.95,0.4,66
7,0.95,0.5,43
7,0.95,0.6,30
7,0.95,0.7,22
7,0.95,0.8,17
7,0.95,0.9,13
7,0.95,1.0,11
7,0.99,0.1,1670
7,0.99,0.2,418
7,0.99,0.3,186
7,0.99,0.4,105
7,0.99,0.5,67
7,0.99,0.6,47
7,0.99,0.7,35
7,0.99,0.8,27
7,0.99,0.9,21
7,0.99,1.0,17
8,0.75,0.1,459
8,0.75,0.2,115
8,0.75,0.3,51
8,0.75,0.4,29
8,0.75,0.5,19
8,0.75,0.6,13
8,0.75,0.7,10
8,0.75,0.8,8
8,0.75,0.9,6
8,0.75,1.0,5
8,0.90,0.1,824
8,0.90,0.2,206
8,0.90,0.3,92
8,0.90,0.4,52
8,0.90,0.5,33
8,0.90,0.6,23
8,0.90,0.7,17
8,0.90,0.8,13
8,0.90,0.9,11
8,0.90,1.0,9
8,0.95,0.1,1096
8,0.95,0.2,274
8,0.95,0.3,122
8,0.95,0.4,69
8,0.95,0.5,44
8,0.95,0.6,31
8,0.95,0.7,23
8,0.95,0.8,18
8,0.95,0.9,14
8,0.95,1.0,11
8,0.99,0.1,1721
8,0.99,0.2,431
8,0.99,0.3,192
8,0.99,0.4,108
8,0.99,0.5,69
8,0.99,0.6,48
8,0.99,0.7,36
8,0.99,0.8,27
8,0.99,0.9,22
8,0.99,1.0,18
9,0.75,0.1,487
9,0.75,0.2,122
9,0.75,0.3,55
9,0.75,0.4,31
9,0.75,0.5,20
9,0.75,0.6,14
9,0.75,0.7,10
9,0.75,0.8,8
9,0.75,0.9,7
9,0.75,1.0,5
9,0.90,0.1,859
9,0.90,0.2,215
9,0.90,0.3,96
9,0.90,0.4,54
9,0.90,0.5,35
9,0.90,0.6,24
9,0.90,0.7,18
9,0.90,0.8,14
9,0.90,0.9,11
9,0.90,1.0,9
9,0.95,0.1,1135
9,0.95,0.2,284
9,0.95,0.3,127
9,0.95,0.4,71
9,0.95,0.5,46
9,0.95,0.6,32
9,0.95,0.7,24
9,0.95,0.8,18
9,0.95,0.9,15
9,0.95,1.0,12
9,0.99,0.1,1764
9,0.99,0.2,441
9,0.99,0.3,196
9,0.99,0.4,111
9,0.99,0.5,71
9,0.99,0.6,49
9,0.99,0.7,36
9,0.99,0.8,28
9,0.99,0.9,22
9,0.99,1.0,18
10,0.75,0.1,513
10,0.75,0.2,129
10,0.75,0.3,57
10,0.75,0.4,33
10,0.75,0.5,21
10,0.75,0.6,15
10,0.75,0.7,11
10,0.75,0.8,9
10,0.75,0.9,7
10,0.75,1.0,6
10,0.90,0.1,890
10,0.90,0.2,223
10,0.90,0.3,99
10,0.90,0.4,56
10,0.90,0.5,36
10,0.90,0.6,25
10,0.90,0.7,19
10,0.90,0.8,14
10,0.90,0.9,11
10,0.90,1.0,9
10,0.95,0.1,1169
10,0.95,0.2,293
10,0.95,0.3,130
10,0.95,0.4,74
10,0.95,0.5,47
10,0.95,0.6,33
10,0.95,0.7,24
10,0.95,0.8,19
10,0.95,0.9,15
10,0.95,1.0,12
10,0.99,0.1,1803
10,0.99,0.2,451
10,0.99,0.3,201
10,0.99,0.4,113
10,0.99,0.5,73
10,0.99,0.6,51
10,0.99,0.7,37
10,0.99,0.8,29
10,0.99,0.9,23
10,0.99,1.0,19"
    )
  )


# Note: there is a random component in the estimation of the `multz` and this
# may lead to a small differences in the calculation of the sample size
# Only test every 10 to decrease the testing time

test_that("ss_best_normal function reproduce Bechhofer table 2.1", {
  for (i in seq(1, nrow(t2_1), 10)) {
    expect_equal(
      ss_best_normal(
          ngroups=t2_1[i, "groups"], 
          dif =t2_1[i, "delta"], 
          1, power = t2_1[i, "power"],
          seed = 12345
      ),
      expected = t2_1[i, "n"],
      tolerance = 1
    )
  }
})

test_that("ss_best_normal rejects invalid group values", {
  expect_error(ss_best_normal(
    ngroups = 1,
    dif = 0.5,
    sd = 1,
    power = 0.8
  ),
  regexp = "ngroups must be >=2")
  expect_error(ss_best_normal(
    ngroups = 0,
    dif = 0.5,
    sd = 1,
    power = 0.8
  ))
  expect_error(ss_best_normal(
    ngroups = -1,
    dif = 0.5,
    sd = 1,
    power = 0.8
  ))
})

test_that("ss_best_normal rejects invalid delta", {
  expect_error(ss_best_normal(
    ngroups = 3,
    dif = 0,
    sd = 1,
    power = 0.8
  ), regexp = "dif must be > 0")
  expect_error(ss_best_normal(
    ngroups = 3,
    dif = -0.1,
    sd = 1,
    power = 0.8
  ))
})

test_that("ss_best_normal rejects invalid standard deviation", {
  expect_error(ss_best_normal(
    ngroups = 3,
    dif = 0.5,
    sd = 0,
    power = 0.8
  ),
  regexp = "sd must be > 0")
  expect_error(ss_best_normal(
    ngroups = 3,
    dif = 0.5,
    sd = -1,
    power = 0.8
  ))
})

test_that("ss_best_normal rejects invalid power values", {
  expect_error(ss_best_normal(
    ngroups = 3,
    dif = 0.5,
    sd = 1,
    power = 0
  ), regexp = "power must be > 0 and < 1")
  expect_error(ss_best_normal(
    ngroups = 3,
    dif = 0.5,
    sd = 1,
    power = 1
  ), regexp = "power must be > 0 and < 1")
  expect_error(ss_best_normal(
    ngroups = 3,
    dif = 0.5,
    sd = 1,
    power = -0.1
  ))
  expect_error(ss_best_normal(
    ngroups = 3,
    dif = 0.5,
    sd = 1,
    power = 1.1
  ))
})

test_that("ss_best_normal rejects invalid seed", {
  expect_error(ss_best_normal(
    ngroups = 3,
    dif = 0.5,
    sd = 1,
    power = 0.8,
    seed = -10
  ),
  regexp = "seed must be > 0")
  expect_error(ss_best_normal(
    ngroup = 3,
    dif = 0.5,
    sd = 1,
    power = 0.8,
    seed = 0
  ))
})

test_that("ss_best_normal and power_best_normal agrees",{
 for(poweri in c(0.5,0.7,0.8,0.9)) {
   for (difi in c(1,0.5,0.25)) {
     for(ngroupsi in c(3,5,8)) {
       for(sdi in c(2,1,0.5)){
        #cat(poweri,difi,ngroupsi,sdi,"\n")
         n <- ss_best_normal(
           power = poweri,
           dif = difi,
           sd = sdi,
           ngroups = ngroupsi
         )
         pbn <- power_best_normal(
           npergroup = n,
           dif = difi,
           sd = sdi,
           ngroups = ngroupsi
         )
         # n is an integer (ceiling), so it must reach at least the requested
         # power, but n - 1 must not (otherwise n would not be minimal)
         expect_gte(pbn, poweri)
         if (n > 1) {
           pbn_minus1 <- power_best_normal(
             npergroup = n - 1,
             dif = difi,
             sd = sdi,
             ngroups = ngroupsi
           )
           # mvtnorm's pmvnorm/qmvnorm target ~0.001 numerical accuracy, so
           # allow that much slack right at the boundary
           expect_lt(pbn_minus1, poweri + 0.001)
         }
       }
     }
   }
 }         
})