###################################################### # Lab_coevolution.R # # # # R script for the assignment on co-evolution # # of alcohol consumption and friendship. # # Written by Christian Steglich and Tom Snijders. # # Version August 14, 2020 # # # ###################################################### # load RSiena package: library(RSiena) # This analysis uses the s50 data, # an excerpt of which is included in RSiena. # For a description, see # http://www.stats.ox.ac.uk/~snijders/siena/s50_data.htm # See short description on the help page: ?s50 # Identify dependent network variable: friendship <- sienaDependent(array(c(s501, s502, s503), dim=c(50, 50, 3))) # Identify dependent behavior variable: drinking <- sienaDependent(s50a, type="behavior") # Bind data together for Siena analysis: CoEvolutionData <- sienaDataCreate(friendship, drinking) # What do we have: CoEvolutionData # Create effects object for model specification: CoEvolutionEffects <- getEffects(CoEvolutionData) # First write descriptive results to file: print01Report(CoEvolutionData, modelname = 'CoEvol_init') # Look at this file! # Specify the model according to the assignment. # First the structural effects that are not yet included by default: CoEvolutionEffects <- includeEffects(CoEvolutionEffects, transTrip, transRecTrip) CoEvolutionEffects <- includeEffects(CoEvolutionEffects, inPop, outAct, inAct, reciAct) CoEvolutionEffects # Now sender, receiver and homophily effects for friendship formation: CoEvolutionEffects <- includeEffects(CoEvolutionEffects, egoX,altX,simX,interaction1="drinking") CoEvolutionEffects # Now the assimilation effect for drinking: CoEvolutionEffects <- includeEffects(CoEvolutionEffects, name="drinking",avSim,interaction1="friendship") CoEvolutionEffects # Note that you need to additionally specify 'name="drinking"' because # the default for 'name' is the network variable (here "friendship"). # Now create an algorithm object: CoEvolutionAlgo <- sienaAlgorithmCreate(projname = 'CoEvol_results') # Estimate the model: CoEvolutionResults <- siena07(CoEvolutionAlgo, data=CoEvolutionData, effects=CoEvolutionEffects) CoEvolutionResults # If the t-ratios for convergence for the non-fixed effects # are not satisfactorily small (all less than 0.1 in absolute value, # overall maximum convergence t-ratio less than 0.25), # run the estimation again. # We used different settings of the algorithm; see Section 6.3.3 of the manual. CoEvolutionAlgo.plus <- sienaAlgorithmCreate(projname = 'CoEvol_results', nsub=1, n2start=3000, n3=2000) (CoEvolutionResults <- siena07(CoEvolutionAlgo.plus, data=CoEvolutionData, effects=CoEvolutionEffects, prevAns=CoEvolutionResults)) # Test the specific parameter of interest: ?Wald.RSiena Multipar.RSiena(CoEvolutionResults, 18) # For a one-sided test, if the parameter is positive, # the p-value can be divided by 2. # Replace Average Similarity by Total Similarity: CoEvolutionEffects2 <- includeEffects(CoEvolutionEffects,include=FALSE, name="drinking",avSim,interaction1="friendship") CoEvolutionEffects2 <- includeEffects(CoEvolutionEffects2, name="drinking",totSim,interaction1="friendship") CoEvolutionEffects2 CoEvolutionResults2 <- siena07(CoEvolutionAlgo, data=CoEvolutionData, effects=CoEvolutionEffects2) CoEvolutionResults2 # If convergence is not adequate, run the estimation again: (CoEvolutionResults2 <- siena07(CoEvolutionAlgo.plus, data=CoEvolutionData, effects=CoEvolutionEffects2, prevAns=CoEvolutionResults2)) Multipar.RSiena(CoEvolutionResults2, 18) ################################################## # Explore 5-parameter model; for efficiency purposes, # use only the network as the dependent variable: drink <- varCovar(s50a) NetworkData <- sienaDataCreate(friendship, drink) NetworkEffects <- getEffects(NetworkData) NetworkAlgo <- sienaAlgorithmCreate(projname = 'Network_results') NetworkEffects <- includeEffects(NetworkEffects, transTrip, transRecTrip) NetworkEffects <- includeEffects(NetworkEffects, inPop, outAct, inAct, reciAct) NetworkEffects <- includeEffects(NetworkEffects, egoX, altX, egoSqX, altSqX, diffSqX, interaction1="drink") NetworkEffects (NetworkResults <- siena07(NetworkAlgo, data=NetworkData, effects=NetworkEffects)) # If convergence is not adequate, run the estimation again: #(NetworkResults <- siena07(NetworkAlgo, data=NetworkData, effects=NetworkEffects, # prevAns=NetworkResults)) # The effects of drinking alter and drinking squared alter seem not to matter: Multipar.RSiena(NetworkResults, 9:10) NetworkEffects <- includeEffects(NetworkEffects, altX, altSqX, interaction1="drink", include=FALSE) (NetworkResults <- siena07(NetworkAlgo, data=NetworkData, effects=NetworkEffects)) NetworkResults0 <- NetworkResults NetworkEffects <- includeEffects(NetworkEffects, egoX, altX, egoSqX, altSqX, egoXaltX, interaction1="drink") NetworkEffects (NetworkResults <- siena07(NetworkAlgo, data=NetworkData, effects=NetworkEffects)) # If convergence is not adequate, run the estimation again: # (NetworkResults <- siena07(NetworkAlgo, data=NetworkData, effects=NetworkEffects, # prevAns=NetworkResults)) # It seems we may drop the effect of drinking squared ego: Multipar.RSiena(NetworkResults, 12) # With this small data set, the decision we make about how to go on # is quite uncertain. We just drop this one effect. NetworkEffects <- includeEffects(NetworkEffects, egoSqX, interaction1="drink", include=FALSE) (NetworkResults <- siena07(NetworkAlgo, data=NetworkData, effects=NetworkEffects)) # If convergence is not adequate, run the estimation again: #(NetworkResults <- siena07(NetworkAlgo, data=NetworkData, effects=NetworkEffects, # prevAns=NetworkResults)) # Drinking squared alter also may be dropped: Multipar.RSiena(NetworkResults, 10) NetworkEffects <- includeEffects(NetworkEffects, altSqX, interaction1="drink", include=FALSE) (NetworkResults <- siena07(NetworkAlgo, data=NetworkData, effects=NetworkEffects)) # If convergence is not adequate, run the estimation again: #(NetworkResults <- siena07(NetworkAlgo, data=NetworkData, effects=NetworkEffects, # prevAns=NetworkResults)) # This exploration leads simply to dropping the two squared effects... # Now go backto the co-evolution model CoEvolutionEffects3 <- getEffects(CoEvolutionData) CoEvolutionEffects3 <- includeEffects(CoEvolutionEffects3, transTrip, transRecTrip) CoEvolutionEffects3 <- includeEffects(CoEvolutionEffects3, inPop, outAct, inAct, reciAct) CoEvolutionEffects3 <- includeEffects(CoEvolutionEffects3, egoX, altX, egoXaltX, interaction1="drinking") CoEvolutionEffects3 <- includeEffects(CoEvolutionEffects3, name="drinking",avAlt,interaction1="friendship") CoEvolutionEffects3 CoEvolutionResults3 <- siena07(CoEvolutionAlgo, data=CoEvolutionData, effects=CoEvolutionEffects3) # If convergence is not adequate, run the estimation again: CoEvolutionResults3 <- siena07(CoEvolutionAlgo.plus, data=CoEvolutionData, effects=CoEvolutionEffects3, prevAns=CoEvolutionResults3) CoEvolutionResults3 Multipar.RSiena(CoEvolutionResults3, 18) CoEvolutionEffects4 <- setEffect(CoEvolutionEffects3,avAlt,include=FALSE, name="drinking",interaction1="friendship") CoEvolutionEffects4 <- includeEffects(CoEvolutionEffects4, name="drinking",totAlt,interaction1="friendship") # Also specify to get the score-type tests for # the indegree and outdegree effects on drinking: CoEvolutionEffects4 <- includeEffects(CoEvolutionEffects4,indeg,outdeg,fix=TRUE,test=TRUE, name="drinking",interaction1="friendship") CoEvolutionEffects4 (CoEvolutionResults4 <- siena07(CoEvolutionAlgo, data=CoEvolutionData, effects=CoEvolutionEffects4)) # If convergence is not adequate, run the estimation again: (CoEvolutionResults4 <- siena07(CoEvolutionAlgo.plus, data=CoEvolutionData, effects=CoEvolutionEffects4, prevAns=CoEvolutionResults4)) Multipar.RSiena(CoEvolutionResults4, 20) # To see the results of the score-type tests for # the indegree and outdegree effects on drinking: score.Test(CoEvolutionResults4,18) score.Test(CoEvolutionResults4,19) save.image('Lab_coevol.RData') # The test for the outdegree effect on drinking comes close to significance, # the indegree effect is clearly not significant. # We can conclude at this point that all four specifications of # social influence result in a positive estimated effect, # with two-sided p-values all between 0.05 and 0.10, # and one-sided p-values correspondingly between 0.025 and 0.05. # This data set is too small to allow simultaneous estimation # of all four influence effects and do a backward elimination. # Choosing the most significant result is methodologically not permitted, # because it amounts to chance capitalization: # the level of significance would not be respected. # The best procedure is to make a choice before seeing the data, # based on theory or prior experience with similar data sets. # The choice can also be made based on the fit of the model, # if results from sienaGOF() differentiate between the models. # This is not elaborated in the current script. # Since the p-values in all four models are between 0.05 and 0.10 # (this will also depend on the random simulations; # results will be more stable with n3=5000 in sienaAlgorithmCreate) # we can conclude that there is a tendency toward evidence # for social influence. # All estimation results are in the file CoEvol_results.txt.