---
title: "Basic RSiena"
author: "Tom Snijders"
date: "March 16, 2021"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{BasicRSiena}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

This is an example of a basic sequence of commands
for estimating a basic model by function `siena07()` of `{RSiena}`. 
With a lot of use of help pages; this can be skipped as you like. 
Note that lines starting with # are comment lines, not commands. 
Sometimes in the following, some R commands are given preceded by #; 
this is to illustrate what could be done without executing it now. 

What is your current working directory?

```{r} 
#getwd()
```

If you wish it to be different, change it by

```{r} 
# setwd()
```

# Define data sets

If you have internet access, you can download the data
from the Siena website ("Data sets" tab)
http://www.stats.ox.ac.uk/~snijders/siena/s50_data.zip
and unzip it in your working directory. 
The data description is at
http://www.stats.ox.ac.uk/~snijders/siena/s50_data.htm

Then you can read the data files by the commands
(this can be replaced by using the internal data set, see below)

```{r}
# friend.data.w1 <- as.matrix(read.table("s50-network1.dat"))
# friend.data.w2 <- as.matrix(read.table("s50-network2.dat"))
# friend.data.w3 <- as.matrix(read.table("s50-network3.dat"))
# drink <- as.matrix(read.table("s50-alcohol.dat"))
# smoke <- as.matrix(read.table("s50-smoke.dat"))
```

But without internet access, the data can be obtained
from within RSiena (see below), because this is an internal data set.

```{r} 
library(RSiena)
```

With
```{r} 
?RSiena
```
you get to see the basic operation.

Now we use the internally available s50 data set.

Look at its description:

```{r} 
?s50
```
3 waves, 50 actors  
Look at the start and end of the first wave matrix:
```{r} 
head(s501)
tail(s501)
```
and at the alcohol variable
```{r} 
s50a
```

Now define the objects with the same names as above
(this step is superfluous if you read the data already).
```{r} 
        friend.data.w1 <- s501
        friend.data.w2 <- s502
        friend.data.w3 <- s503
        drink <- s50a
        smoke <- s50s
```

Now the data must be given the specific roles of variables
in an RSiena analysis.

## Dependent variable

Look at the help page:
```{r} 
?sienaDependent
```
First create a 50 * 50 * 3 array composed of the 3 adjacency matrices
```{r} 
friendshipData <- array( c( friend.data.w1, friend.data.w2, friend.data.w3 ),
           dim = c( 50, 50, 3 ) )
```
and next give this the role of the dependent variable:
```{r} 
	friendship <- sienaDependent(friendshipData)
```

What did we construct?

```{r}
friendship
```

We also must prepare the objects that will be the explanatory variables.

## Actor covariates

We use smoking for wave 1 as a constant actor covariate:

```{r} 
smoke1 <- coCovar( smoke[ , 1 ] )
```
A variable actor covariate is defined for drinking:
```{r} 
alcohol <- varCovar( drink )
```
(That for smoke we only use the first wave is purely for the purpose 
of illustrating constant actor covariates.)

#

Put the variables together in the data set for analysis

```{r} 
?sienaDataCreate
mydata <- sienaDataCreate( friendship, smoke1, alcohol )
```
and check what we have
```{r} 
mydata
```

You can get an outline of the data set with some basic descriptives from
```{r} 
print01Report( mydata, modelname="s50")
```
This creates an output file s50.txt in your current working directory. 
Note that if an existing file has this name, it will be overwritten! 

# Specifying the model

## Effects object

For the model specification we need to create the so-called effects object:
```{r} 
myeff <- getEffects( mydata )
```
This contains all the effects that are available given the structure
of this data set. They can be seen from
```{r} 
effectsDocumentation(myeff)
```
For a precise description of all effects, see Chapter 12 in the RSiena manual.
A basic specification of the structural effects is as follows:
```{r} 
?includeEffects
myeff <- includeEffects( myeff, transTrip, cycle3)
```
and this is put together with some covariate effects:
```{r} 
myeff <- includeEffects( myeff, egoX, altX, simX, interaction1 = "alcohol" )
myeff <- includeEffects( myeff, simX, interaction1 = "smoke1" )
myeff
```

## Algorithm

Now create an object with algorithm settings.
There are a lot of options for the algorithm.
The most important ones are explained in Sections 6.3.3 and 6.4
of the Siena manual.
Often you can just accept the default choices,
but it is good to specify the name for the output file
(which you may replace by any name you prefer)

```{r}
?sienaAlgorithmCreate
myalgorithm <- sienaAlgorithmCreate( projname = 's50' )
```
The file s50.txt already exists (created by print01Report);
the following estimations will append new results to this file.

# Estimate parameters

The function siena07 is the workhorse of the package;
it carries out the estimation.
```{r}
?siena07
ans <- siena07( myalgorithm, data = mydata, effects = myeff)
ans
```

This gives results from a random starting point,
and these will allways be different - hopefully, in unimportant ways.
To use a fixed starting point, so that you will always get
the same results (if also the same vector of initial parameters was used),
use the "seed" parameter:

```{r}
# myalgorithm <- sienaAlgorithmCreate( projname = 's50', seed=435123 )
```

For checking convergence, look at the
'Overall maximum convergence ratio' mentioned above, under the parameter estimates. 

This can also be shown by requesting

```{r}
ans$tconv.max
```

If this is less than 0.25, convergence is good.
If convergence is inadequate, estimate once more,
using the result obtained as the "previous answer"
from which estimation continues:

```{r}
ans <- siena07( myalgorithm, data = mydata, effects = myeff, prevAns=ans)
ans
```
If convergence is good, you can look at the estimates.
More extensive results are given by requesting 
```{r}
summary(ans)
```

Still more extensive results are given in the output file
s50.txt in the current directory.

Note that by putting an R command between parentheses (....),
the result will also be printed to the screen.
Next add the transitive reciprocated triplets effect,
an interaction between transitive triplets and reciprocity,
and estimate again, going on from the previous answer
(that is what `prevAns' means)

```{r}
(myeff <- includeEffects( myeff, transRecTrip))
(ans1 <- siena07( myalgorithm, data = mydata, effects = myeff, prevAns=ans))
```
If necessary, repeat the estimation, starting from the new result:
```{r}
(ans1 <- siena07( myalgorithm, data = mydata, effects = myeff, prevAns=ans1))
```

This might still not have an overall maximum convergence ratio
less than 0.25. If not, you could go on once more.

Inspect the file s50.txt in your working directory
and understand the meaning of its contents.

Parameters can be tested by dividing the estimate by its standard error,
and referring this to a standard normal distribution
(with the usual critical value of 1.96, rounded to 2.0).

This can be done also in the following way, e.g., for homophily
with respect to alcohol consumption:

```{r}
?Multipar.RSiena
Multipar.RSiena(ans1, 9)
```

Multi-parameter test are also possible.
To have a joint test of the three effects of alcohol:
```{r}
?Multipar.RSiena
Multipar.RSiena(ans1, 7:9)
```

Focusing on alcohol similarity, the effect is significant;
diluting the effects of alcohol by also considering ego and alter,
the three effects simultaneously are not significant.

# Assignment 1

## 1a.

Drop the effect of smoke1 similarity and estimate the model again.
Do this by the function setEffects() using the <<include>> parameter.
Give the changed effects object and the new answer object new names,
such as effects1 and ans1, to distinguish them.

## 1b.

Change the three effects of alcohol to the single effect
of alcohol similarity, and estimate again.

# Networks and behavior study

Now we redefine the role of alcohol drinking
as a dependent behaviour variable.

Once again, look at the help file
```{r}
?sienaDependent
```
now paying special attention to the 'type' parameter.
```{r}
drinking <- sienaDependent( drink, type = "behavior" )
```

Put the variables together in the data set for analysis:
```{r}
NBdata <- sienaDataCreate( friendship, smoke1, drinking )
NBdata
```

Define the effects object for this data set,
and give the model specification:
```{r}
NBeff <- getEffects( NBdata )
effectsDocumentation(NBeff)
NBeff <- includeEffects( NBeff, transTrip, transRecTrip )
NBeff <- includeEffects( NBeff, egoX, egoSqX, altX, altSqX, diffSqX,
                         interaction1 = "drinking" )
NBeff <- includeEffects( NBeff, egoX, altX, simX, interaction1 = "smoke1" )
NBeff
```
For including effects also for the dependent behaviour variable, see:
```{r}
?includeEffects
```
We now have to mention the name of the dependent behavior variable;
it did not need to be mentioned before, because the default name
is the first in the data set, which is friendship.
```{r}
NBeff <- includeEffects( NBeff, avAlt, name="drinking",
                         interaction1 = "friendship" )
NBeff
```
Define an algorithm with a new project name
```{r}
myalgorithm1 <- sienaAlgorithmCreate( projname = 's50_NB' )

```
Estimate again, using the second algorithm right from the start:
```{r}
NBans <- siena07( myalgorithm1, data = NBdata, effects = NBeff)
```
You may improve convergence (considering the overall maximum
convergence ratio) by repeated estimation in the same way as above.

Look at results:
```{r}
NBans
```
Make a nicer listing of the results:
```{r}
siena.table(NBans, type="html", sig=TRUE)
```

This produces an html file in your working directory; siena.table can also produce a LaTeX file.

# Assignment 2

## 2a.

Replace the average alter effect by average similarity (avSim)
or total similarity (totSim) and estimate the model again.

## 2b.

Add the effect of smoking on drinking and estimate again.

# Assignment 3

Read Sections 13.3 and 13.4 of the Siena Manual, download scripts
`SelectionTables.R` and `InfluenceTables.R` from the [Siena website](http://www.stats.ox.ac.uk/~snijders/siena/),
and make plots of the selection table and influence table for drinking.
