\documentclass[12pt,a4paper]{article}
\usepackage[pdftex,dvipsnames]{color}
\usepackage[noend]{algorithmic}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage[pdfstartview={}]{hyperref}
\newcommand{\eps}{\varepsilon}
\newcommand{\abso}[1]{\;\mid#1\mid\;}
\renewcommand{\=}{\,=\,}
\newcommand{\+}{\,+\,}
% ----------------------------------------------------------------
\newcommand{\remark}[1]{\par\noindent{\color[named]{ProcessBlue}#1}\par}
\newcommand{\mcc}[2]{\multicolumn{#1}{c}{#2}}
\newcommand{\mcp}[2]{\multicolumn{#1}{c|}{#2}}
\newcommand{\nm}[1]{\textsf{ #1}}
\newcommand{\nnm}[1]{\textsf{\textit{#1}}}
\newcommand{\nmm}[1]{\nnm{#1}}
\newcommand{\ts}[1]{\par{\color[named]{Red}TS: #1}\par}

% no labels in list of references:
\makeatletter
\renewcommand\@biblabel{}
\makeatother

\hyphenation{Snij-ders Duijn DataSpecification dataspecification dependentvariable ModelSpecification}

% centered section headings with a period after the number;
% sans serif fonts for section and subsection headings
\renewcommand{\thesection}{\arabic{section}.}
\renewcommand{\thesubsection}{\thesection\arabic{subsection}}
\makeatletter
 \renewcommand{\section}{\@startsection{section}{1}
                {0pt}{\baselineskip}{0.5\baselineskip}
                {\centering\sffamily} }
 \renewcommand{\subsection}{\@startsection{subsection}{2}
                {0pt}{0.7\baselineskip}{0.3\baselineskip}
                {\sffamily} }
 \renewcommand{\subsubsection}{\@startsection{subsubsection}{3}
                {0pt}{0.5\baselineskip}{0.3\baselineskip}
                {\it\sffamily} }
\makeatother


\renewcommand{\baselinestretch}{1.0} %% For line spacing.
\setlength{\parindent}{0pt}
\setlength{\parskip}{1ex plus1ex}
\raggedright
\newcommand{\sfn}[1]{\textbf{\texttt{#1}}}
\newcommand{\R}{{\sf R }}
\newcommand{\Rn}{{\sf R}}
\newcommand{\rs}{{\sf RSiena}}
\newcommand{\RS}{{\sf RSiena }}
\newcommand{\SI}{{\sf Siena3 }}
\newcommand{\Sn}{{\sf Siena3}}
\begin{document}
\algsetup{indent=2em}

\title{RSiena}
\author{Ruth M.\ Ripley\footnote{University of Oxford}\\
modified by Tom A.B.\ Snijders\footnote{Universities of Oxford and Groningen}}
\maketitle
\date{}


\section{Grand overview of  RSiena}

This overview of \RS sometimes refers to the earlier version \SI to
facilitate understanding by those who have a knowledge
of the code of \Sn. Readers who have no such knowledge
can happily ignore these references.

Two special features of the \RS package are that it uses C++ code for
time-consuming parts, and that a wrapper is available to use the package
seemingly independently of \Rn.

The script file \nnm{sienascript} starts up \R in such a way that it launches a
tcl/tk graphical user interface (gui) resembling the Stocnet interface for \Sn,
and controls all processing thereafter. This interface is accessible on Windows
via the command \nnm{siena01Gui}. (Prior to R 2.12.0 the funciton
\sfn{siena.exe} was provided as an equivalent to \sfn{sienascript}: it has been
removed as it was rarely used and proved difficult to maintain.)

The high-level functions called by the gui, such as \nnm{siena07} described
below, are also accessible within \R with the usual \Rn-type user interface
along the lines of model-fitting functions such as \verb|lm()|. (A formula
interface is still on the wish list, but already now it is relatively
straightforward to use the package without the gui.)

The \R functions call C++ only where speed is critical. From my profiling of
\nnm{siena07}, the estimation function, I think that only the simulation
function (which simply performs one simulation of the complete model for a given
set of parameters and returns the statistics from the simulated networks) needs
to be in C++. The bulk of the time is spent in calculating the contributions to
the effects when simulating.

In \Sn, this simulation function is \nm{FRAN} which, for method of moments
estimation, simply calls \nm{simstats}: in \rs, it is the C++ function
\nnm{model} called by one of the \R functions \nnm{simstats0c} or
\nnm{maxlikec}, which are the two candidates currently available for the element
of the \RS model object named \nnm{FRAN}. In this document I have used the name
\nm{FRAN} to refer to the simulation routine.

\nnm{siena07} is intended to be written in such a way that different simulation
functions could be used within the same Robbins Monro algorithm. In practice the
separation is not quite complete, but it is nearly so.

It would be feasible within a C++ \nm{FRAN} to call \R functions for some
effects or functions if desired, to facilitate adding new ones, although they
would be slow to run!

We use functions from the C part of \R, to provide random numbers within the
C++.1 \nnm{simstats0c} and \nnm{maxlikec} have three types of calls: a initial
one which calls various C++ routines to setup the data, a final one which sets
the C++ data pointers to null to clean up the C++ memory, and multiple
``normal'' ones which call the function \nnm{model} to perform one complete
simulation. (In this it does not correspond exactly to \nm{simstats} in Siena3!)

With this design, we have introduced parallel processing by using multiple \R
processes. In \nnm{simstats0c} we run some of the simulations in each process:
this was trivial to introduce into Phase 1 and Phase 3, but to use it in Phase2
we have altered the algorithm to use the average of more than one simulation at
a time.

For \nnm{maxlikec} we use different processes for each wave.
This is because the chains are carried from one simulation to the next,
and organizing the parallel processes by simulations (update steps of the Robbins
Monro algorithm) would require too much passing of information.

We use the
\R package \verb|parallel| to create and control the multiple processes, and
to provide multiple random number streams.  The term
`cluster' below refers to the cluster of multiple processors. (I vary between
using \emph{processors} and \emph{processes} as it is possible to run \RS with
multiple \emph{processes} on a machine with only one processor.)


\section{Data types}

\RS provides various classes of data objects, designed to interface with the
functions \nnm{robmon} and \nm{simstats}. A brief list:

\begin{description}
\item[\nnm{siena}] Data for a single project
\item[\nnm{sienaGroup}] A list of \nnm{siena} objects, with global attributes,
  used for multi-group projects
\item[\nnm{sienaEffects}] Data frame of effects.
\item[\nnm{sienaGroupEffects}] Data frame of effects for a group object.
\item[\nnm{sienaModel}] Contains the fitting options.
\item[\nnm{sienaNodeSet}] Actor set, used to distinguish nodes in data sets with
  multiple or two-mode networks.
\item[\nnm{sienaDependent}] A single dependent variable, (i.e.\ network or behavior
  variable, all waves)
\item[\nnm{coCovar}] Constant covariate
\item[\nnm{coDyadCovar}] Constant dyadic covariate
\item[\nnm{varCovar}] Varying covariate
\item[\nnm{varDyadCovar}] Varying dyadic covariate
\item[\nnm{sienaCompositionChange}] List of changes, entry for each node.
\item[\nnm{sienaFit}] Currently contains (almost) everything from the
  estimation.
\end{description}

The structure of each is documented in the corresponding \R help file:
\sfn{?classname}.

Data objects of class \nnm{siena} are created by the function
\nnm{sienaDataCreate}.

Effect objects of class \nnm{sienaEffects} are created by the function
\nnm{getEffects}.

The creation functions can be called directly by the user or from the Gui or via
\nnm{sienaDataCreateFromSession}, depending on whether the data is alrady in \R
objects or still on files.

The function \nnm{robmon} requires a \nnm{sienaModel} object
as an argument. One element of this object is named \nm{FRAN} and contains the
name of the required simulation function.

The functions\nnm{simstats0c} or \nnm{maxlikec}, used as an instance of
\nm{FRAN}, require a \nnm{siena} (or \nnm{sienaGroup}) object and a
\nnm{sienaEffects} (or \nnm{sienaGroupEffects}) object as arguments.

\section{sienaDataCreate}

This function has only one named argument: a list of actor (node) sets. The
default is a single set of the required size named \sfn{Actors}. All other
arguments are unnamed and correspond to networks, covariates or composition
change files. The objects are validated and have various attributes
added. For covariates the attributes are added using a \emph{method} for their
class.

\begin{algorithmic}
\STATE Check that objects have names, using the object name if none is given in
the function call.
\IF{no objects}
\STATE stop
\ENDIF
\IF{any duplicate names}
\STATE stop
\ENDIF
\STATE create a list of each type of object, checking that all dependent
variables have the same number of observations. (Stop if not).
\IF{no dependent variables}
\STATE stop
\ENDIF
\IF{no node set argument}
\STATE create a list of nodesets containing a single
nodeset named Actors, with the number of nodes equal to the number of senders of
the first dependent variable
\ENDIF
\FORALL{covariates}
\STATE Appropriate validation and processing (see below).
\ENDFOR
\STATE Process any composition change objects (see section
\ref{sec:compositionChange}) \\
\STATE Process the dependent variables (see section \ref{sec:dependent})
\STATE Check constraints if there are multiple networks. (section
\ref{sec:constraints}).
\STATE Calculate similarity means for alters for each covariate and dependent
network (see section \ref{sec:covarDist2}): dropped in version 1.1-285.
\end{algorithmic}
\subsection{ Covariates}
\subsubsection{Constant Covariate}
\begin{algorithmic}
\STATE Check the nodeset (section \ref{sec:nodeset})
\STATE Create attributes: (a class method)\\
\sfn{ }\nnm{mean} ignore missings\\
\sfn{ }\nnm{range} Extent of range, ignore missings. Make sure is a double.\\
\sfn{ }\nnm{range2} Ends of range, ignore missings\\
\sfn{ }\nnm{moreThan2} TRUE if more than 2 distinct values, ignoring missings\\
\sfn{ }\nnm{vartotal} sum of non-missing values\\
\sfn{ }\nnm{poszvar} TRUE if more than 1 distinct value in the centered values
or any missing\\
\sfn{ }\nnm{simMean} See section \ref{sec:similarityMean}\\
\sfn{ }\nnm{nonMissingCount} count of non missing values\\
\sfn{ }\nnm{name} name of object
\STATE Subtract the mean from the values
\end{algorithmic}
\subsubsection{Changing covariate}
\begin{algorithmic}
\IF{less than 3 waves}
\STATE Stop: changing covariate inappropriate (to reduce confusion among users!)
\ENDIF
\STATE Check the nodeset (section \ref{sec:nodeset})
\IF{less than (number of waves - 1) columns}
\STATE Stop: not enough values
\ENDIF
\IF{more than (number of waves - 1) columns}
\STATE remove the excess, carefully preserving the attributes apart from the
dimensions.
\ENDIF
\STATE Create attributes: (a class method)\\
\sfn{ }\nnm{mean} ignore missings\\
\sfn{ }\nnm{meanp}  mean for each wave, ignore missings\\
\sfn{ }\nnm{range} Extent of range, ignore missings. Make sure a double.\\
\sfn{ }\nnm{rangep} Extent of range for each wave, ignore missings.
Make sure is a double.\\
\sfn{ }\nnm{range2} Ends of range, ignore missings\\
\sfn{ }\nnm{moreThan2} TRUE if more than 2 distinct values, ignoring missings\\
\sfn{ }\nnm{vartotal} sum of non-missing values\\
\sfn{ }\nnm{poszvar} TRUE if more than 1 distinct value in the centered values
or any missing\\
\sfn{ }\nnm{simMean} See section \ref{sec:similarityMean}\\
\sfn{ }\nnm{nonMissingCount} count of non missing values\\
\sfn{ }\nnm{name} name of object
\STATE Subtract the mean from the values
\end{algorithmic}
\subsubsection{Constant dyadic covariate}
\begin{algorithmic}
\STATE Check the nodesets (section \ref{sec:nodeset})
\IF{attribute type is oneMode}
\STATE set diagonal to \emph{missing} so is ignored in mean and range
\ENDIF
\STATE Create attributes: (a class method) \\
\sfn{ }\nnm{mean} ignore missings\\
\sfn{ }\nnm{range} Extent of range, ignore missings. Make sure a double.\\
\sfn{ }\nnm{range2} Ends of range, ignore missings\\
\sfn{ }\nnm{name} name of object
\IF{attribute type is oneMode}
\STATE set diagonal to \emph{zero}
\ENDIF
\end{algorithmic}
\subsubsection{Changing dyadic covariate}
\begin{algorithmic}
\IF{less than 3 waves}
\STATE Stop: changing covariate inappropriate (to reduce confusion among users!)
\ENDIF
\STATE Check the nodesets (section \ref{sec:nodeset})
\IF{less than (number of waves - 1) columns}
\STATE Stop: not enough values
\ENDIF
\IF{more than (number of waves - 1) columns}
\STATE remove the excess, carefully preserving the attributes apart from the
dimensions.
\ENDIF
\IF{attribute type is oneMode}
\STATE set all diagonals to \emph{missing} so are ignored in mean and range
\ENDIF
\STATE Create attributes: (a class method) \\
\sfn{ }\nnm{mean} ignore missings\\
\sfn{ }\nnm{range} Extent of range, ignore missings. Make sure a double.\\
\sfn{ }\nnm{name} name of object
\IF{attribute type is oneMode}
\STATE set all diagonals to \emph{zero}
\ENDIF
\end{algorithmic}
\subsection{Composition change objects}
\label{sec:compositionChange}
\begin{algorithmic}
\STATE Check there are no duplicates in the nodesets: only one change object
per nodeset is allowed.
\FORALL{composition change objects}
\STATE Check the nodeset (section \ref{sec:nodeset})
\STATE Check that the ends of each interval in each object are not less than 1
or greater than the number of waves and that each line has an even number of
digits.
\STATE Generate a data frame of events(section \ref{sec:events}),
a matrix of activeStart flags(section \ref{sec:active}) and a matrix
of actions(section \ref{sec:actions})\\
\STATE Add these to the object as attributes
\ENDFOR
\end{algorithmic}

\subsubsection{Events}
\label{sec:events}
\begin{algorithmic}
\STATE Data frame with columns: \\
\sfn{ }\nnm{event} ``join'' or ``leave'' (a factor)\\
\sfn{ }\nnm{period}\\
\sfn{ }\nnm{actor}\\
\sfn{ }\nnm{time} between 0 and 1
\end{algorithmic}
\subsubsection{ActiveStart Flags}
\label{sec:active}
\begin{algorithmic}
\STATE activeStart matrix has a row per actor and a column per period\\
\sfn{ }TRUE if the actor is active at the start of the period, otherwise FALSE
\end{algorithmic}
\subsubsection{Action}
\label{sec:actions}
\begin{algorithmic}
\STATE Action matrix is same shape as Active flag matrix, with entries \\
\sfn{ }0 Active at start\\
\sfn{ }1 Inactive at start, never previously active\\
\sfn{ }2 Inactive at start, previously active but never active again\\
\sfn{ }3 Inactive at start, previously active and active again\\
\end{algorithmic}
\subsection{Check NodeSet}
\label{sec:nodeset}
\begin{algorithmic}
\IF {Nodeset name in the list and lengths match}
\STATE Valid
\ELSE
\STATE Invalid
\ENDIF
\end{algorithmic}
\subsection{Dependent variables}
\label{sec:dependent}
NB The attributes list tends to change rather quickly and some items may no
longer be required. Some should be used in the C, but are not...
\begin{algorithmic}
\STATE Validate the nodeset(s)
\STATE Create an attribute \nnm{name} with name of the object
\IF{behavior variable}
\STATE Create attributes:\\
\sfn{ }\nnm{distance} sum of absolute differences by period, ignoring missings\\
\sfn{ }\nnm{vals} table of values by period, NA included as a value\\
\sfn{ }\nnm{nval} vector of non-missing counts by period\\
\sfn{ }\nnm{noMissing} vector of missing counts be period\\
\sfn{ }\nnm{range} overall range\\
\sfn{ }\nnm{range2} overall min and max\\
\sfn{ }\nnm{moreThan2} TRUE if number of distinct values more than 2 (includes
missings as a value) ??? inconsistent with covariates?\\
\sfn{ }\nnm{poszvar} TRUE if more than one distinct value or any missing
values.\\
\sfn{ }\nnm{modes} vector of modes of rounded values per period. Might give
multiple results?\\
\sfn{ }\nnm{missing} TRUE if any missing\\
\sfn{ }\nnm{simMean} value of similarity mean (see section
\ref{sec:similarityMean})\\
\sfn{ }\nnm{structural} FALSE Not allowed!\\
\sfn{ }\nnm{balmean} NA\\
\sfn{ }\nnm{structMean} NA\\
\sfn{ }\nnm{uponly} TRUE if all changes increase, ignoring missings\\
\sfn{ }\nnm{downonly} TRUE if all changes decrease, ignoring missings
\ELSE[bipartite or onemode]
\STATE create attributes:\\
\sfn{ }\nnm{distance} count of changes, ignoring missing and structural values
and diagonals if not bipartite.\\
\sfn{ }\nnm{uponly} TRUE if ties are only ever created, never lost\\
\sfn{ }\nnm{downonly} TRUE if ties are only ever lost, never created\\
\IF{one-mode}
\STATE Create attributes:\\
\sfn{ }\nnm{balMean} see section(\ref{sec:balmean})\\
\sfn{ }\nnm{structMean} see section(\ref{sec:structmean})\\
\sfn{ }\nnm{symmetric} TRUE if all waves are symmetric\\
\sfn{ }\nnm{missing} TRUE if any missing values (except on diagonal)\\
\sfn{ }\nnm{structural} TRUE if any 10 or 11\\
\sfn{ }\nnm{vals} table of counts of values by period\\
\sfn{ }\nnm{nval} Counts of non-missing values by period, excluding diagonal\\
\sfn{ }\nnm{range2} Min and max of non-structural values\\
\sfn{ }\nnm{noMissing} Number of missing values by period\\
\sfn{ }\nnm{noMissingEither} Number of missing values at start or finish of
period (excludes final).\\
\sfn{ }\nnm{nonMissingEither} Number of non missing values at start or finish of
period (excludes final).\\
\sfn{ }\nnm{simMean} NA\\
\sfn{ }\nnm{ones} Count of values equal to 1 by period\\
\sfn{ }\nnm{density} Density of network by period\\
\sfn{ }\nnm{degree} Average degree by period\\
\sfn{ }\nnm{averageOutDegree} overall average degree\\
\sfn{ }\nnm{averageInDegree} overall average degree\\
\sfn{ }\nnm{maxObsOutDegree} Maximum observed outdegree by period\\
\sfn{ }\nnm{missings} count of missings by period\\
\ELSIF{bipartite}
\STATE Create attributes:\\
\sfn{ }\nnm{balMean} NA\\
\sfn{ }\nnm{structMean} NA\\
\sfn{ }\nnm{symmetric} FALSE\\
\sfn{ }\nnm{missing} TRUE if any missing values\\
\sfn{ }\nnm{structural} TRUE if any 10 or 11\\
\sfn{ }\nnm{vals} table of counts of values by period\\
\sfn{ }\nnm{nval} Counts of non-missing values by period\\
\sfn{ }\nnm{range2} Min and max of non-structural values\\
\sfn{ }\nnm{noMissing} Number of missing values by period\\
\sfn{ }\nnm{noMissingEither} Number of missing values at start or finish of
period (excludes final).\\
\sfn{ }\nnm{nonMissingEither} Number of non missing values at start or finish of
period (excludes final).\\
\sfn{ }\nnm{simMean} NA\\
\sfn{ }\nnm{ones} Count of values equal to 1 by period\\
\sfn{ }\nnm{density} Density of network by period\\
\sfn{ }\nnm{degree} Average degree by period\\
\sfn{ }\nnm{averageOutDegree} overall average degree\\
\sfn{ }\nnm{averageInDegree} overall average degree\\
\sfn{ }\nnm{missings} count of missings by period\\
\ENDIF
\ENDIF
\end{algorithmic}
\subsection{Similarity mean}
\label{sec:similarityMean}
\begin{algorithmic}
\FORALL[waves]{columns of matrix}
\FORALL[actors]{entries in column}
\STATE in a copy of the column set this entry to NA
\STATE Calculate 1 - abs(this entry - copy column)/ range
\STATE Sum this over the nonmissing entries in this vector
\STATE Count the nonmissing entries in this vector
\ENDFOR
\ENDFOR
\STATE Sum nonmissing entries and counts over the columns
\STATE Calculate the similarity mean as this total sum divided by total count.
\STATE For possible later use, also return the total sum and total count.
\end{algorithmic}
\subsection{Balance mean}
\label{sec:balmean}
\begin{algorithmic}
\STATE In calculations, remove diagonal and replace structural values by the
values represented.\\
\STATE Numerator = sum over all columns of\\
2 * count of non-zero entries * count of non-missing non-nonzero entries\\
\STATE Denominator = sum over all columns of \\
count of non-missing entries times one less than this.
\STATE Mean is numerator divided by denominator.
\end{algorithmic}
\subsection{Structural  mean}
\label{sec:structmean}
\begin{algorithmic}
\STATE In calculations, remove diagonal and replace structural values by the
values represented.\\
\STATE Numerator = sum over all rows of\\
2 * count of non-zero entries * count of non-missing non-nonzero entries\\
\STATE Denominator = sum over all rows of \\
count of non-missing entries times one less than this.
\STATE Mean is numerator divided by denominator.
\end{algorithmic}
\subsection{Constraints between networks}
\label{sec:constraints}
\begin{algorithmic}
\STATE Make a two column matrix containing the names of all possible pairs of
dependent variables, including pairs with themselves.
\STATE Identify any dependent variables that can relate: same type and have the
same node sets.
\STATE Create a list of these possibly relating dependent variables for each of
access later.
\FORALL{row in the matrix of pairs where the two columns are not the same}
\IF{nodeset(s) match and type matches and not behavior and either both networks
 are symmetric or neither is}
\STATE In checks, replaces structurals by 0/1 first, and ignore missings
\STATE Check \emph{higher}: first network always greater than or equal to second
network.
\STATE Check \emph{disjoint}: sum of product of two networks not greater than 0.
\STATE Check \emph{atLeastOne}: sum of two networks never equal 0.
\ENDIF
\ENDFOR
\end{algorithmic}
\subsection{Similarity means at distance 2}
\label{sec:covarDist2}
\begin{algorithmic}
\FORALL{constant covariates, varying covariates}
\FORALL{dependent networks which have the node set of the constant covariate as
  their receivers (not behavior variables)}
\STATE Calculate the corresponding similarity mean: see section
\ref{sec:simdist2}).
\ENDFOR
\ENDFOR
\FORALL{behavior variables}
\FORALL{dependent networks which have the node set of the constant covariate as
  their receivers (not behavior variables)}
\STATE Substract the mean from the behavior variable
\STATE Find the range of the behavior variable
\STATE Calculate the corresponding similarity mean using the centered behavior
variable values (omitting final column) and the calculated range: see section
\ref{sec:simdist2}).
(i.e.\ centering and range are done on complete variable)
\ENDFOR
\ENDFOR
\end{algorithmic}
\subsubsection{Alter similarity calculation}
\label{sec:simdist2}
\begin{algorithmic}
\FORALL{observations except the last}
\STATE Replace structurals by 0/1
\FORALL{rows of network matrix}
\IF{sum of nonmissing entries is 0}
\STATE Set vi for row to 0
\ELSIF{all covariate values corresponding to non zero network entries are
  missing}
\STATE Set vi for row to NA
\ELSE
\STATE Set vi for row to sum of covariate times network row divided by the sum
of the network row, ignoring missings in both cases.
\ENDIF
\ENDFOR
\STATE Call rangeAndSimilarity using vi and the range if passed in (behavior
variables) to obtain the values \nnm{simTotal} and \nnm{simCnt} for this
overvation.
\STATE Accumulate these two values
\ENDFOR
Divide sum of \nnm{simTotal} by sum of \nnm{simCnt} over observations (excluding
the final one).
\end{algorithmic}
\section{getEffects}

This function generates an effects data object corresponding to a Siena Data
object or a Siena Group Object.

In general, effects are driven by selecting rows in the \sfn{allEffects} data
frame for some effect group and then substituting variable names into the spaces
marked by \sfn{xxxxxx} and similar.

For a group object, the effects are created using the first data object plus the
total number of observations. Attributes are first copied from the group level
to the first data object. The only changes required are to fill in the starting
values for the rate effects for the later objects and to adjust the starting
values for density, reciprocity, linear effects.

The function \sfn{networkRateEffects} creates the required number of rate
effects for networks.  \sfn{createEffects} extracts the rows from the effects
data frame for an effect group and calls the function \sfn{substituteNames} to
replace the variable fields by the current variable name. It now creates the
complete effect rows including endowment effect copies.

\subsection{Siena Data Object}
\begin{algorithmic}
\FORALL{dependent variables}
\STATE Create appropriate effects
\STATE Set \nnm{netType} to the value \nnm{oneMode},
      \nnm{bipartite},  \nnm{behavior}, or \nnm{continuous}.
\ENDFOR
\end{algorithmic}
\subsection{OneMode Network Effects}
\begin{algorithmic}
  \STATE Call \sfn{networkRateEffects} to get the rate effects
  \STATE Use \sfn{symmetricObjective} or \sfn{nonSymmetricObjective} effect
  groups to create the basic objective function effects.
  \FORALL{dyadic covariates with first node set matching}
  \STATE Use \sfn{dyadObjective} effect group to add appropriate objective
  function effects
  \ENDFOR
  \FORALL{constant covariates, behavior variables, changing covariates with the
    same node set}
  \STATE Call function \sfn{covarOneModeEff} to add the appropriate effects.\\
  Note \nnm{poszvar} is always TRUE for behavior variables.
  \ENDFOR
  \IF{any covariates or behavior variables}
  \STATE Add \nnm{nintn} rows for user specified interaction effects
  \ENDIF
  \FORALL{distinct dependent network variables with the same node set}
  \IF{oneMode}
  \STATE Use \sfn{nonSymmetricSymmetricObjective} or
  \sfn{nonSymmetricNonSymmetricObjective} effect groups to add appropriate
    effects
  \STATE Use \sfn{tripleNetworkObjective} effect group to add appropriate
    effects for pairs of other dependent networks
    (`other' meaning that they have the role of explanatory variables)
    that either are both one-mode,
    or both are bipartite with the same second node sets
  \ELSIF{bipartite which matches on nodeset 1}
  \STATE Use \sfn{nonSymmetricBipartiteObjective} effect group to add
  appropriate effects
  \ENDIF
  \FORALL{actor covariates or behavior variables with the same node set}
  \STATE Use \sfn{covarNetNetObjective} effect group to add appropriate effects.
  \ENDFOR
  \ENDFOR
  \IF{more than one network}
  \STATE paste the network name at the front of all the objective function
  effects
  \ENDIF
  \STATE Alter the text for endowment effects to start "Lost ties:"
  \STATE Calculate the starting values for the default effects
  (see section \ref{sec:onemodestart})
  \STATE Select the default rate effects by setting \nnm{include} to TRUE for
  the \sfn{basic rate} effects.
  \STATE Add the starting value for the rate to the \nnm{initialValue} column
  of the basic rate effects
  \IF {symmetric}
  \STATE Set \nnm{include} to TRUE for the degree (density) evaluation effect
  and transitive triads evaluation effect.
  \ELSE
  \IF{no period is uponly or downonly}
  \STATE  Set \nnm{include} to TRUE for the degree (density) evaluation effect
  \STATE Add the starting values for the degree(density) evaluation effect
  calculated by \sfn{getNetworkStartingVals} to the \nnm{initialValue}.
  \ELSE
  \STATE Remove both degree (density) effects from the data frame.
  \ENDIF
  \STATE Set \nnm{include} to TRUE for the reciprocity evaluation effect.
  \ENDIF
\end{algorithmic}
\subsection{Behavior Variable Effects}
\begin{algorithmic}
\STATE Use \sfn{behaviorRate} effect group to get the rate effects. Either
remove the second one or duplicate the second and remove the first to match the
number of observations.
\STATE Use \sfn{behaviorObjective} effect group to create the basic objective
function effects.
\FORALL{other dependent variables which match on first node set}
\STATE Use \sfn{behaviorOneModeObjective} or \sfn{behaviorBipartiteObjective} to
add the objective function effects with respect to this network.
\STATE Use \sfn{behaviorOneModeRate} or \sfn{behaviorBipartiteRate} to add
the rate effects with respect to this network.
\ENDFOR
\FORALL{constant covariates, other behavior variables or changing covariates}
\STATE Call \sfn{covBehEff} and \sfn{covBBehEff}to add the interaction effects
\ENDFOR
\FORALL{networks with same node set (first for bipartite)}
\STATE Use \sfn{behaviorOneModeObjective2} or \sfn{behaviorBipartiteObjective2}
effect group to create a second set of objective function effects.
\ENDFOR
\STATE Add  \nnm{behNintn} unspecified behavior interaction effects
  \STATE Create the effects data frame by calling \sfn{createObjEffectList} and
  \sfn{createRateEffectList}.  This creates e.g.\ the evaluation and endowment
  effect copies.
  \STATE Select the default effects by setting \nnm{include} to TRUE for
  \sfn{basic rate} and \sfn{linear shape}
  (if not any period uponly or downonly)
  and \sfn{quadratic shape}
  (if the range of the variable is greater than or equal to 2)
  evaluation effects.
  \IF {any period uponly or downonly}
  \STATE remove the linear effects (evaluation and endowment)
  from the data frame.
  \ENDIF
  \STATE Add the starting values for the default effects calculated by
  \sfn{getBehaviorStartingVals}  (see section \ref{sec:behaviorstart})
 to the \nnm{initialValue} column of the data
  frame.
\STATE Alter the text for endowment effects to start with "dec. beh."
\end{algorithmic}
\subsection{Bipartite Network Effects}
\begin{algorithmic}
  \STATE Call \sfn{networkRateEffects} to get the rate effects
  \STATE Use \sfn{bipartiteObjective} effect group to create the objective
  function effects.
  \FORALL{dyadic covariates with both node sets matching}
  \STATE Use \sfn{dyadObjective} effect group to create the appropriate effects
  \ENDFOR
  \FORALL{constant covariates, behavior variables, changing covariates}
  \STATE Call function \sfn{covarBipartiteEff} to add the appropriate effects\\
  poszvar is always TRUE for behavior variables.
  \ENDFOR
  \IF{any covariates or behavior variables}
  \STATE Add \nnm{nintn} rows for user specified interaction effects
  \ENDIF
  \FORALL{distinct dependent network variables with the same node set}
  \IF{oneMode}
  \STATE Use \sfn{bipartiteSymmetricObjective} or
  \sfn{bipartiteNonSymmetricObjective} effect groups to add the appropriate
the effects
  \ELSIF{bipartite and matches first node set)}
  \STATE Use \sfn{bipartiteBipartiteObjective} effect group to add the
  appropriate  effects
  \ENDIF
  \STATE NB no \sfn{covarNetNetObjective} here?
  \ENDFOR
  \IF{more than one network}
  \STATE paste the network name at the front of all the objective function
  effects
  \ENDIF
  \STATE Create the effects data frame by calling \sfn{createObjEffectList} and
  \sfn{createRateEffectList}.  This creates e.g.\ the evaluation and endowment
  effect copies.
  \STATE Alter the text for endowment effects to start "Lost ties:"
  \STATE Calculate the starting values for the default effects
  (see section \ref{sec:bipartitestart})
  \STATE Select the default rate effects by setting \nnm{include} to TRUE for
  the \sfn{basic rate} effects.
  \STATE Add the starting value for the rate to the \nnm{initialValue} column
  of the basic rate effects
  \IF{no period is uponly or downonly}
  \STATE  Set \nnm{include} to TRUE for the degree (density) evaluation effect
  \STATE Add the starting values for the degree(density) evaluation effect
  calculated by \sfn{getBipartiteStartingVals} to the \nnm{initialValue}.
  \ELSE
  \STATE Remove both degree (density) effects from the data frame.
  \ENDIF
\end{algorithmic}
\subsection{covarOneModeEff}
\begin{algorithmic}
\STATE Use \sfn{covarSymmetricObjective} or \sfn{covarNonSymmetricObjective}
effect group to create the objective function effects
\STATE Use \sfn{covarSymmetricRate} or \sfn{covarNonSymmetricRate} to create the
rate effects
\IF {not poszvar}
\STATE Reduce the new objective function effects to just ``altX'' and ``altSqX'' (symmetric)
or ``egoX'' (non symmetric)
\ENDIF
\IF{ not morethan2}
\STATE Remove the ``altSqX'' effect.
\ENDIF
\end{algorithmic}
\subsection{covarBipartiteEff}
\begin{algorithmic}
\IF {first node set matches}
\STATE Use \sfn{covarBipartiteRate} effect group to create the
rate effects
\ENDIF
\STATE Use \sfn{covarBipartiteObjective} effect group to create the objective
function effects
\IF{first node set matches}
\STATE reduce the objective function effects to ``egoX'', ``altDist2'', and ``totDist2''
\ELSIF {poszvar}
\STATE reduce the new objective function effects to the rows ``altX'' and
``altSqX''
\IF{ not morethan2}
\STATE remove the ``altSqX'' effect
\ENDIF
\ELSE
\STATE no objective function effects
\ENDIF
\end{algorithmic}
\subsection{covBehEff}
\begin{algorithmic}
\STATE Use \sfn{covarBehaviorObjective} effect group to create a potential set
of objective function effects
\IF {covariate and behavior variable are different}
\STATE  Create objective function effects as the first row of potential set
\ENDIF
\FORALL{oneMode dependent variables with the same node set}
\STATE Add an objective function effect from the second row of the potential
set.
\ENDFOR
\IF {any objective function effects}
\STATE Set \nnm{shortName} to \sfn{effFrom}
\ENDIF
\STATE Use \sfn{covarBehaviorRate} effect group to create the
rate effects
\STATE Use \sfn{covarABehaviorBipartiteObjective} effect group to create
objective function effects for bipartite dependent networks combined
with covariates on the first node set
 \end{algorithmic}
\subsection{covBBehEff}
\begin{algorithmic}
\STATE Use \sfn{covarBBehaviorBipartiteObjective} effect group to create
objective function effects for bipartite dependent networks combined
with covariates on the second node set
 \end{algorithmic}
\subsection{covarNetNetEff}
\begin{algorithmic}
\IF {\nnm{poszvar} }
\STATE  Use \sfn{covarNetNetObjective} effect group to create additional
objective function effects if the second network is one-mode;
\STATE  Use \sfn{covarABNetNetObjective} effect group to create additional
objective function effects if the second network is one-mode or two-mode;
\STATE  Use \sfn{covarANetNetObjective} effect group to create additional
objective function effects if the second network is one-mode or
\{two-mode while the covariate is defined for the first mode\};
\STATE  Use \sfn{covarBNetNetObjective} effect group to create additional
objective function effects if the second network is one-mode or
\{two-mode while the covariate is defined for the second mode\}.
\ENDIF
\end{algorithmic}
\subsection{CreateRateEffectList}
Add the \nnm{name} column by duplicating the dependent variable name, and
\nnm{effectFn} and \sfn{statisticFn} as empty lists.
\subsection{CreateObjectEffectList}
\begin{algorithmic}
\STATE
Add the \nnm{name} column by duplicating the dependent variable name, and
\nnm{effectFn} and \sfn{statisticFn} as empty lists.
\STATE Add an endowment effect row if required for each objective function
effect.
\end{algorithmic}
 \subsection{SienaGroupObject}
 \begin{algorithmic}
 \STATE First create the effects for the first data object, but inserting the
 correct number of basic rate effects for the whole group.
\FORALL{other data objects}
 \FORALL{dependent variables}
 \STATE Create the starting values for this dependent variable
 \STATE Insert the rate starting values in the \nnm{initialValue} field of the
 correct effects
 \STATE Combine the starting values to create an overall one for the objective
 function effects:
 \IF{behavior}
 \STATE Add new $d_i, i=1,\ldots,n-1$ to make one long vector of
 difference vectors between $n$ observations
 \IF {rounded range of variable (max-min) is 2 (what happens with range 1)}
 \STATE Add to $n_{min+} = \sum_i n_{i,min+}$ and the others
 \STATE Tendency is
 $$  \log(\frac{(n_{min+} + 2 ) * (n_{max+} + n_{max0} + 4)))}
 {(n_{max-} + 2 ) * (n_{min+} + n_{min0} + 4)))}
 $$
 \IF {abs(tendency) $>2$}
 \STATE trim to $\pm 2$
 \ENDIF
 \ELSE[range less than 2 or greater than 2]
 \STATE Let $\bar{d} = \mathrm{mean}(d_i)$, ignoring missing values
 \STATE Let $\sigma_d^2 = \mathrm{var} (d_i)$, ignoring missing values
 \IF{$\bar{d} < 0.9 * \sigma_d^2$}
 \STATE tendency = $0.5 * \log((\bar{d} + \sigma_d^2)/(\bar{d} - \sigma_d^2))$
 \ELSE
 \STATE tendency = $\bar{d} /(\sigma_d^2+ 1)$
 \ENDIF
 \IF {abs(tendency) greater than 3)}
 \STATE Trim to $\pm 3$
 \ENDIF
 \ENDIF
\ELSIF{onemode}
\STATE
\ELSE[bipartite]
\STATE
\ENDIF
 \ENDFOR
 \ENDFOR
\end{algorithmic}
\subsection{Behavior Starting Values}
\label{sec:behaviorstart}
\begin{algorithmic}
  \STATE Calculate $d_i, i=1,\ldots,n-1$
  difference vectors between $n$ observations
  \IF {rounded range of variable (max-min) is 2 (what happens with range 1)}
  \FORALL{intervals $i$}
  \STATE Let $n_{i,min+} =$ number who start at minimum and go up
  \STATE Let $n_{i,min0} =$ number who start at minimum and stay there
  \STATE Let $n_{i,max-} =$ number who start at maximum and go down
  \STATE Let $n_{i,max0} =$ number who start at maximum and stay there
  \STATE Calculate $$ v = \frac{n_{i, min+} +1}{n_{i, min+} + n_{i,min0} + 2} +
  \frac{n_{max-} +1}{n_{i,max0} + n_{i,max-} + 2}$$
  \IF {$v>0.9$}
  \STATE$v=0.5$
  \ENDIF
  \ENDFOR
  \STATE Starting rate is $-\log(1-v)$
  \STATE Let $n_{min+} = \sum_i n_{i,min+}$ total number who start at minimum
  and go up
  \STATE Let $n_{min0} = \sum_i n_{i,min0}$ total number who start at minimum and
  stay there
  \STATE Let $n_{max-} = \sum_i n_{i,max-}$ total number who start at maximum
  and go down
  \STATE Let $n_{max0} = \sum_i n_{i,max0}$ total number who start at maximum
  and stay there
  \STATE Tendency is
  $$  \log(\frac{(n_{min+} + 2 ) * (n_{max+} + n_{max0} + 4)))}
  {(n_{max-} + 2 ) * (n_{min+} + n_{min0} + 4)))}
  $$
  \IF {abs(tendency) $>2$}
  \STATE trim to $\pm 2$
  \ENDIF
  \ELSE[range less than 2 or greater than 2]
  \FORALL{intervals i}
  \STATE starting rate is $\max(\mathrm{var}(d_i),
  0.1 * \sum_i \mathrm{abs}(d_i)/\mathrm{nactors} $
  \ENDFOR
  \STATE Let $\bar{d} = \mathrm{mean}(d_i)$, ignoring missing values
  \STATE Let $\sigma_d^2 = \mathrm{var} (d_i)$, ignoring missing values
  \IF{$\bar{d} < 0.9 * \sigma_d^2$}
  \STATE tendency = $0.5 * \log((\bar{d} + \sigma_d^2)/(\bar{d} - \sigma_d^2))$
  \ELSE
  \STATE tendency = $\bar{d} /(\sigma_d^2+ 1)$
  \ENDIF
  \IF {abs(tendency) greater than 3)}
  \STATE Trim to $\pm 3$
  \ENDIF
  \ENDIF
\end{algorithmic}
\subsection{One mode network Starting Values}
\label{sec:onemodestart}
\begin{algorithmic}
\STATE Temporarily subtract 10 from structural values
\STATE Let $dif_i$ be the number of differences between start and end of
interval $i$, ignoring missings
\STATE Let $n_{ijk}, j, k = 0,1$ be counts of cells with value $j$ at start
and $k$ at end of interval $i$, ignoring missings
\STATE Let $n_i$ be the number of cells which are not missing at both start
and end of interval $i$.
\STATE Let $d_i$ be the sum of absolute differences by period, ignoring missings
(already calculated and stored in the attribute \nnm{distance}).
\STATE Let $\lambda_i$ be the starting value of basic rate parameter for
interval $i$
\IF{symmetric}
\STATE
$$\lambda_i = \mathrm{nactors} * (0.2 + d_i)/(n_i \%/\% 2 + 1)$$
\ELSE
\STATE
$$\lambda_i = \mathrm{nactors} * (0.2 + 2 * d_i)/(n_i + 1)$$
\ENDIF
\STATE Trim $\lambda_i$ to be between 0.1 and 100.
\IF{symmetric}
\STATE Divide $n_{ijk}$ by 2
\ENDIF
\STATE starting value for degree parameter:
\begin{align*}
\mathrm{Define } \qquad  p_{i01} &=& \begin{cases}
n_{i01}/ (n_{i01} + n_{i00}) &n_{i01} + n_{i00} >= 1  \\
0.5& \mathrm{otherwise}
\end{cases}\\
p_{i10} &=& \begin{cases}
n_{i10}/ (n_{i10} + n_{i11}) &n_{i10} + n_{i11} >= 1  \\
0.5& \mathrm{otherwise}
\end{cases}\\
  p_{i00} &=& \begin{cases}
n_{i00}/ (n_{i01} + n_{i00}) &n_{i01} + n_{i00} >= 1  \\
0.5& \mathrm{otherwise}
\end{cases}\\
p_{i11} &=& \begin{cases}
n_{i11}/ (n_{i10} + n_{i11}) &n_{i10} + n_{i11} >= 1  \\
0.5& \mathrm{otherwise}
\end{cases}
\end{align*}
\STATE Trim $p_{ijk}$ to lie between 0.02 and 0.98
\STATE Calculate $p_i$
\begin{align*}
p_i &=& \begin{cases}
4 / (p_{i00} / n_{i01} + p_{i11} / n_{i10})& n_{i10} * n_{i01} >= 1  \\
1e-6& \mathrm{otherwise}
\end{cases}
\end{align*}
\STATE Starting value for degree parameter is\\
$$ \sum_i 0.5 * \log(p_{i01} / p_{i10}) * p_i / \sum_i p_i $$
\end{algorithmic}

\subsection{Bipartite Starting Values}
\label{sec:bipartitestart}
\begin{algorithmic}
\STATE Temporarily subtract 10 from structural values
\STATE Let $dif_i$ be the number of differences between start and end of
interval $i$, ignoring missings
\STATE Let $n_{ijk}, j, k = 0,1$ be counts of cells with value $j$ at start
and $k$ at end of interval $i$, ignoring missings
\STATE Let $n_i$ be the number of cells which are not missing at both start
and end of interval $i$. (Diagonal included here.)
\STATE Let $d_i$ be the sum of absolute differences by period, ignoring missings
(already calculated and stored in the attribute \nnm{distance}).
\STATE Let $\lambda_i$ be the starting value of basic rate parameter for
interval $i$
\STATE
$$\lambda_i = \mathrm{nsenders} * (0.2 + 2 * d_i)/(n_i + 1)$$
\STATE Trim $\lambda_i$ to be between 0.1 and 100.
\STATE starting value for degree parameter:
\begin{align*}
\mathrm{Define } \qquad  p_{i01} &=& \begin{cases}
n_{i01}/ (n_{i01} + n_{i00}) &n_{i01} + n_{i00} >= 1  \\
0.5& \mathrm{otherwise}
\end{cases}\\
p_{i10} &=& \begin{cases}
n_{i10}/ (n_{i10} + n_{i11}) &n_{i10} + n_{i11} >= 1  \\
0.5& \mathrm{otherwise}
\end{cases}\\
  p_{i00} &=& \begin{cases}
n_{i00}/ (n_{i01} + n_{i00}) &n_{i01} + n_{i00} >= 1  \\
0.5& \mathrm{otherwise}
\end{cases}\\
p_{i11} &=& \begin{cases}
n_{i11}/ (n_{i10} + n_{i11}) &n_{i10} + n_{i11} >= 1  \\
0.5& \mathrm{otherwise}
\end{cases}
\end{align*}
\STATE Trim $p_{ijk}$ to lie between 0.02 and 0.98
\STATE Calculate $p_i$
\begin{align*}
p_i &=& \begin{cases}
4 / (p_{i00} / n_{i01} + p_{i11} / n_{i10})& n_{i10} * n_{i01} >= 1  \\
1e-6& \mathrm{otherwise}
\end{cases}
\end{align*}
\STATE Starting value for degree parameter is\\
$$ \sum_i 0.5 * \log(p_{i01} / p_{i10}) * p_i / \sum_i p_i $$
\end{algorithmic}

\section{sienaGroupCreate}

This function combines a list of siena data objects for common processing.
It is also used in \sfn{initializeFRAN} to convert a single data object to a group one so that all later
processing can have a common argument.

Some validation is performed to check that all the data objects match in terms
of the dependent variables (name, type, nodesets) and covariates (names and
nodesets).

If there is more than one data object, the constant covariates and constant
dyadic covariates must be changed into changing ones. New covariates are created
and the old ones removed from the lists. The attributes are copied over rather
than recalculated (although they are mostly changed later).

Overall values are calculated for the balance mean, and network ranges. For
behavior variables and covariates overall ranges, means and similarity mean
values are calculated.

The following overall values are copied down to the individual objects:
\begin{description}

\item[dependent variables] symmetric, missing, structural, poszvar, range,
moreThan2 (some only if behavior)
\item[changing covariates] range,  poszvar, moreThan2
\item[changing dyadic covariates] range, range2

\end{description}
\subsection{Group attributes}

The group object has various attributes, copied from or combinations of the
attributes on the individual objects.

\begin{description}
\item[netnames] the names of the dependent variables (these must be the same in
  each data object)
\item[symmetric] logical vector indicating whether the corresponding independent
  variable is symmetric (set to FALSE for behavior variables and bipartite
  networks.
\item[structural] logical vector indicating presence or absence of any
  structural values
\item[numberNonMissingNetwork] vector of count of non missing values for non
  behavior variables
\item[numberMissingNetwork] vector of count of missing values for non
  behavior variables
\item[numberNonMissingBehavior] vector of count of non missing values for
  behavior variables
\item[numberMissingBehavior] vector of count of missing values for
  behavior variables
\item[allUpOnly] logical vector indicating that the values for this
  independent variable never decrease over time
\item[allDownOnly] logical vector indicating that the values for this
  independent variable never increase over time
\item[anyUpOnly] logical vector indicating that for one or more of the intervals
  the values for this independent variable do not decrease
\item[anyDownOnly] logical vector indicating that for one or more of the
  intervals the values for this independent variable do not increase.
\item[allHigher] Table of logicals indicating whether \emph{higher} attribute is
  true for each pair of networks in every data object
\item[allDisjoint] Table of logicals indicating whether \emph{disjoint}
  attribute is true for each pair of networks in every data object
\item[allAtLeastOne] Table of logicals indicating whether \emph{atLeastOne}
  attribute is true for each pair of networks in every data object
\item[anyHigher] Table of logicals indicating whether \emph{higher} attribute is
  true for each pair of networks in any data object
\item[anyDisjoint] Table of logicals indicating whether \emph{disjoint}
  attribute is true for each pair of networks in any data object
\item[anyAtLeastOne] Table of logicals indicating whether \emph{atLeastOne}
  attribute is true for each pair of networks in any data object
\item[types] types of the independent variables
\item[observations] A single integer with the total number of periods to
  process.
\item[periodNos] A list of the period numbers (misses out the final one for each
  data object)
\item[groupPeriods] Vector of the number of total number of periods for each
  data object.
\item[netnodesets] A list containing the nodeset(s) for each dependent
  variable
\item[cCovars] Vector of names of constant covariates
\item[vCovars] Vector of names of changing covariates
\item[dycCovars] Vector of names of constant dyadic covariates
\item[dyvCovars] Vector of names of changing dyadic covariates
\item[ccnodesets] A vector containing the nodeset(s) for each constant covariate
\item[cvnodesets] A vector containing the nodeset(s) for each changing covariate
\item[dycnodesets] A list containing the nodeset(s) for each constant dyadic
  covariate
\item[dyvnodesets] A list containing the nodeset(s) for each changing dyadic
  covariate
\item[compositionChange] Logical vector indicating the presence of composition
  change data for any of the data objects
\item[exooptions] Named vector of composition change file options for the named
  node sets. Read from the first data object: assumed all the same
\item[names] Vector of names of the data objects
\item[class] ("sienaGroup", "siena")
\item[balmean] Vector of overall balance means, one for each dependent variable,
  NA for behavior variables and bipartites.
\item[structmean] Vector of overall structural means, one for each dependent
  variable, NA for behavior variables and bipartites.
\item[averageOutDegree] Vector of overall average outdegrees. NA for behavior
  variables.
\item[averageInDegree] Vector of overall average indegrees. NA for behavior
  variables and bipartites.
\item[bRange] Vector of overall ranges for behavior variables: entries
  corresponding to networks are NA
\item[behRange] Matrix with two rows, and column for each independent
  variable. Set to the overall min and max for behavior variables, NA for others
\item[bSim] Overall similarity mean for behavior variables, NA for networks
\item[bPoszvar] logical vector, NA for networks. For behavior variables TRUE if
  more than 1 distinct value in the overall values or any missing (always?)
\item[bMorethan2] logical vector. NA for networks. For behavior variables
  TRUE if more than 2 distinct values in overall variable, ignoring missings
\item[cCovarPoszvar] logical vector for constant covariates indicating presence
  if more than one distinct value or any missing (overall). NB only exist if
  there is only one data object.
\item[cCovarMoreThan2] logical vector for constant covariates indicating
  presence of more than 2 distinct values (missing is counted as a value)
\item[cCovarRange] vector of ranges for constant covariates
\item[cCovarRange2] matrix of min and max for constant covariates
\item[cCovarSim] vector of overall similarity means for constant covariates
\item[cCovarMean] vector of means for constant covariates
\item[vCovarPoszvar] logical vector for changing covariates indicating presence
  if more than one distinct value or any missing (overall).
\item[vCovarMoreThan2] logical vector for changing covariates indicating
  presence of more than 2 distinct values (missing is counted as a value)
\item[vCovarRange] vector of ranges for changing covariates
\item[vCovarSim] vector of overall similarity means for changing covariates
\item[vCovarMean] vector of means for changing covariates
\item[dycCovarRange] vector of ranges for constant dyadic covariates
\item[dycCovarRange2] matrix of min and max for constant dyadic covariates
\item[dycCovarMean] vector of means for constant dyadic covariates
\item[dyvCovarRange] vector of ranges for changing dyadic covariates
\item[dyvCovarRange2] matrix of min and max for constant dyadic covariates
\item[dyvCovarMean] vector of means for changing dyadic covariates
\item[anyMissing] logical vector indicating any missing values in the each
  dependent variable
\item[netRanges] Matrix with two rows and a column for each dependent
  variable. Overall min and max for networks, NA for behavior variables
\end{description}


\section{siena07}

This is a wrapper for the function \nmm{robmon} which performs the processing
that used to be in \nm{polrup} in \Sn.  An optional tck/tk gui is provided,
or progress messages are provided on the console. The choice between these is
made by using \sfn{batch=FALSE} or \sfn{batch=TRUE} respectively.

Details of input and output are on the \R help page. Required input is an object
containing control information for the Robbins-Monro algorithm, and any extra
parameters required by the \nm{FRAN} to be used. As the distinction between the
two parts is not complete, flags \nm{maxlike} and \nm{cconditional} are on the
input object, although not logically relevant to the algorithm.

There is user output written to a file (.txt), together with optional additional
output to the console (suppressed unless \sfn{verbose=TRUE}), which can be
redirected using the \verb|sink()| command.

\nnm{robmon} attempts to duplicate the output of
the \SI procedure \nm{polrup}. It uses
a special function, \nnm{Report} for all output.
This function knows about four files:
\nm{outf}, \nm{lf}, \nm{cf}, \nm{bof}, and can also write to the console.
Currently, all files except \nm{outf} are null, with any other output suppressed
or written to the console. Only \nm{Report} would need altering to alter this
behaviour.  No file connections need to be passed around as parameters.

The object returned from \nnm{siena07} is an object containing everything of
interest from the run, including the estimates of the parameters and the
covariance matrix.  Details of the more useful parts are in the \R
documentation.

\section{User Interrupts}
These are set in callbacks from the \nnm{siena07} tcl/tk gui. When they are
read, a parallel set of flags is used to store the states, so that interrupts
can be processed reliably. All is done using functions, to avoid global
variables, or passing variables around. There are 3 interrupts:
\begin{description}
\item[UserInterrupt] Stop everything, but return the values so far, with (I
hope) some flag to indicate we did not finish.
\item[UserRestart] Go back to the beginning of phase 1 with the current
  parameters
\item[EarlyEndPhase2] Stop the estimation routine and proceed to phase 3 using
  the current parameters.
\end{description}
All six functions, if called with an argument, store the argument as the current
value and if called with no argument, return the value. All values are
booleans. The functions are not exported from the namespace, to avoid burdening
the user with the details of their existence.

\section{Robbins Monro Algorithm---\nnm{robmon}}

The routine \nmm{robmon} contains the Robbins Monro
(stochastic approximation) algorithm. It is the replacement for the
\SI procedure \nm{polrup}. It is not
designed to be called directly by the user, so there will not be an \R
help page for it.

The outline of the algorithm is given in the text \emph{Siena algorithms},
and to understand the description of the code given here
it may be helpful to have read that outline.


\section{ \nnm{robmon}}
\subsection{Input (from \nm{siena07})}
\begin{description}
\item[\nnm{z}] Model fitting object.
\item[\nnm{x}] Input model object, as described in the help page for
 \nnm{sienaModelCreate}.
\item[...] Extra parameters for \nm{FRAN} (including the data!).
\end{description}

It may seem surprising that the data is simply a parameter passed unchanged to
\nm{FRAN}. But this is the point of the separation of the simulation and
estimation routines: \nnm{robmon} could be used to solve the moment equation for
any data: it does not matter whether the data is a network or something else
entirely, as long as a matching \nm{FRAN} is used.

\subsection{Output}
\begin{description}
\item[\nnm{z}]. More or less everything used in the processing.
Details in the help page for \nnm{siena07}.
\end{description}

\subsection{Details}

\subsubsection{Initialize}
 Copy from \nnm{x} everything that we may change
  during the run. Initialize number of iterations, restarted flag,
  force finite difference flag, etc.

\subsubsection{Initial call to \nm{FRAN}}
This call is used to set up the parameters for conditional estimation, and to
set up data in a call to C. The values of the statistics in the observed data
(targets) are returned, along with the addresses of the data objects in C. The
processed data and selected effects is written to a hidden data object within
the function \nnm{FRANstore} from where it can be accessed on later calls, or
passed to other processes.

\subsubsection{Initialise cluster of processes, if required}
This needs to access (but not understand!) the data object created in the call
to \nm{FRAN}. It sets up the processes and random number streams and passes the
data object across. It then does a special call to \nm{FRAN} to create the data
objects in C++ for each process. Later calls to the processes only need minimal
communication, done using cut-down versions of \nnm{x} and \nnm{z}.

\subsubsection{Calculate epsilon}
Used only for MoM estimation in the finite differences option.
Currently 0.1, except for parameters which
must be positive, where it is 0.1 times the parameter starting values.
\ts{Here I would prefer min\{0.1, 0.1 $\times$ starting value\}.
\\Ruth: Easily changed, but we are using typically much bigger values than this:
will it work? Are you sure you don't want max?}

To be
improved, to use prior information on standard error of parameter if available.
  \newpage
\subsubsection{Main loop}
Note \R has no GOTO statement. I use the term \emph{break} to indicate
exit from the current loop only. Interrupts are checked after every iteration
except the first few of phase 3. This documentation does not include all the
details, or it would duplicate the code.

\vspace{1em}

{\small
\begin{algorithmic}
  \REPEAT
  \REPEAT [this one is just to jump out of, only executed once]
  \IF{ all parameters now fixed (2 opportunities to do this in previous loop!)}
      \STATE set a flag to just do Phase 3.\\
  \ENDIF
  \STATE initialize interrupts\\
  announce phase 0 (set up progress bar, calculate iteration min
  and max for phase 2 subphases.)\\
  reset fixed flags\\
  \IF {not just-phase-3 flag set}
      \IF {need to do phase 1}
          \STATE initialize phase 1\\
          run phase 1 iterations 1 to 10\\
          \IF{using finite diffs}
               \STATE check number of changes and
                  change epsilon if necessary
          \ENDIF
          \IF{ user stop or user restart or error or (using finite difference
            and need to repeat with new epsilon)}
               \STATE break
          \ENDIF
          \IF{using finite diffs}
          \STATE fix PARAMETERS with 0 or 1 changes, if any exist\\
          \ENDIF
          \STATE run rest of phase 1 iterations\\
          \IF {user stop or user restart or error}
              \STATE break
          \ENDIF
           \STATE calculate derivative matrix
          \IF {necessary}
              \STATE change the length of phase 1 or force the use of
              finite differences
          \ENDIF
          \IF{ user stop or error or user restart}
              \STATE break
          \ENDIF
          \IF{ necessary}
              \STATE fix some parameters.
          \ENDIF
      \ENDIF
       \STATE Initialise phase 2\\
       run phase 2 subphase 1\\
       \IF{ error or user stop or user restart}
           \STATE break
       \ENDIF
       \STATE run phase 2 subphase 2
       \IF {error or user stop or user restart or we have restarted because of
         epsilon change in phase 1 and not restarted from here
         before!}
       \STATE break
       \ENDIF
       \STATE    run rest of phase 2 subphases
       \IF {error or user stop or user restart}
       \STATE break
       \ENDIF
\ENDIF
      \STATE run phase3
      \IF {not user restart}
      \STATE break
      \ENDIF
    \UNTIL{for ever}
    \IF{do not need to restart because of epsilon or user restart}
\STATE break
\ENDIF
\UNTIL{for ever}
\end{algorithmic}}

\subsubsection{Final processing}
\begin{itemize}
\item
Do a final call to \nm{FRAN}. In conditional estimation, the rescaling
of basic rate parameters will be done here.
\item reset the covariance matrix to 33, 999 as in phase3.

\end{itemize}
\section{Phase 1}
\subsection{Input}
As for \nnm{robmon}.
\subsection{Output}
As for \nnm{robmon}.
\subsection{Details}
\subsubsection{Initialise}
\begin{itemize}
\item
Reset \nnm{SomeFixed} flag (have we fixed any parameters in this
run).
\item Announce phase 1
\item Create arrays to store simulated statistics,
scores and contributions to the derivative matrix from either the
finite difference or maximum likelihood routines. These arrays are
currently redefined in Phase 3, so lost.
\end{itemize}
\subsubsection{Timing}
Timings are calculated between the start of the 2nd iteration and the start of
the 6th iteration.  This is just for determining the frequency of writing
information to the gui.  Write frequency is set to a prettified version of
20/time for 5 iterations, or 5 if elapsed time is very small. For batch mode
this is multiplied by 10, which seems unnecessary in phase 1. If using multiple
clusters the total number of iterations are adjusted to be a multiple of the
number of processes, and the 6th is replaced by the first one greater than or
equal to 6 in the iteration sequence advancing in steps of the number of
processes. I think this is wrong: it should be from 2 to 6 steps... but not very
important!  (I have made some adjustments for those users who have more than 9
processors: we do at least 10 simulations in the first part, and ignore timing
if we do them in too few steps. If there are enough steps in the second part the
timing is done there.)
\subsubsection{An iteration}
\begin{itemize}
\item call \nm{FRAN}. If not OK, return
\item store simulated statistics, scores if present, derivatives if
  present, simulated networks (part!) if present.
\item if required, call finite differences routine and store result
\item check for user interrupts and return if requested
\item Report progress via progress bar or to console
\end{itemize}
\subsubsection{Check epsilons (only for finite differences option)}
\begin{itemize}
\item
If derivatives are being calculated by finite differences, check after
10 iterations that enough different values of the statistics have
occurred. Ideally 5 or more. (The check may not be done immediately after 10
iterations if we have 4 processors, say, but I do only look at 10 of the
results.)
\item
If there are less than 3, \nm{epsilon} is
multiplied by 3 for parameters which must be positive, and by 10
otherwise.
\item
If there are 3 or 4, the multipliers are 2 and $\sqrt{10}$.
\item
If any new values are less than 0.1 times the scale factor, or more
than 100 times it, replace by the bound.
\item
If any are less than 5, we will restart with the new epsilons, unless
we have already done so 4 times.
\item
If we have already restarted 4 times then we continue, after fixing
 the last parameter with only 0 or 1 changes, if there are any such.
\end{itemize}
\subsubsection{End of phase processing}
\begin{itemize}
\item Calculate derivative estimate.
\begin{itemize}
\item
For finite differences or maximum likelihood: the mean of the arrays
returned at each iteration.
\item For the score based method, we need a little notation:\\
  Let $\mathbf{f}_i$ be the simulated deviations from the targets in
  iteration $i$, and $\mathbf{s}_i$ the score function in iteration
  $i$, $N$ be the number of iterations. Then the estimated derivative
  matrix $d_{ij}$ is given by:
$$
D = \sum_i \mathrm{outer}(\mathbf{f}_i,\mathbf{s}_i)/ N -
\mathrm{outer}(\bar{\mathbf{f}}, \bar{\mathbf{s}})/N^2
$$
\item For the score based method, if any of the diagonal values is
  non-positive, we do not continue. First we double the number of
  iterations in phase 1 and start again. Once this number exceeds 200,
  we stop increasing it and force the use of Finite differences. Then
  we don't come through this check!
\item
For either method, if still processing, set the rows and columns related to fixed
parameters to 0's with 1 on the diagonal.
\item
If any diagonal values for non-fixed parameters are not positive, make
them fixed, and set \nnm{newfixed} flags to record which ones have
been fixed.
\item
Calculate the standard deviations of the deviations.
\end{itemize}
\item Invert the derivative matrix:
\begin{itemize}
\item
  Set the rows and columns related to fixed parameters to 0's with 1 on the diagonal.
\item
  Replace any diagonal values less than 1e-8 by 1e-3.
\item
 Do the inversion
\item
If it fails, add 1 to the diagonal and try again
\end{itemize}

\item Quasi-Newton step
\begin{itemize}
\item If inversion of matrix was successful, set \nm{fchange} to 0.5
  times the gain parameter times the matrix product of inverse with
  the mean deviations from targets, otherwise to zero.
\item
Zero the change for any fixed parameters.
\item
Check the jump is not too large: if the maximum absolute value of the
change divided by the corresponding input gain parameter is greater
than 10, divide the changes by this value and multiply them by
10. This caps the maximum ratio to the scale at 10.
\item Check that positive parameters will stay positive. If not, replace the
  change for that parameter by half the current value of the parameter.
\item
If the requested number of subphases in phase 2 is greater than 0, make
the change by subtracting the changes
from the current value of the parameters.
\end{itemize}
\end{itemize}
\section{Phase2}
\subsection{Input}
as for \nnm{robmon}
\subsection{Output}
as for \nnm{robmon}
\subsection{Details}
\subsubsection{Initialize Phase}
Turn off calculation of derivatives. Multiply the gain parameter by 2.
\subsubsection{Process subphase}
\begin{itemize}
\item Initialize
  \begin{itemize}
  \item
    Announce subphase
  \item Extract max and min number of iterations from values stored
    from start. (They are calculated at the start to find the length
    of the progress bar, and it seemed better not to repeat it here!).
  \item Divide the gain parameter by 2
  \item Initialize the sum (which will become the mean) of thetas with
    the current value.
  \item
    Create arrays to store products of successive thetas.
  \end{itemize}
\item
  Perform iterations
  \begin{itemize}
  \item Timing is calculated over the ten iterations 2 to 11. Write frequency is
    set to 20/(time for these 10 iterations) or 20 if time is too small.  It is
    then prettified. Reporting using progress bar or console is done for each of
    the first 10 iterations and then at write frequency.
\item
  Call \nm{FRAN} and store the deviations returned.
  \item The update:
    \begin{itemize}
    \item if only the diagonal of the derivative matrix is to be used
      in the update step (flag \nnm{x\$diag}), calculate \nm{maxrat},
      the maximum ratio of the absolute deviations to their standard
      deviation (as estimated in Phase 1). If this is greater than the
      parameter \nm{maxmaxrat}, set \nm{maxrat} to
      \nm{maxmaxrat}/\nm{maxrat}, and record that the values were
      truncated. \label{sec:maxrat}
    \item
      if \nmm{x\$diag},
      update is current gain * current deviations * \nm{maxrat} /
      diagonal of derivative
    \item otherwise, update is current gain * matrix product of
      current deviations with inverse of derivative matrix.
    \item For parameters which must be positive and would not remain
      so, replace change by 0.5 times current value of parameter,
    \item Zero the change for any fixed parameters.
    \item update theta by subtracting the change.
    \item Add new value of theta to sum of thetas.
    \end{itemize}
  \item
    After each pair of iterations, add the product of the two deviations
  and the square of the most recent to accumulators, and calculate the
  ratio of the former to the latter, \nm{ac}.
  \item
    Check for user interrupts
  \item
    Stop the subphase when either the minimum number of iterations has
    been reached, and the maximum of \nm{ac} is less then 1e-10 or the
    maximum number of iterations has been reached;\\
    or (the next condition is a rare occurrence but
    helps a lot when it occurs) at least 50 iterations have
    been done and the minimum \nm{ac} is $< -0.8$ and we have not done the
    maximum number of subphase repeats already (or user interrupt
    or error)
\item Repeat the subphase if we stopped because of the minimum value of \nm{ac},
  unless we have already done the maximum number (set to 4).
  \end{itemize}
\item
End of subphase. Replace parameters by the average in the subphase. Report details.
\end{itemize}
\section{Phase3}
\subsection{Input}
as for \nnm{robmon}
\subsection{Output}
as for \nnm{robmon}
\subsection{Details}
\subsubsection{Initialize}
\begin{itemize}
\item
Initialize arrays for deviations, scores, derivatives.
\item Divide write frequency by the number of parameters for finite differences,
  2 for score derivatives. Leave unchanged for maximum likelihood. Then
  re-prettify.
\item
 Announce Phase
\end{itemize}
\subsubsection{Iterations}
\begin{itemize}
\item Update progress via console or progress bar for each of first 5
  iterations, then at the 10th and then at the write frequency.
\item Call \nm{FRAN} and store the deviations, scores, derivative
  contributions, and simulation values returned.
\item If using finite differences (we revert to the method requested by the user
  here, even if we altered it in phase 1 because the score method did not work),
  call FiniteDifferences routine and store the resulting differences.
\item After the 10th iteration check for user interrupts each time.
\end{itemize}

\subsubsection{End of Phase 3}
\begin{itemize}
\item
Calculate derivative matrix \nm{dfra}. Formulae as in Phase 1.
\item
Create a flag \nm{diver} for each parameter, set to true if all the
parameter values are fixed and the absolute value of the corresponding
diagonal of the derivative matrix is less than 1e-6. \label{sec:diver}.
\item
Calculate the covariance matrix of the simulated deviations, and
the autocorrelations between them. (Done in the CalculateDerivative3
routine).
\item
Report\ldots
\item Calculate $t$-values for deviations from targets, as mean
  deviation divided by sqrt of diagonal entry of covariance
  matrix. Use 0 for small values of deviations and 999 for small
  values of variance.
\item
Report $t$-values, and comments on their values.
\item For maximum likelihood, report autocorrelations
\item Calculate \nm{cov}, the covariance matrix of the estimates:
\begin{itemize}
\item For maximum likelihood, inverse of \nm[dfra] - variance matrix of
  deviations (here=scores!). (0 the rows and columns coresponding to
  fixed parameters, and put 1 on the diagonal, first.)
\item For others, matrix product of \nm{dinv}, the adjusted covariance
  matrix of the deviations and the transpose of \nm{dinv}.
\item Try to invert \nm{cov}. Report\ldots
\item Update the flag \nm{diver} to mean: Not all parameters are fixed
  and this one is fixed or \nm{diver} was true before (section
  \ref{sec:diver}) or the corresponding entry in the diagonal of the
  covariance matrix of the estimates is less than 1e-9. \label{sec:diver2}
\item
  set entries in \nm{cov} to 33* sqrt(diagonal) off diagonal and 999
  on diagonal for parameters with \nm{diver} true.
\end{itemize}
\item Do ScoreTests if any have been requested.
\end{itemize}
\section{FiniteDifferences}
\subsection{Input}
as for \nnm{robmon} plus \nnm{fra}, the simulated value of the
targets.
\subsection{Output}
as for \nnm{robmon}
\subsection{Details}
For each parameter,
\begin{itemize}
\item Call \nm{FRAN} with epsilon added to this parameter only
\item Calculate the deviations from the simulated statistics with no epsilon.
\item In first 10 iterations of phase 1: Record if difference is
  greater than 1e-06
\item Store and return these deviations divided by epsilon.
\end{itemize}.
\section{Score Tests}
\subsection{Control}
\begin{itemize}
\item Do general test: call \nm{EvaluateTestStatistic} with the complete
  arrays \nm{dfra}, \nm{msf}, the covariance matrix of the deviations,
  and \nm{fra} the mean deviations.
\item The values returned are a chi-squared test and, if the degrees of
  freedom are 1, a one-sided test.
\item If only one test was requested, the two values returned
  correspond to the results required.
\item If more than one test was requested, call \nm{EvaluateTestStatistic}
  with data from which all but one of the parameters for which tests
  are required have been removed. Repeat for each parameter for which
  tests were requested.
\end{itemize}
\subsection{EvaluateTestStatistic}
\begin{itemize}
\item
Partition \nm{dfra} into four: \R code is clear enough, I hope
(\nm{drop=FALSE} just retains the matrix class even if one of the
dimensions is 1)
\begin{verbatim}
d11 <- dfra[!test,!test,drop=FALSE]
d22 <- dfra[test,test,drop=FALSE]
d21 <- dfra[test,!test,drop=FALSE]
d12 <- t(d21)
\end{verbatim}
\item
Similarly create $\Sigma 11$,  $\Sigma 22$, $\Sigma 12$ and $\Sigma 21$
 from \nm{msf}, and \texttt{z1} and \texttt{z2} from
\nm{fra}. Then
\begin{align*}
\intertext{For maximum likelihood}
\mathtt{ov} &= -\mathtt{z2}\\
 \mathtt{vav} &= \left(\mathtt{d22} - \mathtt{d21} \
   \mathtt{d11}^{-1}\ \mathtt{d12}\right) ^{-1}
\intertext{Otherwise}
\mathtt{ov} &= \mathtt{z2} - \mathtt{d21}  \ \mathtt{d11}^{-1} \ \mathtt{z1}\\
\mathtt{vav} &= \left(\Sigma22 -\mathtt{d21}  \  \mathtt{d11}^{-1} \Sigma12 -
(\Sigma21 - \mathtt{d21}  \  \mathtt{d11}^{-1} \Sigma11)\
\mathtt{d11}^{-T} \ \mathtt{d21}^T\right )^{-1}
\intertext{then}
\text{test statistic} &= \mathtt{ov}^T  \ \mathtt{vav}\  \mathtt{ov}\\
\intertext{ and the one-sided one, if appropriate,}
&\mathtt{ov} \sqrt{\mathtt{vav}}
\end{align*}
as \texttt{vav} is then a scalar.
\end{itemize}
\end{document}
