tests/testthat/test-ABCSMC.R

In tjmckinley/SimBIID: Simulation-Based Inference Methods for Infectious Disease Models

test_that("ABCSMC works", {
    
    ## set up SIR simulationmodel
    transitions <- c(
        "S -> beta * S * I -> I",
        "I -> gamma * I -> R"
    )
    compartments <- c("S", "I", "R")
    pars <- c("beta", "gamma")
    model <- mparseRcpp(
        transitions = transitions,
        compartments = compartments,
        pars = pars
    )
    model <- compileRcpp(model)

    ## generate function to run simulators
    ## and return summary statistics
    simSIR <- function(pars, data, tols, u, model) {

        ## run model
        sims <- model(pars, 0, data[2] + tols[2], u)

        ## this returns a vector of the form:
        ## completed (1/0), t, S, I, R (here)
        if(sims[1] == 0) {
            ## if simulation rejected
            return(NA)
        } else {
            ## extract finaltime and finalsize
            finaltime <- sims[2]
            finalsize <- sims[5]
        }

        ## return vector if match, else return NA
        if(all(abs(c(finalsize, finaltime) - data) <= tols)){
            return(c(finalsize, finaltime))
        } else {
            return(NA)
        }
    }

    ## set priors
    priors <- data.frame(
        parnames = c("beta", "gamma"),
        dist = rep("gamma", 2),
        stringsAsFactors = FALSE
    )
    priors$p1 <- c(10, 10)
    priors$p2 <- c(10^4, 10^2)

    ## define the targeted summary statistics
    data <- c(
        finalsize = 30,
        finaltime = 76
    )

    ## set initial states (1 initial infection
    ## in population of 120)
    iniStates <- c(S = 119, I = 1, R = 0)

    ## set initial tolerances
    tols <- c(
        finalsize = 50,
        finaltime = 50
    )
    
    set.seed(50)

    ## run 2 generations of ABC-SMC
    ## setting tolerance to be 50th
    ## percentile of the accepted
    ## tolerances at each generation
    post <- ABCSMC(
        x = data,
        priors = priors,
        func = simSIR,
        u = iniStates,
        tols = tols,
        ptol = 0.2,
        ngen = 2,
        npart = 50,
        model = model
    )
    ## the first run always succeeds, but warns
    expect_snapshot_output(post)

    ## run one further generation
    post <- ABCSMC(post, ptols = 0.5, ngen = 1)
    ## the first run always succeeds, but warns
    expect_snapshot_output(post)
    
    ## the first run always succeeds, but warns
    expect_snapshot_output(summary(post))
    
    ## try to run 2 generations of ABC-SMC using fixed tolerances
    tols <- matrix(c(50, 50, 50, 50), 2, 2, byrow = TRUE)
    colnames(tols) <- c("finalsize", "finaltime")
    expect_error(ABCSMC(
        x = data,
        priors = priors,
        func = simSIR,
        u = iniStates,
        tols = tols,
        ngen = 2,
        npart = 50,
        model = model
    ))
    tols <- matrix(c(50, 50, 51, 20), 2, 2, byrow = TRUE)
    colnames(tols) <- c("finalsize", "finaltime")
    expect_error(ABCSMC(
        x = data,
        priors = priors,
        func = simSIR,
        u = iniStates,
        tols = tols,
        ngen = 2,
        npart = 50,
        model = model
    ))
    
    ## run 2 generations of ABC-SMC using fixed tolerances
    tols <- matrix(c(50, 50, 50, 20), 2, 2, byrow = TRUE)
    colnames(tols) <- c("finalsize", "finaltime")
    post <- ABCSMC(
        x = data,
        priors = priors,
        func = simSIR,
        u = iniStates,
        tols = tols,
        ngen = 2,
        npart = 50,
        model = model
    )
    ## the first run always succeeds, but warns
    expect_snapshot_output(post)
    
    ## run one further generation
    newtols <- c(finalsize = 20, finaltime = 15)
    post <- ABCSMC(post, tols = newtols)
    ## the first run always succeeds, but warns
    expect_snapshot_output(post)
    
    ## run one further generation
    expect_error(ABCSMC(post, tols = newtols))
    newtols <- c(finalsize = 21, finaltime = 15)
    expect_error(ABCSMC(post, tols = newtols))
    
    ## run further generation using ptols
    expect_error(ABCSMC(post, tols = newtols, ptols = 0.8, ngen = 1))
    post <- ABCSMC(post, ptols = 0.8, ngen = 1)
    ## the first run always succeeds, but warns
    expect_snapshot_output(post)
    
    ## run 2 generations of ABC-SMC using fixed tolerances
    tols <- matrix(c(50, 50), 1, 2, byrow = TRUE)
    colnames(tols) <- c("finalsize", "finaltime")
    post <- ABCSMC(
        x = data,
        priors = priors,
        func = simSIR,
        u = iniStates,
        tols = tols,
        npart = 50,
        model = model
    )
    
    ## run two further generations
    newtols <- matrix(c(50, 50, 40, 40), 2, 2, byrow = TRUE)
    colnames(newtols) <- c("finalsize", "finaltime")
    expect_error(ABCSMC(post, tols = newtols))
    
    ## run two further generations
    newtols <- matrix(c(45, 45, 45, 45), 2, 2, byrow = TRUE)
    colnames(newtols) <- c("finalsize", "finaltime")
    expect_error(ABCSMC(post, tols = newtols))
    
    ## run two further generations
    newtols <- matrix(c(45, 45, 40, 40), 2, 2, byrow = TRUE)
    colnames(newtols) <- c("finalsize", "finaltime")
    post <- ABCSMC(post, tols = newtols)
    
    ## test minimum tolerances in various situations
    tols <- matrix(c(50, 50), 1, 2, byrow = TRUE)
    colnames(tols) <- c("finalsize", "finaltime")
    expect_error(ABCSMC(
        x = data,
        priors = priors,
        func = simSIR,
        u = iniStates,
        tols = tols,
        npart = 50,
        model = model,
        mintols = c(40, 40)
    ))
    mintols <- matrix(c(40, 40), 1, 2, byrow = TRUE)
    colnames(mintols) <- c("finalsize", "finaltime")
    expect_error(ABCSMC(
        x = data,
        priors = priors,
        func = simSIR,
        u = iniStates,
        tols = tols,
        npart = 50,
        model = model,
        mintols = mintols
    ))
    mintols <- c(finalsize = 40, finaltime = 40)
    post <- ABCSMC(
        x = data,
        priors = priors,
        func = simSIR,
        u = iniStates,
        tols = tols,
        npart = 50,
        model = model,
        mintols = mintols
    )
    mintols <- c(finalsize = 50, finaltime = 40)
    post <- ABCSMC(
        x = data,
        priors = priors,
        func = simSIR,
        u = iniStates,
        tols = tols,
        npart = 50,
        model = model,
        mintols = mintols
    )
    mintols <- c(finalsize = 51, finaltime = 50)
    expect_error(ABCSMC(
        x = data,
        priors = priors,
        func = simSIR,
        u = iniStates,
        tols = tols,
        npart = 50,
        model = model,
        mintols = mintols
    ))
    
    ## run two further generations
    newtols <- matrix(c(40, 40, 38, 38), 2, 2, byrow = TRUE)
    colnames(newtols) <- c("finalsize", "finaltime")
    expect_error(ABCSMC(post, tols = newtols))
    
    ## run two further generations
    mintols <- c(40, 40)
    names(mintols) <- c("finalsize", "finaltime")
    expect_error(ABCSMC(post, tols = newtols, mintols = mintols))
    
    ## run two further generations
    mintols <- c(35, 35)
    names(mintols) <- c("finalsize", "finaltime")
    post <- ABCSMC(post, tols = newtols, mintols = mintols)
    
    ## run two further generations
    post <- ABCSMC(post, ptols = 0.5, ngen = 1)
    expect_error(ABCSMC(post, ptols = 0.5, ngen = 1))
    mintols <- c(30, 30)
    names(mintols) <- c("finalsize", "finaltime")
    post <- ABCSMC(post, ptols = 0.5, ngen = 1, mintols = mintols)
    
    ## check reproducibility in serial
    
    tols <- c(finalsize = 50, finaltime = 50)
    set.seed(50)
    post <- ABCSMC(
        x = data,
        priors = priors,
        func = simSIR,
        u = iniStates,
        tols = tols,
        npart = 10,
        model = model
    )
    set.seed(50)
    post1 <- ABCSMC(
        x = data,
        priors = priors,
        func = simSIR,
        u = iniStates,
        tols = tols,
        npart = 10,
        model = model
    )
    
    expect_identical(post, post1)
    
    ## check reproducibility in parallel
    
    if(.Platform$OS.type != "windows" & requireNamespace("parallel", quietly = TRUE)) {
        tols <- c(finalsize = 50, finaltime = 50)
        set.seed(50)
        post <- ABCSMC(
            x = data,
            priors = priors,
            func = simSIR,
            u = iniStates,
            tols = tols,
            npart = 10,
            model = model,
            parallel = TRUE,
            mc.cores = 2
        )
        set.seed(50)
        post1 <- ABCSMC(
            x = data,
            priors = priors,
            func = simSIR,
            u = iniStates,
            tols = tols,
            npart = 10,
            model = model,
            parallel = TRUE,
            mc.cores = 2
        )
        
        expect_identical(post, post1)
        
        ## check for changing seeds in parallel
        
        tols <- c(finalsize = 50, finaltime = 50)
        set.seed(50)
        post <- ABCSMC(
            x = data,
            priors = priors,
            func = simSIR,
            u = iniStates,
            tols = tols,
            npart = 10,
            model = model,
            parallel = TRUE,
            mc.cores = 2
        )
        post1 <- ABCSMC(
            x = data,
            priors = priors,
            func = simSIR,
            u = iniStates,
            tols = tols,
            npart = 10,
            model = model,
            parallel = TRUE,
            mc.cores = 2
        )
        
        expect_false(identical(post, post1))
    }
    
})