tests/measles.R

In spatPomp: Inference for Spatiotemporal Partially Observed Markov Processes

library(spatPomp)
set.seed(1)

## a default measles model with fixed shared initial value parameters (IVPs)
## compiled into the object

m_U <- 2 
m_model <- measles(U=m_U)

range(time(m_model))

m <- window(m_model,1950,1955)

m_params <- c(R0=56.8,amplitude=0.554,gamma=30.4,sigma=28.9,mu=0.02,sigmaSE=0.02,rho=0.488,psi=0.116,g=100,alpha=1,iota=0,cohort=0.1)

m_pf <- pfilter(m_model,Np=20,params=m_params)
logLik(m_pf)

## setting IVPs

m2_model <- measles(U=m_U,fixed_ivps=FALSE)
m2_model <- window(m2_model,1950,1955)

m2_shared_params <- c(R0=56.8,amplitude=0.554,gamma=30.4,sigma=28.9,mu=0.02,sigmaSE=0.02,rho=0.488,psi=0.116,g=100,alpha=1,iota=0,cohort=0.1)
m2_unit_ivp_names <- c('S','E','I')
m2_statenames <- paste0(rep(m2_unit_ivp_names,each=m_U),1:m_U)
m2_ivp_names <- paste0(m2_statenames,"_0")
m2_ivps <- rep(c(0.032,0.00005,0.00004),each=m_U)
names(m2_ivps) <- m2_ivp_names
m2_params <- c(m2_shared_params,m2_ivps)

set.seed(1)
m2_pf <- pfilter(m2_model,Np=20,params=m2_params)
logLik(m2_pf)

##
## Note: the measles skeleton is correct only when there is no cohort effect,
## i.e., cohort=0.
##

## A call to igirf using the moment-based guide function can test
## compiled code for eunit_measure, munit_measure, vunit_measure,
## dunit_measure, runit_measure, rprocess, skeleton, rinit and partrans. 
## It also tests igirf (method=moment) with covariates

m3_params <- m_params
m3_params["cohort"] <- 0
m3_igirf_lookahead <- 1
m3_igirf_ninter <- 2
m3_igirf_np <- 5
m3_igirf_nguide <- 5
m3_igirf_ngirf <- 2

suppressWarnings(
  m3_igirf_out <- igirf(m_model, Ngirf = m3_igirf_ngirf,
    params=m3_params,
    rw.sd=rw_sd(g=0.02),
    cooling.type = "geometric",
    cooling.fraction.50 = 0.5,
    Np=m3_igirf_np,
    Ninter = m3_igirf_ninter,
    lookahead = m3_igirf_lookahead,
    Nguide = m3_igirf_nguide,
    kind = 'moment',
    verbose = FALSE
  )
)
logLik(m3_igirf_out)

## another call to igirf tests igirf (method=bootstrap) with covariates
suppressWarnings(m3_igirf_bootstrap_out <- igirf(m3_igirf_out,kind='bootstrap'))
logLik(m3_igirf_bootstrap_out)

## test error message
try(measles(U=1000))


## --------------------------------------------------------------------
## using measles to test spatPomp() covariate replacement functionality
## ____________________________________________________________________

m_covrep <- spatPomp(m,
  covar=spatPomp::measlesUK[
    measlesUK$city %in% m@unit_names,
    c("year","city","pop","births")]
)
for(slt in slotNames(m)) if(!identical(slot(m,slt),slot(m_covrep,slt))) print(slt)
## covariates are not the same since this replacement does not deal with
## birth delay