simulate.tergm() now returns $model and $monitor statistics as elemen…

…ts of a list $stats. fixes #66
statnet · Jun 22, 2022 · 32214a3 · 32214a3
1 parent ccff6eb
commit 32214a3
Show file tree

Hide file tree

Showing 6 changed files with 27 additions and 33 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: tergm
-Version: 4.1-2430
-Date: 2022-06-03
+Version: 4.1-2431
+Date: 2022-06-12
 Title: Fit, Simulate and Diagnose Models for Network Evolution Based on Exponential-Family Random Graph Models
 Authors@R: c(
         person("Pavel N.", "Krivitsky", role=c("aut","cre"), email="[email protected]", comment=c(ORCID="0000-0002-9101-3362")),

diff --git a/R/simulate.stergm.R b/R/simulate.stergm.R
@@ -227,13 +227,13 @@ simulate.network <- function(object, nsim=1, seed=NULL,
 
   if(output != "stats") {
     attributes(rv) <- c(attributes(rv), list(coef.form = coef.form, coef.diss = coef.diss))
-    stats.gen <- attr(rv, "stats.gen")
+    stats.gen <- attr(rv, "stats")$model
     if(NCOL(stats.gen) > 0) {
       attr(rv, "stats.form") <- stats.gen[,grepl("^Form~", colnames(stats.gen)) | grepl("^offset\\(Form~", colnames(stats.gen)),drop=FALSE]
       attr(rv, "stats.diss") <- stats.gen[,grepl("^Persist~", colnames(stats.gen)) | grepl("^offset\\(Persist~", colnames(stats.gen)),drop=FALSE]    
     }
   } else {
-    stats.gen <- rv$stats.gen
+    stats.gen <- rv$stats$model
     if(NCOL(stats.gen) > 0) {
       rv$stats.form <- stats.gen[,grepl("^Form~", colnames(stats.gen)) | grepl("^offset\\(Form~", colnames(stats.gen)),drop=FALSE]
       rv$stats.diss <- stats.gen[,grepl("^Persist~", colnames(stats.gen)) | grepl("^offset\\(Persist~", colnames(stats.gen)),drop=FALSE]    
@@ -270,13 +270,13 @@ simulate.networkDynamic <- function(object, nsim=1, seed=NULL,
 
   if(output != "stats") {
     attributes(rv) <- c(attributes(rv), list(coef.form = coef.form, coef.diss = coef.diss))
-    stats.gen <- attr(rv, "stats.gen")
+    stats.gen <- attr(rv, "stats")$model
     if(NCOL(stats.gen) > 0) {
       attr(rv, "stats.form") <- stats.gen[,grepl("^Form~", colnames(stats.gen)) | grepl("^offset\\(Form~", colnames(stats.gen)),drop=FALSE]
       attr(rv, "stats.diss") <- stats.gen[,grepl("^Persist~", colnames(stats.gen)) | grepl("^offset\\(Persist~", colnames(stats.gen)),drop=FALSE]    
     }
   } else {
-    stats.gen <- rv$stats.gen
+    stats.gen <- rv$stats$model
     if(NCOL(stats.gen) > 0) {
       rv$stats.form <- stats.gen[,grepl("^Form~", colnames(stats.gen)) | grepl("^offset\\(Form~", colnames(stats.gen)),drop=FALSE]
       rv$stats.diss <- stats.gen[,grepl("^Persist~", colnames(stats.gen)) | grepl("^offset\\(Persist~", colnames(stats.gen)),drop=FALSE]    

diff --git a/R/simulate.tergm.R b/R/simulate.tergm.R
@@ -126,8 +126,8 @@
 #' @return Depends on the \code{output} argument:
 #'  \item{"stats"}{If \code{stats == FALSE}, an \code{\link{mcmc}} matrix with
 #'    monitored statistics, and if \code{stats == TRUE}, a
-#'    list containing elements \code{stats} for statistics specified in the
-#'    \code{monitor} argument, and \code{stats.gen} for the model statistics.
+#'    list containing elements \code{monitor} for statistics specified in the
+#'    \code{monitor} argument, and \code{model} for the model statistics.
 #'    If \code{stats == FALSE} and no monitored statistics are specified,
 #'    an empty list is returned, with a warning.
 #'    When \code{nsim>1}, an \code{\link{mcmc.list}} (or list of them) of
@@ -149,7 +149,7 @@
 #'    \describe{
 #'      \item{\code{formula}, \code{monitor}:}{Model
 #'    and monitoring formulas used in the simulation, respectively.}
-#'      \item{\code{stats}, \code{stats.gen}:}{Network statistics as above.}
+#'      \item{\code{stats}:}{Network statistics as above.}
 #'      \item{\code{coef}:}{Coefficients used in the simulation.}
 #'      \item{\code{changes}:}{A four-column matrix summarizing the changes in the
 #'    \code{"changes"} output. (This may be removed in the future.)}
@@ -374,8 +374,7 @@ simulate_formula.network <- function(object, nsim=1, seed=NULL,
                nwd <- to.networkDynamic.lasttoggle(nw)
                nwd <- networkDynamic.apply.changes(nwd,changes)
                attributes(nwd) <- c(attributes(nwd), # Don't clobber existing attributes!
-                                    list(stats.gen = stats.gen,
-                                         stats = stats.mon,
+                                    list(stats = list(model = stats.gen, monitor = stats.mon),
                                          coef = coef,
                                          changes = changes))
                nwd <- .add.net.obs.period.spell(nwd, start-1+time.offset, time.slices)
@@ -388,15 +387,14 @@ simulate_formula.network <- function(object, nsim=1, seed=NULL,
                # assume that simulate.stergm has added +1 to all the time values, so subtract 1 for an offset of 0
                changes[,1]<-changes[,1]-1+time.offset
                attributes(changes) <- c(attributes(changes), # Don't clobber existing attributes!
-                                        list(stats.gen = stats.gen,
-                                             stats = stats.mon,
+                                        list(stats = list(model = stats.gen, monitor = stats.mon),
                                              coef = coef,
                                              start = nw%n%"time" + 0,
                                              end = nw%n%"time" + time.slices))
                changes
              },
              stats = {
-               list(stats.gen = stats.gen, stats = stats.mon)
+               list(model = stats.gen, monitor = stats.mon)
              },
              final = {
                changes <- z$changed
@@ -407,8 +405,7 @@ simulate_formula.network <- function(object, nsim=1, seed=NULL,
                newnw <- as.network(z$newnetwork)
                newnw <- .add.net.obs.period.spell(newnw, start-1+time.offset, time.slices)
                attributes(newnw) <- c(attributes(newnw), # Don't clobber existing attributes!
-                                      list(stats.gen = stats.gen,
-                                           stats = stats.mon,
+                                      list(stats = list(model = stats.gen, monitor = stats.mon),
                                            coef = coef,
                                            start = nw%n%"time" + 0,
                                            end = nw%n%"time" + time.slices,
@@ -423,8 +420,7 @@ simulate_formula.network <- function(object, nsim=1, seed=NULL,
                changes[,1]<-changes[,1]-1+time.offset
                newnw <- z$newnetwork
                attributes(newnw) <- c(attributes(newnw), # Don't clobber existing attributes!
-                                      list(stats.gen = stats.gen,
-                                           stats = stats.mon,
+                                      list(stats = list(model = stats.gen, monitor = stats.mon),
                                            coef = coef,
                                            start = nw%n%"time" + 0,
                                            end = nw%n%"time" + time.slices,
@@ -557,8 +553,7 @@ simulate_formula.networkDynamic <- function(object, nsim=1, seed=NULL,
   }
 
   attributes(object) <- c(attributes(object), # Don't clobber existing attributes!
-                          list(stats.gen = attr(sim, "stats.gen"),
-                               stats = attr(sim, "stats"),
+                          list(stats = attr(sim, "stats"),
                                coef = coef,
                                changes = rbind(attr(object,"changes"),matrix(c(sim), nrow=nrow(sim),ncol=ncol(sim),dimnames=list(rownames(sim),colnames(sim))))
                                ))

diff --git a/tests/basic_tests.R b/tests/basic_tests.R
@@ -31,8 +31,7 @@ stopifnot(all(summary(nw~edges+Form(~edges)+Persist(~edges))==c(16,16,16)))
 
 # Now, let's start at Time 2, add edge (1,2) at time 4, remove it at time 5, and then delete edge (1,4) at time 7 and re-add it at time 8.
 o <- tergm.godfather(nw~edges+Form(~edges)+Persist(~edges), toggles=rbind(c(4,1,2),c(5,1,2),c(7,1,4),c(8,1,4)),start=2, end=10, end.network=TRUE)
-attr(o, "stats") # Statistics are appropriate: note how both formation and dissolution "lag":
-stopifnot(all(attr(o,"stats")==structure(c(16,17,16,16,15,16,16,16,
+stopifnot(all(attr(o,"stats")$monitor==structure(c(16,17,16,16,15,16,16,16,
                                            16,17,17,16,16,16,16,16,
                                            16,16,16,16,15,15,16,16),
                                          .Dim=c(8L,3L))))

diff --git a/tests/sim_gf_sum.R b/tests/sim_gf_sum.R
@@ -45,17 +45,17 @@ simtest <- function(S, edges, dur, n, dir=FALSE, bip=0){
   # Simulate. Starting from an ordinary network:
   dynsim<-simulate(g1 ~ Form(~edges) + Persist(~edges),coef=c(coef.form,coef.diss),time.slices=S,verbose=TRUE,monitor=~edgecov("dc")+edgecov.ages("dc"), dynamic=TRUE)
 
-  sim.stats <- attr(dynsim, "stats")
+  sim.stats <- attr(dynsim, "stats")$monitor
 
   # Resuming from a networkDynamic:
   dynsim2<-simulate(dynsim ~ Form(~edges) + Persist(~edges),coef=c(coef.form,coef.diss),time.slices=S,verbose=TRUE,monitor=~edgecov("dc")+edgecov.ages("dc"), dynamic=TRUE)
 
-  sim.stats <- rbind(sim.stats, attr(dynsim2, "stats"))
+  sim.stats <- rbind(sim.stats, attr(dynsim2, "stats")$monitor)
 
   # Resuming from a resumed networkDynamic:
   dynsim3<-simulate(dynsim2 ~ Form(~edges) + Persist(~edges),coef=c(coef.form,coef.diss),time.slices=S,verbose=TRUE,monitor=~edgecov("dc")+edgecov.ages("dc"), dynamic=TRUE)
 
-  sim.stats <- rbind(sim.stats, attr(dynsim3, "stats"))
+  sim.stats <- rbind(sim.stats, attr(dynsim3, "stats")$monitor)
 
   # Replay the simulation using a networkDynamic:
   gf1.stats <- as.matrix(tergm.godfather(dynsim3~edgecov("dc")+edgecov.ages("dc"), start=0, end=S*3))

diff --git a/tests/testthat/test-simulate.R b/tests/testthat/test-simulate.R
@@ -202,12 +202,12 @@ test_that("simulate.network behaves reasonably", {
   expect_equal(old_nw6, new_nw6, check.attributes = FALSE)
 
   expect_identical(attr(old_stats_nw, "stats"), attr(new_stats_nw, "stats"))
-  expect_identical(attr(old_stats_nw, "stats.gen"), attr(new_stats_nw, "stats.gen"))
-  expect_identical(cbind(attr(old_stats_nw, "stats.form"), attr(old_stats_nw, "stats.diss")), as.matrix(attr(new_stats_nw, "stats.gen")))
+  ## expect_identical(attr(old_stats_nw, "stats.gen"), attr(new_stats_nw, "stats.gen"))
+  expect_identical(cbind(attr(old_stats_nw, "stats.form"), attr(old_stats_nw, "stats.diss")), as.matrix(attr(new_stats_nw, "stats")$model))
 
   expect_identical(old_stats$stats, new_stats$stats)
-  expect_identical(old_stats$stats.gen, new_stats$stats.gen)
-  expect_identical(cbind(old_stats$stats.form, old_stats$stats.diss), as.matrix(new_stats$stats.gen))
+  ## expect_identical(old_stats$stats.gen, new_stats$stats.gen)
+  expect_identical(cbind(old_stats$stats.form, old_stats$stats.diss), as.matrix(new_stats$stats)$model)
 })
 
 test_that("simulate.networkDynamic behaves reasonably", {
@@ -264,12 +264,12 @@ test_that("simulate.networkDynamic behaves reasonably", {
   expect_equal(attr(old_nwD8, "stats"), attr(new_nwD8, "stats"))
 
   expect_equal(attr(old_stats_nwD, "stats"), attr(new_stats_nwD, "stats"))
-  expect_equal(attr(old_stats_nwD, "stats.gen"), attr(new_stats_nwD, "stats.gen"))
-  expect_equal(cbind(attr(old_stats_nwD, "stats.form"), attr(old_stats_nwD, "stats.diss")), as.matrix(attr(new_stats_nwD, "stats.gen")))
+  ## expect_equal(attr(old_stats_nwD, "stats.gen"), attr(new_stats_nwD, "stats.gen"))
+  expect_equal(cbind(attr(old_stats_nwD, "stats.form"), attr(old_stats_nwD, "stats.diss")), as.matrix(attr(new_stats_nwD, "stats")$model))
 
   expect_equal(old_stats$stats, new_stats$stats)
-  expect_equal(old_stats$stats.gen, new_stats$stats.gen)
-  expect_equal(cbind(old_stats$stats.form, old_stats$stats.diss), as.matrix(new_stats$stats.gen))
+  ## expect_equal(old_stats$stats.gen, new_stats$stats.gen)
+  expect_equal(cbind(old_stats$stats.form, old_stats$stats.diss), as.matrix(new_stats$stats$model))
 
   expect_equal(old_nwD_constr, new_nwD_constr, check.attributes = FALSE)
   expect_equal(old_nwD_constr2, new_nwD_constr2, check.attributes = FALSE)