###############################################################################
##################    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 18, 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")

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

# 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=algo))
# 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=algo,
                        prevAns=ans, thetaBound=100))
# This is explored using a score-type test:
effs2a <- set_effect(effs, outIso, fix=TRUE, test=TRUE)
(ans2a <- siena(school13e, effects=effs2a, control_algo=algo, prevAns=ans))
test_parameter(ans2a, tested=TRUE, method="score", onestep=TRUE)
# The score-type test shows that the effect is significant.
# Fixing the parameter at -50 leads to convergence
# (I repeated estimation to lower the value of the maximum convergence ratio):
effs2b <- set_effect(effs, outIso, fix=TRUE, test=TRUE, initialValue=-50)
(ans2b <- siena(school13e, effects=effs2b, control_algo=algo))
(ans2b <- siena(school13e, effects=effs2b, control_algo=algo, prevAns=ans2b))
# Fixing the parameter at -10 leads also to convergence:
effs2c <- set_effect(effs, outIso, fix=TRUE, test=TRUE, initialValue=-10)
(ans2c <- siena(school13e, effects=effs2c, control_algo=algo))
(ans2c <- siena(school13e, effects=effs2c, control_algo=algo, prevAns=ans2c))
# The result of ans2c can be used,
# because the  maximum convergence ratio is small enough,
# and the score-type test for the value -10 is not significant.
# The score-type test for the value -50 is also not significant.
# This means that
# the parameter estimate for the out-isolate 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.

# Compare the estimates of the models with parameters -10 and -50:
round(cbind(coef(ans2b, shortenNames=FALSE), coef(ans2c)), 3)
# and the standard errors:
round(cbind(ans2b$se , ans2c$se), 3)
# How do the differences between the estimates
# compare to the average standard errors?
round((coef(ans2b) - coef(ans2c))/((ans2b$se + ans2c$se)/2), 3)
# All are less than 0.12. These differences are quite small, but not very small.

# Let us try to get more precise estimates.
# 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))
# Make the comparison again:
round(cbind(coef(ans2b1, shortenNames=FALSE), coef(ans2c1)), 3)
round(cbind(ans2b1$se , ans2c1$se), 3)
round((coef(ans2b1) - coef(ans2c1))/((ans2b1$se + ans2c1$se)/2), 3)
# Now all relative differences are less than 0.03, which is very small indeed.
# This confirms that, for plausible results with this specification, 
# the estimate for the parameter of the out-isolate effect 
# could be -10 as well as -50, but not 0.


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