###############################################################################
## File g160_Onemode_Twomode.r                                               ##
## Analyses the Glasgow 160 data                                             ##
## One-mode - two-mode coevolution                                           ##
## By Tom Snijders                                                           ##
## Version: July 17, 2026                                                    ##
###############################################################################

library(RSiena)


# The data can be downloaded as
download.file("https://www.stats.ox.ac.uk/~snijders/siena/Glasgow_data.zip",
                                               destfile="Glasgow_data.zip")
unzip("Glasgow_data.zip")

# Read data sets:
load("Glasgow-friendship.RData")  # friendship networks
load("Glasgow-demographic.RData") # for gender data
load("Glasgow-various.RData") # for covariate data
load("Glasgow-lifestyle.RData") # for lifestyle data
load("Glasgow-selections.RData") # For selections into 129 and 50 data

# What do we have now:
ls()
dim(leisure1)
str(leisure1,1)
dimnames(leisure1)[[2]]

# Item 15 is not an activity, and will not be used.

# at home: activities 1, 3, 7, 8
# organized activities: 5, 9, 12
# going out to activities: 4, 10, 11, 14
# unsupervised: 6
# bored: 15
# requires money: 4, 8, 9, 10, 11, 14

students <- as_nodeset_rsiena(160, nodeSetName="students")
activities <- as_nodeset_rsiena(13, nodeSetName="activities", names=dimnames(leisure1)[[2]][2:14])

table(leisure1)
table(leisure2)
table(leisure3)

# Explore:
leis1 <- leisure1
leis1[leis1 <= 2 ] <- 1
leis1[leis1 >= 3 ] <- 0
which(is.na(leisure1))
colSums(leis1, na.rm=TRUE)

# Recode to binary:
leis1a <- leisure1
leis1a[leis1a <= 1 ] <- 1
leis1a[leis1a >= 2 ] <- 0

colSums(leis1a, na.rm=TRUE)

# This shows that almost all listen daily to tapes and CDs.
# This seems not to differentiate.
# These activities differ naturally in whether one would do them
# daily or almost so, or less than that.
# Daily or almost: 1, 2, 3, 5, 6, 7, 8, 14

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

# Some exploration

tl1 <- apply(leisure1,2, tabulate)
tl2 <- apply(leisure2,2, tabulate)
tl3 <- apply(leisure3,2, tabulate)

colSums(tl1 + tl2 + tl3)
# Note: these are almost constant, so missings are not a problem

# For Table T21.1
xtable(t(tl1 + tl2 + tl3))

# new thresholds, depending on frequencies:

threshold <- c(1,1,1,2,1,1,1,1,2,2,3,3,1,3)
threshold

leis <- array(0, dim=c(160,14,3))
for (i in 1:160){for (j in 1:14) {leis[i,j,1] <- ifelse(leisure1[i,j] <= threshold[j],1,0)}}
for (i in 1:160){for (j in 1:14) {leis[i,j,2] <- ifelse(leisure2[i,j] <= threshold[j],1,0)}}
for (i in 1:160){for (j in 1:14) {leis[i,j,3] <- ifelse(leisure3[i,j] <= threshold[j],1,0)}}

dim(leis)
leisure <- as_dependent_rsiena(leis[,2:14,], nodeSet=c("students", "activities"))
# Note that 15, "I do nothing much (am bored)" is dropped;
# also, the first is dropped, because it is too common.
# This leads in function as_dependent_rsiena to a bipartite network even without mentioning it,
# because two node sets are mentioned.
leisure

table(money)
money[money<0] <- NA
money[money> 99 ] <- NA
table(money)


# Define non-centered gender variables, values 0-1:
girls <- as_covariate_rsiena(sex.F-1, nodeSet="students", centered=FALSE)
boys <- as_covariate_rsiena(2-sex.F, nodeSet="students", centered=FALSE)
ages <- as_covariate_rsiena(age, nodeSet="students")
p.mon <- as_covariate_rsiena(money, type='monadic', nodeSet="students")

# How are leisure activities related to gender?
colSums(girls * leisure[,,1], na.rm=TRUE) / sum(girls)
colSums(boys * leisure[,,1], na.rm=TRUE)/ sum(boys)
colSums(girls * leisure[,,1], na.rm=TRUE) / colSums(leisure[,,1], na.rm=TRUE)
colSums(girls * leisure[,,1], na.rm=TRUE) / colSums(leisure[,,1], na.rm=TRUE)
# Used for columns 'proportion F' in Table 21.1:
(gl <- apply(leisure,3,function(x){colSums(girls*x, na.rm=TRUE)/colSums(x,na.rm=TRUE)}))
matrix(round(rowMeans(gl),2),13,1)

# the following activities require money: 4, 8, 9, 10, 11, 14 (old order)
payment <- rep(0, 13)
payment[c(3, 7, 8,  9, 10, 13)] <- 1
paying <- as_covariate_rsiena(payment, nodeSet="activities", centered=FALSE)
notpaying <- as_covariate_rsiena(1-payment, nodeSet="activities", centered=FALSE)
pay.act <- as_covariate_rsiena(payment, nodeSet="activities", centered=TRUE)

# Time trend:
timetrend <- as_covariate_rsiena(matrix(c(0, 1), 160, 2, byrow = TRUE),
             type='monadic', nodeSet="students")
# Composition change:
comp.ch <- as_composition_file_rsiena("cc.txt", nodeSet ="students",
                            option=2)
							
# Recode valued friendship to binary friendship:
friendship.1[friendship.1==2] <- 1
friendship.2[friendship.2==2] <- 1
friendship.3[friendship.3==2] <- 1
# Keep structural zeros because of use of sienaGOF later on.

# Identify dependent network variable:
friendship <- sienaDependent(array(c(friendship.1, friendship.2,friendship.3),
						dim=c(160, 160, 3)), nodeSet="students")


(meanF <- apply(leis[sex.F==2,,], c(2,3), mean, na.rm=TRUE))
(meanM <- apply(leis[sex.F==1,,], c(2,3), mean, na.rm=TRUE))
cbind(dimnames(leisure1)[[2]][1:14], round(meanF,2), round(meanM,2))

# Conclusion with respect to gender differences:
# somewhat more done by F:   "I listen to tapes or CDs" ,  "I read comics, mags or books" ,
#               "I hang round in the streets", "I go to dance clubs or raves"
# much more done by M:  "I go to sport matches"
# somewhat more done by M:  "I take part in sports", "I play computer games"    ,
#   "I spend time on my hobby (eg art, an instrument)" ,    "I go to something like B.B., Guides or Scouts"

# Construct data set:
G160_leisure_data <- make_data_rsiena(leisure, girls, boys, ages,
        p.mon, timetrend, paying, notpaying, pay.act,
        comp.ch, nodeSets=list(students, activities))
G160_leisure_data

write_report(G160_leisure_data, outputName = "G160_leisure")

save.image(file="G160_leisure0.RData")

#############################################################################
# Next section: co-evolution one-mode --- two-mode

load(file="G160_leisure0.RData")

# Create 1mode-2mode data set:

G160_netleis_Data <- make_data_rsiena(friendship, leisure,girls, boys, ages,
        p.mon, timetrend,  paying, notpaying, pay.act,
        comp.ch, nodeSets=list(students, activities))
G160_netleis_Data

# Specify model:
G160_netleis_eff <- make_specification(G160_netleis_Data)

G160_netleis_eff <- set_effect(G160_netleis_eff, list(outAct, inPop, reciAct),
             depvar="friendship")
G160_netleis_eff <- set_effect(G160_netleis_eff, inAct, depvar="friendship")
G160_netleis_eff <- set_effect(G160_netleis_eff, egoX, depvar="friendship",
             covar1="girls")
G160_netleis_eff <- set_effect(G160_netleis_eff, altX, depvar="friendship",
             covar1="girls")
G160_netleis_eff <- set_effect(G160_netleis_eff, sameX, depvar="friendship",
             covar1="girls")
G160_netleis_eff <- set_effect(G160_netleis_eff, gwespFF, type="creation",
             depvar="friendship", parameter=69)
G160_netleis_eff <- set_effect(G160_netleis_eff, gwespFF, type="endow",
             depvar="friendship", parameter=69)

G160_netleis_eff <- set_effect(G160_netleis_eff, list(outAct, inPop),
             depvar="leisure")
G160_netleis_eff <- set_effect(G160_netleis_eff, outInAss, depvar="leisure",
             parameter=1)
G160_netleis_eff <- set_effect(G160_netleis_eff, cycle4, depvar="leisure",
             parameter=1)

# What do we have now:
G160_netleis_eff


save.image(file="G160netleis0.RData")

# First estimate a baseline model without inter-network dependencies:
# Keep out-degrees for friendship to <= 6, like in the data:
alg_model <- set_model_saom(MaxDegree=c(friendship=6,leisure=20))
alg_alg <- set_algorithm_saom(seed=12345)
alg_algr <- set_algorithm_saom(seed=123456, 
				nsub=1, firstg=0.05, n2start=3000, n3=5000)
alg_alg20 <- set_algorithm_saom(seed=1234567,
				nsub=1, firstg=0.01, n2start=3000, 
				n3=20000)
alg_out0 <- set_output_saom(lessMem=TRUE)
				
(G160_netleis_1  <- siena(data=G160_netleis_Data, effects=G160_netleis_eff,
             nbrNodes=8, control_model=alg_model, control_algo=alg_alg))
(G160_netleis_1  <- siena(data=G160_netleis_Data, effects=G160_netleis_eff,
             nbrNodes=8, control_model=alg_model, control_algo=alg_algr, 
			 control_out=alg_out0, prevAns=G160_netleis_1)) 
			 
save(G160_netleis_1, file="G160netleis1.RData")

# Add some multivariate effects;

G160_netleis_eff2 <- set_effect(G160_netleis_eff, outActIntn, depvar="leisure",
             covar1="friendship", parameter=1)
G160_netleis_eff2 <- set_effect(G160_netleis_eff2, outActIntn,
             depvar="friendship", covar1="leisure", parameter=1)
G160_netleis_eff2 <- set_effect(G160_netleis_eff2, inActIntn, depvar="leisure",
             covar1="friendship", parameter=1)
G160_netleis_eff2 <- set_effect(G160_netleis_eff2, outPopIntn,
             depvar="friendship", covar1="leisure", parameter=1)
G160_netleis_eff2 <- set_effect(G160_netleis_eff2, to, depvar="leisure",
             covar1="friendship", parameter=1)
G160_netleis_eff2 <- set_effect(G160_netleis_eff2, from, depvar="friendship",
             covar1="leisure", parameter=1)

(G160_netleis_2  <- siena(data=G160_netleis_Data, effects=G160_netleis_eff2,
             nbrNodes=8, control_model=alg_model, control_algo=alg_alg, 
			 prevAns=G160_netleis_1))

(G160_netleis_2  <- siena(data=G160_netleis_Data, effects=G160_netleis_eff2,
             nbrNodes=8, control_model=alg_model, control_algo=alg_algr, 
			 control_out=alg_out0,
			 prevAns=G160_netleis_2))
			 			 
save(G160_netleis_2, file="G160netleis2.RData")


# Add the gender interaction effects that we found in the two-mode-only analysis:
G160_netleis_eff3 <- set_effect(G160_netleis_eff2, egoX, depvar="leisure",
             covar1="girls")
G160_netleis_eff3 <- set_interaction(G160_netleis_eff3, list(egoX, outAct),
             depvar="leisure", covar1=c("girls", ""))
G160_netleis_eff3 <- set_interaction(G160_netleis_eff3, list(egoX, totInDist2),
             depvar="leisure", covar1=c("girls", "girls"))
G160_netleis_eff3 <- set_interaction(G160_netleis_eff3, list(egoX, totInDist2),
             depvar="leisure", covar1=c("boys", "boys"))
G160_netleis_eff3

(G160_netleis_3  <- siena(data=G160_netleis_Data, effects=G160_netleis_eff3,
             nbrNodes=8, prevAns=G160_netleis_2, control_model=alg_model,
             control_algo=alg_alg))
			 			 
(G160_netleis_3  <- siena(data=G160_netleis_Data, effects=G160_netleis_eff3,
             prevAns=G160_netleis_3, control_model=alg_model,
             control_algo=alg_algr))

G160_netleis_3

write_result(G160_netleis_2, sig=TRUE, d=4)
write_result(G160_netleis_3, sig=TRUE, d=4)

# Add another two-mode homophily effect:
G160_netleis_eff5 <- set_effect(G160_netleis_eff3, sameXCycle4,
             depvar="leisure", covar1="girls", parameter=1)
(G160_netleis_5  <- siena(data=G160_netleis_Data, effects=G160_netleis_eff5,
             nbrNodes=5, prevAns=G160_netleis_3, control_model=alg_model,
             control_algo=alg_alg))
(G160_netleis_5  <- siena(data=G160_netleis_Data, effects=G160_netleis_eff5,
             nbrNodes=5, prevAns=G160_netleis_5, control_model=alg_model,
             control_algo=alg_algr))

# This is the model reported in Onemode_twomode_s.pdf.