add quantile 95% CI in summary

R/summary.R

@@ -9,6 +9,7 @@
 #' \item \code{mean}: (t-dimensional vector) Mean of simulated counts of the compartment;
 #' \item \code{sd}: (t-dimensional vector) Standard deviantion of simulated counts of the compartment;
 #' \item \code{CI}: (t by 2 matrix) 95\% confidence intervals of simulated counts of the compartment (assumed Gaussian),
+#' \item \code{qntCI}: (t by 2 matrix) 95\% quantile confidence intervals,
 #' }
 #' where t is the total number of time points in the simulation (i.e. object$nstep).
 #'
@@ -20,12 +21,13 @@ summary_seiqhrf <- function(object){
     
     #### extract mean and sd
     df_mean <- as.data.frame(object, out = "mean")
-    df_mean <- df_mean[, names(df_mean) %in% prev_names]     ### subset columns
-    df_mean <- df_mean[, order(match(names(df_mean), prev_names))]  ### order columns
+    sub_cols <- names(df_mean)[ order(match(names(df_mean), prev_names))][1:length(prev_names)]
+    df_mean <- df_mean[, sub_cols]  
     
-    df_sd <- as.data.frame(object, out = "sd")
-    df_sd <- df_sd[, names(df_sd) %in% prev_names]    
-    df_sd <- df_sd[, order(match(names(df_sd), prev_names))]
+    df_sd <- as.data.frame(object, out = "sd")[, sub_cols]
+
+    df_qnt_left <- as.data.frame(object, out = "qnt", qval = .025)[, sub_cols]
+    df_qnt_right <- as.data.frame(object, out = "qnt", qval = .975)[, sub_cols]
     
     #### prepare result 
     res_names <- c("Susceptible", "Exposed", "Infectious", "Quarantine", "Hospitalisation", "Recovered", "Fatal")
@@ -36,13 +38,16 @@ summary_seiqhrf <- function(object){
         prev_mean <- df_mean[[prev_no]]
         prev_sd <- df_sd[[prev_no]]
         prev_ci <- cbind(prev_mean - 1.96*prev_sd, prev_mean + 1.96*prev_sd)
+        prev_qci <- cbind(df_qnt_left[[prev_no]], df_qnt_right[[prev_no]])
         names(prev_mean) <- c(1:nsteps)
         names(prev_sd) <- c(1:nsteps)
         rownames(prev_ci) <- c(1:nsteps)
+        rownames(prev_qci) <- c(1:nsteps)
         
         res[[prev_no]]$"mean" <- prev_mean
         res[[prev_no]]$"sd" <- prev_sd
         res[[prev_no]]$"CI" <- prev_ci
+        res[[prev_no]]$"qntCI" <- prev_qci
        
     }
 

man/merge.seiqhrf.icm.Rd

 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/icm_seiqhrf.R
+% Please edit documentation in R/seiqhrf.R
 \name{merge.seiqhrf.icm}
 \alias{merge.seiqhrf.icm}
 \title{Merge icm}

man/seiqhrf.Rd

 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/icm_seiqhrf.R
+% Please edit documentation in R/seiqhrf.R
 \name{seiqhrf}
 \alias{seiqhrf}
 \title{SEIQHRF function}

man/summary_seiqhrf.Rd

@@ -15,6 +15,7 @@ A list of compartments, each compartment is a list of three:
 \item \code{mean}: (t-dimensional vector) Mean of simulated counts of the compartment;
 \item \code{sd}: (t-dimensional vector) Standard deviantion of simulated counts of the compartment;
 \item \code{CI}: (t by 2 matrix) 95\% confidence intervals of simulated counts of the compartment (assumed Gaussian),
+\item \code{qntCI}: (t by 2 matrix) 95\% quantile confidence intervals,
 }
 where t is the total number of time points in the simulation (i.e. object$nstep).
 }

tests/testthat/test-icm-seiqhrf.R → tests/testthat/test-seiqhrf.R

-test_that("Identical output as Churches' original function: icm_seiqhrf", {
+test_that("Identical output as Churches' original function: seiqhrf", {
     
     full_params <- set_param(s.num = 1000, nsteps = 10, arec.rate = 0)
     param <- full_params$param