R Under development (unstable) (2025-09-25 r88874 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2025 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.

> # Load the EpiILM package
> library("EpiILM")
Loading required package: coda
Loading required package: adaptMCMC
Loading required package: parallel
Loading required package: Matrix
Loading required package: LaplacesDemon
> set.seed(101)
> 
> # sample size
> n <- 500
> 
> # Generating an undirected binary contact network:
> contact <- matrix(0, n, n)
> for(i in 1:(n-1)) {
+ 	contact[i, ((i+1):n)] <- rbinom((n-i), 1, 0.05)
+ 	contact[((i+1):n), i] <- contact[i, ((i+1):n)]
+ 	}
> 
> # Generating the susceptibility and transmissibilty covariates:
> X1 <- round(rexp(n, 1/100))
> X2 <- round(rgamma(n, 50, 0.5))
> # Simulate epidemic form SIR network-based ILM
> 
> infp <- rep(3, n)
> SIR.net <- epidata(type = "SIR", n = 500, tmax = 50,
+      sus.par = c(0.003, 0.002), trans.par = c(0.0003, 0.0002),
+      contact = contact, infperiod = infp,
+      Sformula = ~ -1 + X1 + X2, Tformula = ~ -1 + X1 + X2)
> 
> # epimcmc function to estimate the model parameters:
> t_end <- max(SIR.net$inftime)
> prior_par <- matrix(rep(1, 4), ncol = 2, nrow = 2)
> 
> # This took 305.7 seconds on an Apple MacBook Pro with i5-core Intel 2.4 GHz
> # processors with 8 GB of RAM.
> mcmcout_SIR.net <- epimcmc(SIR.net, tmax = t_end, niter = 100,
+         Sformula = ~-1 + X1 + X2, Tformula = ~-1 + X1 + X2,
+         sus.par.ini = c(0.003, 0.01), trans.par.ini = c(0.01, 0.01),
+         pro.sus.var = c(0.0, 0.1), pro.trans.var = c(0.1, 0.1),
+         prior.sus.dist = c("gamma", "gamma"), prior.trans.dist = c("gamma", "gamma"),
+         prior.sus.par = prior_par, prior.trans.par = prior_par,
+         adapt = TRUE, acc.rate = 0.5)
  generate 100 samples 
> 
> # Summary of MCMC results
> summary(mcmcout_SIR.net)
Model: SIR network-based discrete-time ILM 
Method: Markov chain Monte Carlo (MCMC) 

Iterations = 1:100
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 100 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

            Mean        SD  Naive SE Time-series SE
alpha.2 0.009838 0.0011424 1.142e-04      1.555e-04
phi.1   0.010086 0.0006051 6.051e-05      8.237e-05
phi.2   0.010466 0.0032753 3.275e-04      4.458e-04

2. Quantiles for each variable:

        2.5%  25%  50%  75% 97.5%
alpha.2 0.01 0.01 0.01 0.01  0.01
phi.1   0.01 0.01 0.01 0.01  0.01
phi.2   0.01 0.01 0.01 0.01  0.01

> 
> # Posterior predictive forecasting
> # Posterior prediction starting from time point 5
> pred.SIR.net.point.5 <- pred.epi(SIR.net, xx = mcmcout_SIR.net, tmin = 5,
+          Sformula = ~-1 + X1 + X2, Tformula = ~-1 + X1 + X2,
+          criterion = "newly infectious",
+          n.samples = 50)
  generate 50 epidemics 
> 
> # Posterior prediction starting from time point 14
> pred.SIR.net.point.8 <- pred.epi(SIR.net, xx = mcmcout_SIR.net, tmin = 14,
+         Sformula = ~-1 + X1 + X2, Tformula = ~-1 + X1 + X2,
+         criterion = "newly infectious",
+         n.samples = 50)
  generate 50 epidemics 
> 
> 
> proc.time()
   user  system elapsed 
   2.37    0.51    2.89