R Under development (unstable) (2024-02-17 r85935 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library(spatPomp) Loading required package: pomp > set.seed(22) > > #model_type <- "he10" > model_type <- "mostly shared" > parNames <- c("alpha","R0","g","sigma","gamma","amplitude","cohort","sigmaSE","S_0","E_0","I_0","rho","psi","iota","mu") > # he10 defaults to alpha=1, cohort=0, which means the usual transformations are undefined. > # here, we don't estimate either > if(model_type == "mostly fixed"){ + sharedParNames <- c("R0","psi") + unitParNames <- c("rho","S_0") + estParNames <- c(sharedParNames,unitParNames) + fixedParNames <- setdiff(parNames,estParNames) + } else if(model_type == "mostly shared"){ + sharedParNames <- c("R0","psi","g","sigma","gamma","amplitude","sigmaSE") + unitParNames <- c("rho","S_0","E_0","I_0") + estParNames <- c(sharedParNames,unitParNames) + fixedParNames <- setdiff(parNames,estParNames) + } else if(model_type == "plausible parameters shared"){ + # parameters are shared when that makes mechanistic sense. + sharedParNames <- c("R0","g","sigma","gamma","amplitude") + unitParNames <- c("sigmaSE","S_0","E_0","I_0","rho","psi") + estParNames <- c(sharedParNames,unitParNames) + fixedParNames <- setdiff(parNames,estParNames) + } else if(model_type == "all unit-specific"){ + # all parameters estimated except life expecancy + # and immigration, which should not be needed when there is coupling + fixedParNames <- c("mu","iota") + sharedParNames <- NULL + unitParNames <- setdiff(parNames,fixedParNames) + estParNames <- c(sharedParNames,unitParNames) + } else if(model_type == "he10"){ + # all the parameters estimated by He et al (2010) Table 2 + fixedParNames <- c("mu","g") + sharedParNames <- NULL + unitParNames <- setdiff(parNames,fixedParNames) + estParNames <- c(sharedParNames,unitParNames) + } > > # Note: here we assume that there are no unestimated unit-specific > # parameters. That could readily be accommodated if needed. > > h_model <- he10(U=2,dt=4/365,Tmax=1950.5, + expandedParNames=estParNames) > > coef(h_model) alpha1 iota1 cohort1 mu1 R01 R02 psi1 1.0e+00 0.0e+00 0.0e+00 2.0e-02 3.0e+01 3.0e+01 1.5e-01 psi2 g1 g2 sigma1 sigma2 gamma1 gamma2 1.5e-01 4.0e+02 4.0e+02 5.2e+01 5.2e+01 5.2e+01 5.2e+01 amplitude1 amplitude2 sigmaSE1 sigmaSE2 rho1 rho2 S_01 5.0e-01 5.0e-01 1.5e-01 1.5e-01 5.0e-01 5.0e-01 3.2e-02 S_02 E_01 E_02 I_01 I_02 3.2e-02 5.0e-05 5.0e-05 4.0e-05 4.0e-05 > > paste("bpfilter logLik for he10 model:",logLik(bpfilter(h_model,Np=10,block_size=1))) [1] "bpfilter logLik for he10 model: -3827.56515459983" > > > h_U <- length(unit_names(h_model)) > > ivpParNames <- c("S_0","E_0","I_0") > ivpEstParNames <- intersect(ivpParNames,estParNames) > regEstParNames <- setdiff(estParNames,ivpParNames) > > estParNames_expanded <- unlist(lapply(estParNames,function(x)paste0(x,1:h_U))) > regEstParNames_expanded <- unlist(lapply(regEstParNames,function(x)paste0(x,1:h_U))) > ivpEstParNames_expanded <- unlist(lapply(ivpEstParNames,function(x)paste0(x,1:h_U))) > fixedParNames_expanded <- paste0(fixedParNames,1) > > > reg_rw.sd <- rep(list(0.02),times=length(regEstParNames_expanded)) > names(reg_rw.sd) <- regEstParNames_expanded > if("alpha"%in%estParNames) reg_rw.sd[paste0("alpha",1:h_U)] <- 0.005 > > ivp_rw.sd <- lapply(ivpEstParNames_expanded,function(x)expression(ivp(0.05))) > names(ivp_rw.sd) <- ivpEstParNames_expanded > h_rw.sd <- do.call(rw_sd,c(reg_rw.sd,ivp_rw.sd)) > > all_units = seq_len(length(unit_names(h_model))) > nblocks = 2 > block_list = split(all_units, sort(all_units %% nblocks)) > block_list <- lapply(block_list, as.integer) > > h_ibpf <- ibpf(h_model, + params=coef(h_model), + sharedParNames=sharedParNames, + unitParNames=unitParNames, + Nbpf=2, + spat_regression=0.1, + Np=10, + rw.sd=h_rw.sd, + cooling.fraction.50=0.5, + block_list=block_list + ) > > paste("ibpf logLik for he10 model:",logLik(bpfilter(h_ibpf,Np=10,block_size=1))) [1] "ibpf logLik for he10 model: -6589.70922701501" > > coef(h_ibpf) alpha1 iota1 cohort1 mu1 R01 R02 1.000000e+00 0.000000e+00 0.000000e+00 2.000000e-02 3.194959e+01 3.194959e+01 psi1 psi2 g1 g2 sigma1 sigma2 1.489545e-01 1.489545e-01 4.104272e+02 4.104272e+02 4.821791e+01 4.821791e+01 gamma1 gamma2 amplitude1 amplitude2 sigmaSE1 sigmaSE2 4.722066e+01 4.722066e+01 4.736288e-01 4.736288e-01 1.461788e-01 1.461788e-01 rho1 rho2 S_01 S_02 E_01 E_02 5.423109e-01 5.009793e-01 3.257674e-02 3.188634e-02 4.835306e-05 4.617649e-05 I_01 I_02 3.749471e-05 4.049693e-05 > > ## test error message when munit_measure is undefined > try(girf(h_model,kind="moment", + Np=10,Ninter=2,Nguide=10,lookahead=1,tol=1e-5)) Error in .local(object, ...) : girf with kind = 'moment' requires munit_measure > > # Create second ibpfd_spatPomp object with different chain length, to test error > h_ibpf2 <- ibpf(h_model, + params=coef(h_model), + sharedParNames=sharedParNames, + unitParNames=unitParNames, + Nbpf=3, + spat_regression=0.1, + Np=10, + rw.sd=h_rw.sd, + cooling.fraction.50=0.5, + block_list=block_list + ) > > # Should correctly make ibpfList object > is(c(h_ibpf, h_ibpf), "ibpfList") [1] TRUE > > # Throws error because they have different chain lengths > try(c(h_ibpf, h_ibpf2)) Error in validObject(.Object) : invalid class "ibpfList" object: error in 'c': to be combined, 'ibpfd_spatPomp' objects must have chains of equal length > > > proc.time() user system elapsed 3.09 0.14 6.84