###############################################################################
##################    RscriptOnestep.R           ##############################
###############################################################################
#
# This is an R script to illustrate the use of the one-step estimator in RSiena.
# It requires RSiena version at least 1.6.7.
#
# R script written by Tom A.B. Snijders.
# Version July 13, 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,  gwespFF))
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.
# However, we want to get better precision.
# We use 2 parallel processes to speed up estimation.
algo05 <- set_algorithm_saom(seed=4937, nsub=1, n2start=2000, n3=5000)
(ans0 <- siena(mydata, effects=myeff, control_algo=algo05, prevAns=ans0,
                    nbrNodes=2))

# Suppose we wish to investigate further degree effects.
# Two effects are added to the model as fixed-and-tested effects:

myeff1 <- set_effect(myeff, inAct, parameter=1, fix=TRUE, test=TRUE)
myeff1 <- set_effect(myeff1, reciAct, 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)
(ans1d <- siena(mydata, effects=myeff, control_algo=algo0, prevAns=ans0))

# The print of ans1 shows that the score-type test is not significant.
# The MoM estimate of ans0 and the one-step estimate after ans0
# are almost the same, especially in view of the standard errors:
(theta0 <- estimate_onestep(ans0))
round(cbind(coef(ans0), theta0, ans0$se), 3)

# Look again at ans1:
ans1
# The print of ans1 shows that score-type test is significant.
# This is also shown by
test_parameter(ans1, method="score")
# and for the separate effects by
test_parameter(ans1, tested=6, method="score")
test_parameter(ans1, tested=7, method="score")
# Both indegree-activity and reciprocal degree-activity 
# are significant in the score-type test,
# which does not control either effect for the other.

# The one-step estimate for freeing both fixed parameters is
cbind((ones0 <- estimate_onestep(ans1, fixed=0)))
# and for freeing each individually they are
(ones6 <- estimate_onestep(ans1, fixed=7))
(ones7 <- estimate_onestep(ans1, fixed=6))
# Compare them with a bit of rounding to have a better view:
ans1_coef <- coef(ans1)
round(cbind(ans1_coef, ones6, ones7, ones0), 3)
# We see that for the other 6 parameters,
# the one-step estimates ones7, where reciprocal degree-activity is freed,
# are quite a lot closer to the one-step estimates ones0, where both are freed.

# Let us do a forward model selection procedure step by step,
# adding first those effects
# that are most strongly significant in the score-type test.
# From the score-type tests above, we see that we have to start
# with freeing reciprocal degree-activity.

myeff2 <- set_effect(myeff1, reciAct, parameter=1)

(ans2 <- siena(mydata, effects=myeff2, control_algo=myalgo, prevAns=ans0))
# Again we want to achieve higher precision:
(ans2 <- siena(mydata, effects=myeff2, control_algo=algo05, prevAns=ans2,
                    nbrNodes=2))
# Now the score-type test is no longer significant:
test_parameter(ans2, method="score")

# The forward stepwise model selection now can stop.

# Compare the MoM estimate in ans2,
# the one-step estimate after ans1 if only inAct is freed,
# the one-step estimate after ans1 if both effects are freed,
# and the MoM estimate of ans1:
round(cbind(ans2$theta,  estimate_onestep(ans1, fixed=6),
            estimate_onestep(ans1, fixed=0), ans1$theta), 3)
# Especially for the correlated effect outdegree activity,
# which is included in all models,
# and which loses its significance when going from ans1 to ans2,
# the difference is quite important.

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

# 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=10, 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))
(ans <- siena(school13e, effects=effs, control_algo=myalgo, prevAns=ans))
siena_table(ans, sig=TRUE)

# How well does the model fit?
(ans_d <- siena(school13e, effects=effs, control_algo=algo0, prevAns=ans, returnDeps=TRUE))
(gof_d <- test_gof(ans_d, OutdegreeDistribution, levls=0:9, cumulative=FALSE, varName="fri"))
# p = 0
plot(gof_d)
# We start improving the fit for the outdegrees equal to 0.

# 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,
                                        prevAns=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.
# (chi-squared = 4.90, not tremendously high, but 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 "keepUnchanged=0"):
estimate_onestep(ans2a, fixed=0)
# Compare with MoM estimate for fixed specification:
round(cbind(ans2a$theta, estimate_onestep(ans2a, fixed=0)), 3)
# Estimate for outIso is  -5.244, quite high;
# others remain at moderate values.
(effs2b <- update_theta(effs2a, ans2a, onestep=TRUE, keepUnchanged=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))
siena_table(ans2b, sig=TRUE)
test_parameter(ans2b, tested=TRUE, method="score", onestep=TRUE)
estimate_onestep(ans2b, fixed=0)
# Now the out-isolate parameter is not significantly different
# from its new value at  -5.2444.
# However, the one-step estimate for the out-isolate parameter
# now has changed to  -9.238.
# what happens if we follow the recipe one more?
(effs2c <- update_theta(effs2b, ans2b, onestep=TRUE, keepUnchanged=0))
(ans2c <- siena(school13e, effects=effs2c, control_algo=myalgo))
estimate_onestep(ans2c, fixed=0)
test_parameter(ans2c, tested=TRUE, method="score", onestep=TRUE)
siena_table(ans2c, sig=TRUE)
# Again the out-isolate parameter is not significantly different
# from its new value at  -9.238.
# The other parameters are quite similar for ans2b and ans2c.
round(cbind( ans2a$theta, ans2b$theta, ans2c$theta, ans2c$se), 3)
round(cbind( (ans2b$theta- ans2c$theta)/ ans2c$se), 3)

# One more step....:
(effs2d <- update_theta(effs2c, ans2c, onestep=TRUE, keepUnchanged=0))
(ans2d <- siena(school13e, effects=effs2d, control_algo=myalgo))
(ans2d <- siena(school13e, effects=effs2d, control_algo=myalgo, prevAns=ans2d))
estimate_onestep(ans2d, fixed=0)
# This gives a warning and remains the same for outIso.
# With other andom number seeds, it is possible
# that the estimation algorith will diverge;
# probably becaues of numerical under/overflow)
# because of taking exponentials of values that are very large
# in absolute value.

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

# It can now be concluded that the out-isolate parameter is
# significantly different from 0 (ans2a, with p=0.027)
# but not from  -5.2444,  -9.238, 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, ans2e 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 ans2e,
# it is better to obtain even smaller maximum convergence ratios.

# 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".
# To make a moderate time gain (with a factor of about 2)
# I'm using parallel proceesing here,
# with only 2 parallel processes; you could add more,
# if you want to, and if your machine supports it.
algo1 <- set_algorithm_saom(seed=9492, nsub=1, n2start=10000,
                            n3=10000, firstg=0.02)
(ans2b1 <- siena(school13e, effects=effs2b, control_algo=algo1, prevAns=ans2b,
                    nbrNodes=2))
(ans2c1 <- siena(school13e, effects=effs2c, control_algo=algo1, prevAns=ans2c,
                    nbrNodes=2))
(ans2d1 <- siena(school13e, effects=effs2d, control_algo=algo1, prevAns=ans2d,
                    nbrNodes=2))
(ans2e1 <- siena(school13e, effects=effs2d, control_algo=algo1, prevAns=ans2e,
                    nbrNodes=2))
# Compare the estimates of the models with parameters -10 and -50:
estimates <- round(cbind(coef(ans2b1, shortenNames=FALSE),
                    coef(ans2c1),  coef(ans2d1), coef(ans2e1)), 3)
colnames(estimates) <- c("onestep", "twostep",  "threestep", "-50")
print(estimates)
# and the standard errors:
ses <- round(cbind(ans2b1$se , ans2c1$se, ans2d1$se, ans2e1$se), 3)
print(round(cbind(estimates, rowMeans(ses)),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+ ans2e1$se)/4), 3)
round((coef(ans2c1) - coef(ans2d1))/((ans2b1$se + ans2c1$se +
                    ans2d1$se+ ans2e1$se)/4), 3)
round((coef(ans2d1) - coef(ans2e1))/((ans2b1$se + ans2c1$se +
                    ans2d1$se+ ans2e1$se)/4), 3)
# All relative differences are less than 0.08 in absolute value.
# All relative differences between the estimates for ans2c1, ans2d1, and ans2e1
# are less than 0.04, 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 -5 or more strongly negative, but not 0;
# the values for ansc1, ansd1, or anse2, which could be rounded to -10 or lower,
# are slightly better than those for ansb1.


(ans_e <- siena(school13e, effects=effs2d, control_algo=algo0, prevAns=ans2e1, returnDeps=TRUE))
(gof_e <- test_gof(ans_e, OutdegreeDistribution, levls=0:9, cumulative=FALSE, varName="fri"))
# p = 0.016
plot(gof_e)
# The fit still is not great, but better thn the earlier p=0.

save.image("OneStep3.RData")

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