Merge pull request #70 from snlab-nl/develop

Develop

DESCRIPTION

@@ -2,8 +2,8 @@ Encoding: UTF-8
 Package: RSiena
 Type: Package
 Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.3.19
-Date: 2023-02-07
+Version: 1.3.20
+Date: 2023-04-22
 Authors@R: c(person("Tom A.B.", "Snijders", role = c("cre", "aut"), email = "tom.snijders@nuffield.ox.ac.uk", comment = c(ORCID = "0000-0003-3157-4157")),
               person("Ruth", "Ripley", role = "aut"),
               person("Krists", "Boitmanis", role = c( "aut","ctb")),
@@ -36,11 +36,8 @@ Description: The main purpose of this package is to perform simulation-based
    multivariate networks, which can be directed, non-directed, or two-mode;
    and associated actor variables.
    There are also functions for testing parameters and checking goodness of fit.
-   An overview of these models is given in Tom A.B. Snijders (2017), "Stochastic
-   Actor-Oriented Models for Network Dynamics", Annual Review of Statistics and
-   Its Application, 4, 343-363 <doi:10.1146/annurev-statistics-060116-054035>.
-   An extensive manual, scripts, and much further information is at the Siena
-   website <https://www.stats.ox.ac.uk/~snijders/siena/>.
+   An overview of these models is given in Snijders (2017),
+   <doi:10.1146/annurev-statistics-060116-054035>.
 License: GPL-2 | GPL-3 | file LICENSE
 LazyData: yes
 Biarch: yes

NEWS.md

@@ -1,3 +1,41 @@
+# RSiena 1.3.20
+
+## 
+
+2023-04-22
+
+## Changes in RSiena:  
+### Coding:
+### New effects:
+### Corrections:
+  * `updateSpecification` (in `effectsMethods`) now also updates 
+    internal parameter values.
+  * In `TriadCensus`, the empty network will not lead to an error
+    but be reported with the correct triad census.
+  * For `reciAct`, check whether internal parameter ==2 replaced by check whether
+    absolute difference from 2 is less than 0.001.
+  * In `phase2.r`, `z$sd` is calculated using `sqrt(pmax(..., 0))` to avoid the
+    extremely rare case of a negative calculated variance.
+  * In `sienaDataCreate`, handling of structurally determined values 
+    in `checkConstraints` corrected (thanks to issue raised by Jos Elkink).
+### Improvements of functionality:
+  * The keyword `parameter` in `includeInteraction` was dropped because it did
+    not have any consequences. The help page for `includeInteraction` now 
+    explains how internal effect parameters for user-defined interactions
+    are determined.
+  * The column `dimnames` of the `Simulations` array returned by `sienaGOF` 
+    are set to the names of the elements of the auxiliary function.
+  * Standard deviations added to output of `descriptives.sienaGOF`.
+  * Improved error message in `initializeFRAN` in the case of mismatch between
+    effects objects.
+  * Warning in `sienaAlgorithmCreate` if `(maxlike && (!is.null(MaxDegree)))`.
+    This is now also mentioned in the help page for `sienaAlgorithmCreate`.
+### Documentation:
+  * Reference about score-type test added to `Wald.Rd`.
+  * In the help page for `sienaDependent`, it is mentioned that if there are 
+    one-mode as well as two-mode dependent networks,
+    the one-mode networks should come first.
+
 # RSiena 1.3.19
 
 ## 
@@ -18,7 +56,6 @@
     non-missing value, but will stop with an error message 
     if all values are missing.
 
-
 # RSiena 1.3.18
 
 ## 

R/effectsMethods.r

@@ -306,7 +306,7 @@ updateSpecification <- function(effects.to, effects.from, name.to=NULL, name.fro
     effects.to$include[use] <- TRUE
     effects.to$fix[use] <- prevEffects$fix[correspondence][use]
     effects.to$test[use] <- prevEffects$test[correspondence][use]
-    effects.to$parameter[use] <- prevEffects$parameter[correspondence][use]
+    effects.to$parm[use] <- prevEffects$parm[correspondence][use]
     effects.to$randomEffects[use] <- prevEffects$randomEffects[correspondence][use]
     effects.to
 }

R/initializeFRAN.r

@@ -113,7 +113,7 @@ initializeFRAN <- function(z, x, data, effects, prevAns=NULL, initC,
 		{
 			stop("not valid siena data object")
 		}
-		## check the effects object		
+		## check the effects object
 		defaultEffects <- getEffects(data)
 		if (is.null(effects))
 		{
@@ -137,7 +137,9 @@ initializeFRAN <- function(z, x, data, effects, prevAns=NULL, initC,
 				print(userlist[bad])
 				cat("invalid effect requested: see above; \n")
 				cat("there seems to be a mismatch between data set and effects object.\n")
-				stop("Perhaps the effects object must be created from scratch.")
+				cat("This may have been caused by the use of different versions of RSiena")
+				cat("for creating the effects object and now running siena07.\n")
+				stop("Try creating the effects object with the current version of RSiena.")
 			}
 		}
 		if (!inherits(effects, "data.frame"))

R/phase2.r

@@ -50,17 +50,17 @@ phase2.1<- function(z, x, ...)
 #	z$sd <- sqrt(diag(msf))
 # Instead of the preceding line,
 # the following is used for equality with earlier versions.
-    z$sd <- sqrt(apply(z$sf, 2, function(x) sum(x^2) / nrow(z$sf) - mean(x)^2))
+	z$sd <- sqrt(pmax(apply(z$sf, 2, function(x) sum(x^2) / nrow(z$sf) - mean(x)^2),0))
 	if (!z$gmm) 
 	{
    	z$sd[z$fixed] <- 0
 	  z$standardization <-
-						1/sqrt(diag(as.matrix(z$dinvv %*% msf %*% t(z$dinvv))))
+						1/sqrt(pmax(diag(as.matrix(z$dinvv %*% msf %*% t(z$dinvv))),0))
 	}
 	else
 	{
 	  z$sd[which(z$fixed & !z$gmmEffects)] <- 0
-	  z$standardization <- 1/sqrt(diag(as.matrix(z$dinvv %*% msf %*% t(z$dinvv))))
+	  z$standardization <- 1/sqrt(pmax(diag(as.matrix(z$dinvv %*% msf %*% t(z$dinvv))),0))
 	}
     Report(paste("standardization = ", round(z$standardization,4)), cf)
 	if (sum(z$fixed) < z$pp)

R/print07Report.r

@@ -263,12 +263,6 @@ PrintReport <- function(z, x)
 			}
 			else
 			{
-				if (any(z$fixed))
-				{
-					Report(c('(Values of the covariance matrix of estimates\n',
-							' are meaningless for the fixed parameters.)\n\n'),
-						outf)
-				}
 				Report(c("Covariance matrix of estimates",
 						"(correlations below diagonal):\n"), outf)
 				covcor <- z$covtheta

R/sienaDataCreate.r

@@ -1101,8 +1101,12 @@ checkConstraints <- function(z)
 					{
 						var2 <- depvar2[,, obs]
 					}
-					var1[var1 %in% c(10, 11)] <- var1[var1 %in% c(10, 11)] - 10
-					var2[var2 %in% c(10, 11)] <- var2[var2 %in% c(10, 11)] - 10
+#					var1[var1 %in% c(10, 11)] <- var1[var1 %in% c(10, 11)] - 10
+#					var2[var2 %in% c(10, 11)] <- var2[var2 %in% c(10, 11)] - 10
+					var1[var1==10] <- 0
+					var1[var1==11] <- 1
+					var2[var2==10] <- 0
+					var2[var2==11] <- 1
 					## higher
 					if (any(var1 - var2 < 0, na.rm=TRUE))
 					{

R/sienaGOF.r

@@ -133,6 +133,7 @@ sienaGOF <- function(
 					function (i){auxiliaryFunction(i, sienaFitObject$f,
 						sienaFitObject$sims, j, groupName, varName, ...)})
 				simStatsByPeriod <- matrix(simStatsByPeriod, ncol=iterations)
+				dimnames(simStatsByPeriod)[[1]] <-	plotKey
 				dimnames(simStatsByPeriod)[[2]] <-	1:iterations
 				t(simStatsByPeriod)
 				}))
@@ -165,6 +166,7 @@ sienaGOF <- function(
 					flush.console()
 					simStatsByPeriod <-
 							matrix(simStatsByPeriod, ncol=iterations)
+					dimnames(simStatsByPeriod)[[1]] <-	plotKey
 					dimnames(simStatsByPeriod)[[2]] <-	1:iterations
 					t(simStatsByPeriod)
 					})
@@ -862,9 +864,10 @@ descriptives.sienaGOF <- function (x, center=FALSE, scale=FALSE,
 		key <- key[screen]
 	}
 
-	sims.themin <- apply(sims, 2, min)
-	sims.themax <- apply(sims, 2, max)
-	sims.mean <- apply(sims, 2, mean)
+	sims.themin <- apply(sims, 2, min, na.rm=TRUE)
+	sims.themax <- apply(sims, 2, max, na.rm=TRUE)
+	sims.mean <- apply(sims, 2, mean, na.rm=TRUE)
+	sims.sd <- apply(sims, 2, sd, na.rm=TRUE)
 	sims.min <- pmin(sims.themin, obs)
 	sims.max <- pmax(sims.themax, obs)
 
@@ -889,6 +892,7 @@ descriptives.sienaGOF <- function (x, center=FALSE, scale=FALSE,
 		sims.mean <- sims.mean/sims.range
 		sims.min <- sims.min/sims.range
 		sims.max <- sims.max/sims.range
+		sims.sd <- sims.sd/sims.range
 	}
 
 	screen <- sapply(1:ncol(obs),function(i){
@@ -913,18 +917,19 @@ descriptives.sienaGOF <- function (x, center=FALSE, scale=FALSE,
 				sort(sims[,i])[ind.upper]  )
 	ypg <- sapply(1:ncol(sims), function(i)	mean(sims[,i] > obs[1,i]))
 	ypp <- sapply(1:ncol(sims), function(i)	mean(sims[,i] >= obs[1,i]))
-    violins <- matrix(NA, 9, ncol(sims))
+    violins <- matrix(NA, 10, ncol(sims))
 	violins[1,] <- sims.themax
 	violins[2,] <- yperc.upper
 	violins[3,] <- sims.mean
 	violins[4,] <- yperc.mid
 	violins[5,] <- yperc.lower
 	violins[6,] <- sims.themin
-	violins[7,] <- obs
-    violins[8, ] <- ypg
-    violins[9, ] <- ypp
+	violins[7,] <- sims.sd
+	violins[8,] <- obs
+    violins[9, ] <- ypg
+    violins[10, ] <- ypp
     rownames(violins) <- c("max", "perc.upper", "mean", "median",
-        "perc.lower", "min", "obs", "p>", "p>=")
+        "perc.lower", "min", "sd", "obs", "p>", "p>=")
 	colnames(violins) <- key
 	violins
 }
@@ -1414,25 +1419,27 @@ TriadCensus <- function (i, obsData, sims, period, groupName, varName, levls = 1
           "300"  = 0)
 
   # iterate through all non-empty dyads (from lower to higher ID)
-  for(i in 1:N){
-    for(j in neighborsHigher[[i]]){
-      # set of nodes that are linked to i and j
-      third <- setdiff( union(neighbors[[i]], neighbors[[j]]),
-                    c(i, j) )
-      # store triads with just one tie
-      triadType <- ifelse(matReciprocal[i,j] == 2, 3, 2)
-      tc[triadType] <- tc[triadType] + N - length(third) - 2
-      for (k in third){
+  if (length(neighborsHigher) > 0){ # else mat is the zero matrix
+	for(ii in 1:N){
+		for(j in neighborsHigher[[ii]]){
+      # set of nodes that are linked to ii and j
+		third <- setdiff( union(neighbors[[ii]], neighbors[[j]]),
+                    c(ii, j) )
+		# store triads with just one tie
+		triadType <- ifelse(matReciprocal[ii,j] == 2, 3, 2)
+		tc[triadType] <- tc[triadType] + N - length(third) - 2
+		for (k in third){
         # only store triads once
-        if(j < k || ( i < k && k < j && !(k %in% neighbors[[i]]) ) ){
-          t1 <- matDirected[i,j]
-          t2 <- matDirected[j,k]
-          t3 <- matDirected[i,k]
-          triadType <- lookup[t1, t2, t3]
-          tc[triadType] <- tc[triadType] + 1
-        }
-      }
-    }
+			if(j < k || ( ii < k && k < j && !(k %in% neighbors[[ii]]) ) ){
+				t1 <- matDirected[ii,j]
+				t2 <- matDirected[j,k]
+				t3 <- matDirected[ii,k]
+				triadType <- lookup[t1, t2, t3]
+				tc[triadType] <- tc[triadType] + 1
+				}
+			}
+		}
+	}
   }
   # assign residual to empty triad count
   tc[1] <- 1/6 * N*(N-1)*(N-2) - sum(tc[2:16])

R/sienaModelCreate.r

@@ -28,14 +28,17 @@ sienaModelCreate <- function(fn,
 {
 	model <- NULL
 	checking <- any(grepl("_R_CHECK", names(Sys.getenv())))
+	if (maxlike && (!is.null(MaxDegree)))
+	{
+		warning("maxlike and MaxDegree are incompatible")
+	}	
 	if (is.null(projname) | checking)
 	{
 		model$projname <- tempfile("Siena")
 		if (checking)
 		{
 	cat('If you use this algorithm object, siena07 will create/use an output file',
 				paste('Siena','.txt',sep=''),'.\n')
-
 		}
 		else
 		{

R/sienaTimeTest.r

@@ -1180,7 +1180,6 @@ sienaTimeFix <- function(effects, data=NULL, getDocumentation=FALSE)
 								interaction1= c(effect$interaction1,
 									dname),
 								interaction2=effect$interaction2,
-								parameter=effect$parameter,
 								name=depvar, verbose=FALSE)
 					} else {
 						newEffects <-
@@ -1191,7 +1190,6 @@ sienaTimeFix <- function(effects, data=NULL, getDocumentation=FALSE)
 								interaction1= c(effect$interaction1,
 									dname),
 								interaction2=effect$interaction2,
-								parameter=effect$parameter,
 								name=depvar, verbose=FALSE)
 					}
 					## find the row altered

R/sienaeffects.r

@@ -93,7 +93,7 @@ includeEffects <- function(myeff, ..., include=TRUE, name=myeff$name[1],
 includeInteraction <- function(myeff, ...,
 				include=TRUE, name=myeff$name[1],
 				type="eval", interaction1=rep("", 3), interaction2=rep("", 3),
-				parameter=NULL, fix=FALSE, test=FALSE, random=FALSE,
+				fix=FALSE, test=FALSE, random=FALSE,
 				initialValue=0,
 				character=FALSE, verbose=TRUE)
 {
@@ -228,7 +228,6 @@ includeInteraction <- function(myeff, ...,
 			c(effect1, effect2, effect3)
 		myeff[intn, "fix"] <- fix
 		myeff[intn, "test"] <- test
-		if (!is.null(parameter)) {myeff[intn, "parm"] <- parameter}
 		myeff[intn, "randomEffects"] <- random
 		myeff[intn, "initialValue"] <- initialValue
 	}

README.md

@@ -100,3 +100,25 @@ please install the `{remotes}` package from CRAN and then enter into the console
 The development version of RSiena can be similarly installed as:
 
 `remotes::install_github("snlab-nl/rsiena@develop")`
+
+## Citation
+
+To cite the RSiena package in publications use:
+
+> Ruth M. Ripley, Tom A. B. Snijders, Zsofia Boda, Andras Voros, and Paulina Preciado(2023). Manual
+> for Siena version 4.0. R package version 1.3.19.
+> https://www.cran.r-project.org/web/packages/RSiena/.
+
+A BibTeX entry for LaTeX users is
+
+```bib
+@TechReport{,
+  title = {Manual for {Siena} version 4.0},
+  author = {Ruth M. Ripley and Tom A. B. Snijders and Zsofia B'{o}da and Andr'{a}s V"{o}r"{o}s and Paulina Preciado},
+  year = {2023},
+  institution = {Oxford: University of Oxford, Department of Statistics; Nuffield College},
+  note = {R package version 1.3.19. https://www.cran.r-project.org/web/packages/RSiena/},
+}
+```
+
+For more references, see https://www.stats.ox.ac.uk/~snijders/siena/.