###############################################################################
##################    RscriptSienaScoreTest.R    ##############################
###############################################################################
#
# This is an R script to illustrate the use of score-type tests in RSiena.
# It requires RSiena version at least 1.6.7.
#
# R script written by Tom A.B. Snijders.
# Version May 29, 2026.

###############################################################################
###############################################################################

library(RSiena)


###############################################################################
###############################################################################

# The first example illustrates how to use the score-type test
# for explorative forward model selection.

# We use the s50 data set in RSiena

?s50

# A basic model is specified:

mynet <- as_dependent_rsiena(array(c(s501, s502), dim=c(50, 50, 2)))
mydata <- make_data_rsiena(mynet)
myeff <- make_specification(mydata)
myeff <- set_effect(myeff, list(outAct, inPop, inAct, gwespFF, reciAct))
myeff <- set_interaction(myeff, list(recip, gwespFF))
myeff

# and estimated

myalgo <- set_algorithm_saom(seed=44003)
(ans0 <- siena(mydata, effects=myeff, control_algo=myalgo))

# Convergence is good.
# Suppose we wish to investigate
# whether a degree assortativity effect should be added.
# Look in the manual for "degree assortativity" to see what this means
# and for "assortativity" to see what are the shortNames.
# Two effects are added to the model as fixed-and-tested effects:

myeff1 <- set_effect(myeff, inInAss, parameter=1, fix=TRUE, test=TRUE)
myeff1 <- set_effect(myeff1, outOutAss, parameter=1, fix=TRUE, test=TRUE)

# Since only some fixed effects were added to the model specification,
# we do not need to re-estimate but use prevAns=ans0
# and an algorithm with 0 subphases in Phase 2,
# which will skip estimation entirely.

algo0 <- set_algorithm_saom(seed=4937, nsub=0)
(ans1 <- siena(mydata, effects=myeff1, control_algo=algo0, prevAns=ans0))

# The print of ans1 shows that the score-type test is not significant.

# This is also shown by
test_parameter(ans1, tested=TRUE, method="score")
# and for the separate effects by
test_parameter(ans1, tested=8, method="score")
test_parameter(ans1, tested=9, method="score")

# The score-type test also can be used to produce one-step estimates
# (this is possible since RSiena version 1.6.7):

test_parameter(ans1, tested=TRUE, method="score", onestep=TRUE)

# The one-step estimates are quick-and-simple, but imprecise,
# approximations to what could be obtained as "real" estimates.

myeff2 <- set_effect(myeff1, inInAss, parameter=1)
myeff2 <- set_effect(myeff2, outOutAss, parameter=1)
(ans2 <- siena(mydata, effects=myeff2, control_algo=myalgo, prevAns=ans0))
# This has not converged, so we repeat:
(ans2 <- siena(mydata, effects=myeff2, control_algo=myalgo, prevAns=ans2))

# A comparison between the estimates of ans2 and the one-step estimates
# shows that the parameter estimates are quite similar.
# This is not always the case!
# The similarity is greater according as
# the significance of the effects is smaller and there is less collinearity
# between the estimated and the fixed parameters.
# The one-step estimates can be put into the effects object
# and the fixed effects can be freed in the following way:

(myeff3 <- update_theta(myeff1, ans1, onestep=TRUE, freed=8:9))
# Look at
?update_theta
# to see what is done here (also available since version 1.6.7)
# Let us inspect how well this satisfies the moment equations.
# For this purpose we use algo0 and myeff3, which skips estimation,
# and no prevAns:
(ans30 <- siena(mydata, effects=myeff3, control_algo=algo0))
# Not very good!
# The maximum convergence ratio is larger than 1.
# It becomes a little bit less if we only take half a step:
(myeff31 <- update_theta(myeff1, ans1, onestep=TRUE, freed=8:9, r=0.5))
(ans31 <- siena(mydata, effects=myeff31, control_algo=algo0))

# However, the one-step estimates can also be used as initial values
# for the estimation (with myalgo, not algo0):
(ans3 <- siena(mydata, effects=myeff3, control_algo=myalgo))
# Note the use of myalgo, which performs estimation,
# and note that we just use myeff3 and no prevAns.
# Thanks to the use of the one-step estimates,
# convergence takes place in one run of siena,
# whereas the use of prevAns=ans0 for ans2 required two runs.

# We now have two estimates of the same model,
# with maximum convergence ratios less than 0.25.
# Compare the estimates:
round( cbind(coef(ans2, shortenNames=FALSE), coef(ans3)), 3)
# and the standard errors:
round(cbind(ans2$se , ans3$se), 3)
# How do the differences between the estimates
# compare to the average standard errors?
round((coef(ans2) - coef(ans3))/((ans2$se + ans3$se)/2), 3)
# All are less than 0.2. This is small, but not very small.
# It would be better to use estimates with a lower maximum convergence ratio
# and obtained with a higher value of n3
# to obtain more precise standard errors.

###############################################################################
###############################################################################

# The second example illustrates how the score-type test is used
# to investigate an effect that does not perform very well
# when added to the model.

(myeff4 <- set_effect(myeff2, gwespFB))
myeff4 <- set_interaction(myeff4, list(recip, gwespFB))
(ans4 <- siena(mydata, effects=myeff4, control_algo=myalgo,
                                prevAns=ans2))
(ans4 <- siena(mydata, effects=myeff4, control_algo=myalgo,
                                prevAns=ans4))
# This does not converge; the estimates become larger and larger.
# Let us try to fix and test.
myeff41 <- set_effect(myeff4, gwespFB, fix=TRUE, test=TRUE)
myeff41 <- set_interaction(myeff41, list(recip, gwespFB), fix=TRUE, test=TRUE)
(ans41 <- siena(mydata, effects=myeff41, control_algo=algo0, prevAns=ans2))
test_parameter(ans41, tested=TRUE, method="score", onestep=TRUE)
test_parameter(ans41, tested=4, method="score")
test_parameter(ans41, tested=12, method="score")

# It appears that the effects are not significant.
# We can conclude that it is superfluous
# to add gwespFB and gwespFB*recip to this model.

################################################################################
################################################################################

# For the third example we use a school from the data set
# collected by Andrea Knecht for her dissertation
# Knecht, A., 2008. Friendship Selection and Friends' Influence.
# Dynamics of Networks and Actor Attributes in Early Adolescence.
# PhD dissertation, University of Utrecht.
# Various articles were published about this dissertation, starting with
# Andrea Knecht, Tom A.B. Snijders, Chris Baerveldt, Christian E.G. Steglich,
# and Werner Raub (2010. Friendship and Delinquency: Selection and Influence
# Processes in Early Adolescence, Social Development, 19, 494-514.
# http://dx.doi.org/10.1111/j.1467-9507.2009.00564.x.
# The school used here is "13e".
# The data were slightly transformed in the same way as in the book chapter
# "Stochastic Network Modeling as Generative Social Science",
# Chapter 5 (p. 73-99) in
# "Handbook of Rigorous Theoretical and Empirical Sociology" (2022),
# edited by Klarita Gerxhani, Nan Dirk de Graaf, and Werner Raub.
# https://doi.org/10.4337/9781789909432

# The data can be downloaded as
download.file("https://www.stats.ox.ac.uk/~snijders/siena/school13e.RData",
                                               destfile="school13e.RData")
unzip("schools13e.RData")
load("school13e.RData")
# What do we have now:
school13e

################################################################################
################################################################################

# The third example illustrates score-type tests for a parameter that is
# difficult to estimate, but significant, and may be regarded as infinite.
# We use school 13e. This has 25 actors,
# of whom 2 did not take part in the data collection;
# the others were 8 boys and 15 girls.

# A basic model is
effs <- make_specification(school13e)
effs <- set_effect(effs, list(gwespFF,inPopSqrt,outActSqrt))
effs <- set_effect(effs, reciAct, parameter=2)
effs <- set_effect(effs, sameX, covar1="gnd")
effs <- set_effect(effs, egoX, covar1="wave3")
(ans <- siena(school13e, effects=effs, control_algo=myalgo))
(ans <- siena(school13e, effects=effs, control_algo=myalgo, prevAns=ans))
# When the out-isolate effect is added,
effs2 <- set_effect(effs, outIso)
# estimation diverges with a very large parameter estimate
# (I used thetaBound=100 because an earlier run showed this is necessary):
(ans2 <- siena(school13e, effects=effs2, control_algo=myalgo,
                                        sprevAns=ans, thetaBound=100))
# This can be explored using a score-type test:
effs2a <- set_effect(effs, outIso, fix=TRUE, test=TRUE)
(ans2a <- siena(school13e, effects=effs2a, control_algo=myalgo, prevAns=ans))
test_parameter(ans2a, tested=TRUE, method="score", onestep=TRUE)
# The score-type test shows that the effect is significant.

# Let us take one step from the estimate in ans2a,
# but keep the fixed parameter fixed at its new value
# (this is done by "freed=0"):
(effs2b <- update_theta(effs2a, ans2a, onestep=TRUE, freed=0))
# This is not recommended if there are more than one fixed parameter
# with very large one-step estimates. But here we have only one.
# For estimating by siena, we now should NOT use prevAns,
# because that would employ
# the parameter values in prevAns rather than those in effs2b.
(ans2b <- siena(school13e, effects=effs2b, control_algo=myalgo))
test_parameter(ans2b, tested=TRUE, method="score", onestep=TRUE)
# Now the out-isolate parameter is not significantly different
# from its new value at  -5.2382.
# However, the one-step estimate for the out-isolate parameter
# now has decreased to  -10.049.
# what happens if we follow the recipe one more?
(effs2c <- update_theta(effs2b, ans2b, onestep=TRUE, freed=0))
(ans2c <- siena(school13e, effects=effs2c, control_algo=myalgo))
# Again the out-isolate parameter is not significantly different
# from its new value at  -10.049.
# The other parameters are quite similar for ans2b and ans2c.

# We now take a big step and fix the parameter at -50.
effs2d <- set_effect(effs, outIso, fix=TRUE, test=TRUE, initialValue=-50)
(ans2d <- siena(school13e, effects=effs2d, control_algo=myalgo))
# (I repeated estimation to lower the value of the maximum convergence ratio):
(ans2d <- siena(school13e, effects=effs2d, control_algo=myalgo, prevAns=ans2d))

# It can now be concluded that the out-isolate parameter is
# significantly different from 0 (with p=0.024)
# but not from  -5.2382,  -10.049, or -50.
# A strong or rather strong negative value for this parameter
# is in agreement with this data set, but a 0 value is not.
# To understand this we look at the data
# and when we do this, we see that no actors with outdegrees 1 or 2,
# who may have been be at risk for becoming out-isolates (outdegree 0),
# had outdegree 0 at the next wave. This means that
# it seems that actors avoided being out-isolates.
# There were no actors who had outdegree 0
# (except for structural zeros).

# The result of any of these models ans2b, ans2c, ans2d can be used,
# because the  maximum convergence ratio is small enough,
# and the score-type test for the negative value
# in any of these models is not significant.
# This means that the parameter estimate
# for the out-isolate effect is significantly lower than 0,
# with a p-value obtained from ans2a, which can be rounded to p=0.03;
# the estimate has a large negative value but how large is not determined.
# The standard error can be reported as NA.

# To compare  models ans2b, ans2c, and ans2d,
# it is better to obtain even smaller maximum convergence ratios.

# Use an algorithm that will help better if the estimate is rather close
# (see the manual, Section 7.4: "Use of the algorithm arguments":
algo1 <- set_algorithm_saom(seed=384, nsub=1, n2start=5000,
                            n3=5000, firstg=0.05)
(ans2b1 <- siena(school13e, effects=effs2b, control_algo=algo1, prevAns=ans2b))
(ans2c1 <- siena(school13e, effects=effs2c, control_algo=algo1, prevAns=ans2c))
(ans2d1 <- siena(school13e, effects=effs2d, control_algo=algo1, prevAns=ans2d))

# Let us try to get more precise estimates to make a comparison
# between the fixed values of the
# Use an algorithm that will help better if the estimate is rather close
# (see the manual, Section 7.4: "Use of the algorithm arguments":
algo1 <- set_algorithm_saom(seed=384, nsub=1, n2start=10000,
                            n3=10000, firstg=0.02)
(ans2b1 <- siena(school13e, effects=effs2b, control_algo=algo1, prevAns=ans2b))
(ans2c1 <- siena(school13e, effects=effs2c, control_algo=algo1, prevAns=ans2c))
(ans2d1 <- siena(school13e, effects=effs2d, control_algo=algo1, prevAns=ans2d))
# Compare the estimates of the models with parameters -10 and -50:
round(cbind(coef(ans2b1, shortenNames=FALSE), coef(ans2c1), coef(ans2d1)), 3)
# and the standard errors:
round(cbind(ans2b1$se , ans2c1$se, ans2d1$se), 3)
# How do the differences between the estimates
# compare to the average standard errors?
round((coef(ans2b1) - coef(ans2c1))/((ans2b1$se + ans2c1$se + ans2d1$se)/3), 3)
round((coef(ans2c1) - coef(ans2d1))/((ans2b1$se + ans2c1$se + ans2d1$se)/3), 3)
# All relative differences are less than 0.10, which is very small.
# This confirms that, for plausible results with this specification,
# the estimate for the parameter of the out-isolate effect
# could be -5 or more strongly negative, but not 0.


###############################################################################
###############################################################################