University of Groningen
Modeling the dynamics of networks and continuous behavior
Niezink, Nynke Martina Dorende
IMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from it. Please check the document version below.
Document Version
Publisher's PDF, also known as Version of record
Publication date: 2018
Link to publication in University of Groningen/UMCG research database
Citation for published version (APA):
Niezink, N. M. D. (2018). Modeling the dynamics of networks and continuous behavior. University of Groningen.
Copyright
Other than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).
Take-down policy
If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim.
Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the number of authors shown on this cover page is limited to 10 maximum.
5
Standard error accuracy
5.1 Introduction
In this chapter, we study a problem arising from dependence in the context of the stochastic actor-oriented model. When parameters in this model are esti-mated in the method of moments framework, the standard way of estimating standard errors is by the procedure proposed by Schweinberger and Snijders (2007). In some cases, repeated estimations of the standard errors have been observed to strongly di↵er. Standard error estimates can be exceptionally large. Because of the scaling of the e↵ects in stochastic actor-oriented models, standard errors usually range roughly between 0.05 and 2. However, inflated standard errrors can be as large as ten times the estimated parameter, and occur seem-ingly at random, depending on the particular Monte Carlo simulations used for their estimation (Schweinberger and Snijders, 2007; Ripley et al., 2018). A re-estimation of the model often results in smaller standard errors. This haphazard behavior occurs especially for models including parameters that are hard to es-timate for the data set under study (what is meant by ‘hard to eses-timate’ will be explicated later in the chapter). This has puzzled researchers applying the method and has led to worries about the accuracy of standard error estimates. Stochastic actor-oriented models are used to study the social mechanisms that drive network dynamics, and the interdependent dynamics of social networks and the attributes (e.g., performance, attitudes) of network actors (Snijders, 2001; Snijders et al., 2007). Reciprocity and popularity are two examples of so-cial mechanisms that can play a role in network dynamics. If an unreciprocated friendship tie of a popular person becomes reciprocated, this can be a indica-tion of reciprocity, but also one of popularity. This simple example illustrates the dependencies that give rise to the problem sketched above. Di↵erent social mechanisms sometimes lead to similar outcomes. A classical example of this –
on a more aggregate level – is network autocorrelation, a measure of the simi-larity of actors who are adjacent in a social network. Network autocorrelation can be the consequence of both peer selection based on similarity and of social influence (Steglich et al., 2010).
The problem sketched above may remind one of the (multi)collinearity prob-lem in multiple regression. This probprob-lem has been extensively studied for the regression model, and its extensions, like the multilevel regression or the struc-tural equation model (e.g., Shieh and Fouladi, 2003; Jagpal, 1982). Dependence among predictor variables makes the estimates of regression coefficients sensi-tive to small changes in the model or the data. It also leads to inflated standard errors of the estimates. Increased standard errors may result in the failure to reject false null hypotheses (type II errors). Moreover, regression estimates are sensitive to model misspecification when there is collinearity (Winship and Western, 2016).
Possibly even more so than in regression models, in stochastic actor-oriented models it is neither possible nor desirable to get rid of the dependencies. In the first example discussed above, both reciprocity and popularity are important predictors of friendship formation and we cannot just discard either of these mechanisms. In the second example, selection and influence processes are often the main focus of a study.
This chapter presents the problem of inflated standard errors in the stochastic actor oriented model using concepts borrowed from collinearity in multiple linear regression. We adapt a diagnostic developed in the context of the latter to a diagnostic for standard error inflation. In Section 5.2, we discuss how standard errors are estimated for parameters estimated using the method of moments and Monte Carlo simulations. In Section 5.3, we discuss collinearity in linear regression and the condition number as a collinearity diagnostic. We show how this diagnostic can be applied for the estimation of standard errors in the stochastic actor-oriented model. Section 5.5 explores the e↵ect a particular Monte Carlo simulation may have on standard error estimates. Section 5.6 presents alternative standard error estimators and assesses their performance. Section 5.7 concludes with possible directions for future research and practical recommendations.
Sections 5.2 and 5.3 are illustrated by the simulation example presented here-after, in Section 5.1.1. Sections 5.5 and 5.6 are illustrated by the empirical example introduced in Section 5.4.
5.1.1 Simulation example
We present a simple example to illustrate the way in which dependencies can be a nuisance for standard error estimation. The example does not reflect the situation discussed in the introduction, where standard errors are highly inflated. Nevertheless, the study provides a good idea of possible characteristics of problematic cases.
Figure 5.1 (right) shows a network X2 among 20 actors consisting of 100 ties, randomly selected such that 70 ties are reciprocated. Each actor is part of a group of 10 actors, but group membership and the existence of ties are unrelated. Network X1 is created by the random removal of 20 of the ties in X2 and the random addition of 20 previously non-existing ties (Figure 5.1, left).
Figure 5.1: Simulated networks X1 and X2. Black indicates reciprocated and gray non-reciprocated ties. Vertex color represents group membership.
We study the process of network evolution from X1to X2in four models. Model 1 only includes a density and a reciprocity e↵ect. These e↵ects reflect processes that we know to have played a role in the formation of X2. Model 2 adds the e↵ect of being part of the same group. Model 3 instead adds the outdegree activity e↵ect, the e↵ect of having outgoing ties on the tendency to form ties. Neither e↵ect has played a systematic role in the creation of the data. Model 4 combines all e↵ects. The objective function corresponding to model 4 is
fi( , x) = 1 X j xij+ 2 X j xijxji+ 3 X j xijI(vi= vj) + 4 X j xij 2 , (5.1)
where vi2 {0, 1} denotes the group affiliation of actor i.
Table 5.1 shows parameter estimates and standard errors for the four models. The results for Model 1 show that tie formation is generally unlikely, unless it is an act of reciprocation. This reflects the process behind the generation of X1 and X2 well. The additional same group e↵ect (Model 2) does not change the estimates. If we instead add the outdegree activity e↵ect (Model 3) both
Table 5.1: Four stochastic actor-oriented models estimated on the longitu-dinal network data depicted in Figure 5.1 (dens=density, recip=reciprocity, same=same group, outd=outdegree activity).
Model 1 Model 2 Model 3 Model 4
rate 2.51 (0.43) 2.50 (0.42) 2.49 (0.43) 2.48 (0.43) dens 1 1.34 (0.22) 1.32 (0.28) 0.52 (0.75) 0.53 (0.80) recip 2 2.23 (0.39) 2.24 (0.39) 2.26 (0.40) 2.27 (0.40)
same 3 0.042 (0.34) 0.078 (0.33)
outd 4 0.083 (0.072) 0.069 (0.074)
the estimate and standard error of the density parameter change considerably, while the reciprocity parameter remains una↵ected. Adding both e↵ects simul-taneously does not change this picture.
In conclusion, a redundant e↵ect that is highly correlated with other e↵ects in the model may obscure the true nature of the e↵ects it is correlated with. A redundant e↵ect that is uncorrelated does not have this e↵ect.
5.2 Standard error estimation
Before we discuss standard error estimation, we first briefly review the method of moments framework to estimate the parameters ✓ in the stochastic actor-oriented model. In this framework, for each parameter a statistic is selected that can be evaluated on a given (simulated or observed) data set and that is sensitive to changes in the parameter. Let S denote the vector of all statistics for the parameters ✓. The method of moments estimator ˆ✓ is given by the solution to the moment equation E✓(S) = sobs, where sobs denotes the value of the statistics in the observed data. The moment estimator ˆ✓ cannot be determined analytically, and is estimated using a multivariate Robbins-Monro stochastic approximation algorithm (Robbins and Monro, 1951; Kushner and Yin, 2003).
For estimating the standard errors of ˆ✓, it is inconvenient to use bootstrap (or resampling) methods, because each of the estimation runs required by such methods would be time-consuming. Therefore, Snijders (1996) proposed to approximate the covariance matrix cov(ˆ✓) of the moment estimator by
cov(ˆ✓)⇡ ˆD 1⌃ ( ˆˆ D 1)>, (5.2) where the Jacobian matrix of the statistics D = @ E(S)/@✓ and the the co-variance matrix of the statistics ⌃ = cov(S) are evaluated at the estimator ˆ
The general idea behind the approximation can be sketched as follows. Consider the first-order Taylor expansion of the statistics S around the true parameter ✓,
S(ˆ✓) = S(✓) + D (ˆ✓ ✓) +O(kˆ✓ ✓k2), (5.3) where the Jacobian D is evaluated at ✓. Taking covariances at both sides yields cov(S(ˆ✓))⇡ cov(D ˆ✓) = D cov(ˆ✓)D>. (5.4) Assuming that D is a non-singular matrix, we can rearrange this expression to obtain
cov(ˆ✓)⇡ D 1cov(S(ˆ✓))(D>) 1. (5.5) Given a non-singular estimator ˆD, the covariance matrix cov(ˆ✓) can be es-timated by (5.2). If a model contains e↵ects that are similarly sensitive to changes in the data, D will be close to being singular. We will formally define ‘nearness to singularity’ in Section 5.3.
5.2.1 Monte Carlo estimation
Neither the Jacobian matrix D nor the covariance matrix of the statistics ⌃ can be expressed in closed form. Therefore, they are determined using Monte Carlo estimation, i.e., based on a large number of data sets simulated under the estimated model with estimated parameters ˆ✓. Monte Carlo estimation of ⌃ is trivial. It is the sample covariance matrix of the values of the statistics S in the simulated data sets.
The estimation of the Jacobian D requires more e↵ort. Snijders (2001) esti-mated the Jacobian by Monte Carlo methods based on finite di↵erences with common random numbers. Schweinberger and Snijders (2007) proposed three new estimators related to the score function method of derivative estimation, which improved on the finite di↵erences method both on theoretical and prac-tical grounds. The new estimators are unbiased and reduce the computation time by an order of the number of parameters. They are based on the idea that the Jacobian can be rewritten as
D =@S @✓ = E✓ ✓ S@ log p(Y ) @✓0 ◆ (5.6)
(Schweinberger and Snijders, 2007, Lemma 2). Here Y denotes the complete data corresponding to a stochastic actor-oriented model analysis, that is, the observed data, the holding times of the Markov process and the sequence of all changes in the modeled time interval. The complete data likelihood is denoted by p(Y ).
Based on N Monte Carlo simulations of the data, Y1, . . . , YN, given the moment estimator ˆ✓, we can estimate (5.6) by
DI = 1 N N X i=1 Si @ log p(Yi) @✓0 = 1 N N X i=1 SiJi, (5.7) where Si denotes the statistics evaluated on simulated data Yi and Ji denotes the complete data score @ log p(Yi)/@✓0. The complete data scores can be com-puted using the expression derived by Schweinberger and Snijders (2007) and in Chapter 3 of this dissertation (see Appendix B of Chapter 3). Using variance reduction methods based on control variates, Schweinberger and Snijders (2007) also proposed the estimator
DII= 1 N N X i=1 Si 1 N N X i=1 Si ! Ji. (5.8)
The variance of DIIis considerably smaller than that of DI.
The statistic S is the sum over statistics Sm computed separately based on the periods m between measurements at time tm and tm+1 (see Section 3.4). Therefore, we can rewrite the Jacobian (5.6) as
D = @(S 1+ . . . + SM 1) @✓ = E✓ M 1X m=1 Sm@ log p(Y m) @✓0 ! , (5.9)
where Ymis that part of the simulated data related to period m. The estimator based on a similar decomposition of (5.8),
1 N N X i=1 M 1X m=1 ⇣ Sim 1 N N X i=1 Sim⌘Jim ! , (5.10)
is nowadays the default in the estimation of standard errors in the stochastic actor-oriented model (Snijders, 2017). We will refer to this estimator as the score function estimator.
5.2.2 Application to the simulation example
To understand the results in Section 5.1.1 we consider the matrices ˆ⌃ (upper triangle entries contain correlations) and ˆD corresponding to the first three models: ˆ ⌃1= 38.4 0.75 31.1 44.3 ! ˆ ⌃2= 0 B @ 40.6 0.78 0.70 34.3 47.3 0.55 19.6 16.6 19.4 1 C A ˆ⌃3= 0 B @ 30.7 0.72 0.97 24.5 37.2 0.72 310.7 251.5 3313.7 1 C A , ˆ D1= 34.8 13.3 26.2 21.4 ! ˆ D2= 0 B @ 36.9 15.1 18.1 29.4 23.0 14.6 18.1 7.5 18.3 1 C A ˆD3= 0 B @ 28.2 10.7 283.9 21.5 18.7 222.1 285.6 111.5 3061.2 1 C A .
The covariance matrices ⌃ show that, as expected, the correlation between the simulated statistics of the density parameter 1 and those of the outdegree parameter 3 is very high. This is also reflected in the linear dependencies among the rows of the Jacobian matrix ˆD3, which is especially clear if we scale the rows of ˆD to have unit length:
ˆ Ds 1= 0.93 0.36 0.77 0.63 ! ˆ Ds 2= 0 B @ 0.84 0.34 0.41 0.73 0.57 0.36 0.68 0.28 0.68 1 C A ˆDs 3= 0 B @ 0.099 0.039 0.994 0.096 0.083 0.992 0.093 0.036 0.995 1 C A .
The first and last row of ˆDs
3are very similar, i.e., the statistics of 1and 3are similarly sensitive to changes in the parameters. Jacobian matrices ˆDs
1 and ˆDs2 do not exhibit such a pattern.
5.3 Diagnosing standard error inflation
The rows of Jacobian matrix D are often approximately linearly dependent, as can be the predictors in a linear regression model. In this section we first discuss collinearity in linear regression and present the condition number of the design matrix (the matrix containing the values of all predictors for all subjects) as a collinearity diagnostic. We then show that the condition number of the Jacobian is a good diagnostic for standard error inflation. We illustrate this by application to the Jacobian matrices estimated in the previous section.
5.3.1 Collinearity in linear regression
Consider the standard linear regression model
y = 0+ 1x1+ . . . + kxk+ ✏, (5.11) with ✏ ⇠ N (0, 2). Suppose that n observations on the variables y, x
1, . . . , xk are available. These n realizations can be represented in matrix notation by
Y = X + e, (5.12)
where the random variables ✏1, . . . , ✏nin vector e are independent. Matrix X is called the design matrix. The maximum likelihood estimates of and 2 are given by ˆ = (X>X) 1X>Y (5.13) and ˆ2= 1 n(Y X ˆ) >(Y X ˆ). (5.14)
The covariance matrix of the parameters estimates ˆ is estimated by
cov( ˆ) = ˆ2(X>X) 1. (5.15) Matrix X>X is, up to a factor 1/n, the covariance matrix of the predictor variables. The standard errors of the parameter estimates ˆ are large if the predictor variables are nearly collinear.
Many procedures have been employed to detect collinearity in regression models (see, e.g., Belsley, Kuh, and Welsch, 1980). One way is by examining the correlation matrix of the predictor variables. However, while a high correlation between two variables can indeed point to a potential collinearity problem, the absence of high correlations does not ensure the absence of this problem. It is possible for three or four variables to be collinear while no pair of variables alone is highly correlated.
A diagnostic that is able to indicate near collinearity among two or more vari-ables is the condition number of design matrix X. If a matrix has a high condition number, this is an indication that calculations based on this matrix may yield difficulties and that the entries of its inverse may be large. The condi-tion number of a matrix with a perfect collinearity is infinite. In the following, we define the condition number of a matrix by considering its use in solving systems of linear equations.
5.3.2 The condition number
To calculate regression parameters we solve the equation (X>X) ˆ = X>Y . For their standard erors we solve (X>X) cov( ˆ) = ˆ2I. In general, if we are interested in solving a linear system of equations Az = b, with A a nonsingular square matrix, we can ask how much the solution vector z would change ( z) for a small change A in the elements of A (e.g., Lange, 2010, p.87). Let kzk represent a vector norm andkAk the induced matrix norm, that is,
kAk = sup kzk6=0
kAzk
kzk . (5.16)
It can be shown that
k zk
kz + zk kAkkA
1kk Ak
kAk. (5.17)
This inequality is sharp (Lange, 2010, p.90).
The condition number of A is defined as kAkkA 1k and is denoted by (A).1 It provides a bound on the relative change in the norm of the solution vector z
induced by a relative change in the norm of A. If we take vector b to be the kth unit vector, then z = A 1b contains the kth column of A 1. Therefore, inequal-ity (5.17) implies that the entries of A 1may be large, if A has a high condition number. If A is singular, it is said to have an infinite condition number. A matrix with a very large condition number is said to be ill-conditioned. If we let A = X>X, z = ˆ and b = X>Y , inequality (5.17) shows that the condition number of X>X is a good indicator of ˆ being sensitive to small changes in the data. Moreover, if X>X has a high condition number, the entries of (X>X) 1 and thus the entries of cov( ˆ) may be large.
In the above exposition we have not selected any particular norm. In fact, each matrix norm induces its own condition number. If we let kzk denote the Euclidean norm of z, the induced matrix norm is the spectral norm, i.e., kAk = max(A), the largest singular value of A. The condition number in this case is equal to
(A) = max(A)/ min(A). (5.18) Singular values are defined as follows. The singular value decomposition of an m⇥n matrix A is a factorization of the form UDV>, where U is an m⇥m unitary matrix (its conjugate transpose is also its inverse), D is an m⇥ n rectangular diagonal matrix with non-negative real numbers on the diagonal, and V is an n⇥ n unitary matrix. The diagonal entries of D are called the singular values of A and can also be computed as the square roots of the eigenvalues of A>A. Singular values are often used in collinearity diagnostics in linear regression, that are based on the ideas of conditioning presented in this section. Belsley et al. (1980, Appendix 3A) show that the condition number (X), as defined in (5.18), provides an upper bound on the sensitivity of the diagonal elements of (X>X) 1 to any element of X. This bound thus directly a↵ects the standard errors of ˆ.
For each exact linear dependency among its columns, matrix X has one zero singular value. Also near linear dependencies can be linked to the singular values of X. Let max(X)/ k(X) define the kth condition index of X. There are as many near linear dependencies in X as there are high condition indices (or, equivalently, low singular values) (Silvey, 1969). Based on numerical ex-periments, Belsley et al. (1980) indicate that weak dependencies are associated with condition indices around 5 or 10, and moderate to strong relations with condition indices of 30 to 100. See Silvey (1969) and Belsley et al. (1980, p.105) for a diagnostic procedure to identify the variables involved in a near depen-dency and to assess the degree to which that dependepen-dency degrades the variance of each of the regression estimates. The procedure uses the fact that the
esti-mated variance of regression coefficients can be decomposed into a sum of terms each of which is associated with a singular value of X.
5.3.3 Standard error inflation
In the previous section, we assessed the sensitivity of z in Az = b to changes in A and b. We here assess the sensitivity of cov(✓) in the equation
D cov(✓)D>= ⌃ (5.19)
to changes in ⌃ or D. The resulting inequalities can be applied to understand the sensitivity of cov(ˆ✓) in (5.2) to changes in ˆ⌃ or ˆD. The proofs of the following two propositions are given in Appendix A. The propositions hold for any submultiplicative matrix norm.
Proposition 1. Suppose an approximation of ⌃ by ⌃ + ⌃ would result in an error C in cov(✓). This error satisfies the inequality
k Ck
k cov(✓)k (D)(D >)k ⌃k
k⌃k . (5.20)
Proposition 2. Suppose an approximation of D by D + D would result in an error C in cov(✓). This error satisfies the inequality
k Ck k cov(✓) + Ck (D) k Dk kDk + (D >)k D>k kD>k + (D)(D >)k Dk kDk k D>k kD>k. (5.21)
For norms that satisfy kAk = kA>k, such as the spectral norm, the above inequalities simplify greatly. The inequalities (5.20) and (5.21) respectively reduce to k Ck k cov(✓)k (D) 2k ⌃k k⌃k (5.22) and k Ck k cov(✓) + Ck 2(D) k Dk kDk + (D) 2 ✓ k Dk kDk ◆2 . (5.23)
All these results are inequalities, and hence show the maximum possible error inflation. The actual inflation may be considerably smaller than this maximum. The inequalities motivate the use of the condition number of ˆD to measure the extent of potential standard error problems.
One may wonder whether it is possible to determine the extent to which near dependencies, indicated by high condition indices, degrade the variance of the parameters in a way analogous to the procedure discussed by Belsley et al. (1980,
p.105). A meaningful decomposition of D analogous to the decomposition of the variance of the regression coefficients, which was referred to in Section 5.3.2, however does not exist. The procedure to identify the variables involved in a near dependency discussed by Belsley et al. (1980) is therefore not applicable. However, a leave-one-out method, where statistics are left out one by one and the Jacobian and its condition number are re-computed, may yield insight into which e↵ect a↵ects the conditioning most.
Finally, as touched upon in the previous section, the condition number of D also indicates how sensitive the inverse of D is to changes in D. A proof of the following proposition is given in Appendix A.
Proposition 3. An approximation of D by D + D would result in an error in the inverse of D that satisfies the inequality
kD 1 (D + D) 1k
k(D + D) 1k (D) k Dk
kDk . (5.24)
As a consequence,kD 1k is sensitive to changes in D for D with a high condition number. Since the condition number is defined as kDkkD 1k, the condition number itself in this case is also sensitive to changes in D.
5.3.4 Application to the simulation example
The condition numbers, induced by the spectral norm, of the Jacobian matrices in the example in Section 5.1.1 are given below.
Model 1 Model 2 Model 3 Model 4 ( ˆD) 6.2 10.7 1895.5 2125.3
The condition numbers of Models 3 and 4, which include the outdegree activity e↵ect, are much larger than those of Models 1 and 2. The high condition numbers indicate that there is a near linear dependency among the rows of Jacobian ˆD in these models that was not present in Models 1 and 2. The high correlation between the density and outdegree activity statistics suggests that these e↵ects are involved in the near dependency leading to the increase in the condition number. The condition number of the Jacobian ˆD is a concise way to describe the existence of dependencies.
5.4 Empirical example: friendship and body mass index
In this section, we present a stochastic actor-oriented model of the friendship network and body mass index data studied in the empirical example in Chapter 3, focusing on the standard error estimates. This analysis is used as a running example in Sections 5.5 and 5.6.
In Chapter 3, we studied the co-evolution of friendship and body mass index (BMI), based on the friendship networks and the BMI of 156 high school stu-dents, measured at four occasions. Table 5.2 shows the estimation of Model 2, as specified in Section 3.4, based on only the first two waves of data. The parameter estimates in this model have converged: all convergence t-ratios are smaller than 0.1 in absolute value and the overall maximum convergence ratio of the model is 0.16.
Table 5.2 shows the standard errors based on 1000 and 50,000 Monte Carlo simulations. Some of the standard errors based on the 1000 simulations, such as the ones corresponding to the parameters related to BMI similarity, are somewhat large, but generally the size of these standard error estimates would not immediately raise concern. However, they di↵er considerably from the standard errors based on the 50,000 Monte Carlo simulations, both in size and in terms of conclusions drawn on significance testing. For example, the female alter e↵ect is not significant based on the N = 1000 simulations standard error estimate, but based on the N = 50,000 simulations estimate, it is.
The condition number of the Jacobian ˆD corresponding to the N = 1000 sim-ulations estimates is approximately 140,000, indicating that small deviations in ˆD or ˆ⌃ could be leading to errors in the estimated value of cov(ˆ✓). The condition number of the Jacobian corresponding to the N = 50,000 simulations estimates is approximately 93,000.
Table 5.2 only presents the standard error estimates based on N = 1000 and N = 50,000 simulations. In Section 5.4.1 we provide a more detailed picture by considering the evolution of the standard error estimates over 100,000 Monte Carlo simulations. Of course, these estimates are based on a particular draw of N simulations. A di↵erent draw would have resulted in di↵erent estimates. In Section 5.4.2, we therefore consider the distribution of the standard error estimates by bootstrap methods. Finally, in Section 5.4.3, we visualize the de-pendencies among the moment statistics corresponding to the model parameters in Table 5.2.
Table 5.2: Stochastic actor-oriented model of friendship and BMI dynamics. Standard errors are estimated based on 1000 and 50,000 Monte Carlo simula-tions.
N = 1000 N = 50,000
estimate (s.e.) (s.e.) Friendship dymamics rate 7.76 (1.15) (0.78) outdegree 3.57 (0.50)⇤ (0.25)⇤ reciprocity 3.41 (0.34)⇤ (0.32)⇤ transitivity (gwesp) 2.66 (0.24)⇤ (0.22)⇤ cyclicity (gwesp) 1.50 (0.72)⇤ (0.31)⇤ transitivity (gwesp)⇥ reciprocity 0.53 (0.37) (0.19)⇤ indegree popularity 0.070 (0.038) (0.030)⇤ outdegree activity 0.045 (0.017)⇤ (0.016)⇤
female alter 0.48 (0.42) (0.18)⇤
female ego 0.42 (0.21)⇤ (0.17)⇤
same gender 0.85 (0.56) (0.15)⇤
same home group 0.75 (0.25)⇤ (0.18)⇤
same home group⇥ reciprocity 0.79 (0.76) (0.33)⇤
bmi alter 0.005 (0.039) (0.023)
bmi ego 0.005 (0.026) (0.021)
bmi similarity 0.66 (1.27) (0.92)
bmi similarity⇥ reciprocity 1.49 (3.29) (1.58) Distress dymamics feedback a 0.016 (0.026) (0.020) intercept b0 0.33 (0.07)⇤ (0.07)⇤ error g 0.59 (0.03)⇤ (0.06)⇤ average alter b1 0.087 (0.048) (0.050) ⇤p-value < 0.05.
5.4.1 Convergence of the standard error estimates
The convergence of the standard errors is largely dependent on the convergence of the inverse of Jacobian estimator ˆD 1, and thus of the condition number of
ˆ
D. Figure 5.2a shows the evolution of the condition number of ˆD, plotted on a logarithmic scale, over the first 1000 Monte Carlo simulations. The value at sample k represents the condition number of ˆD based on the first k simulations. The condition number is consistently very high and has not converged yet. Note that the higher the true condition number of D, the more the condition number is sensitive to changes in D, and the more difficult it is to reach convergence.
● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●●● ●●●●●●●●●●●● ●●●●●●●● ●●●●●● 0 200 400 600 800 1000 5.0 5.2 5.4 5.6 5.8 6.0 6.2 samples 10log( condition n umber )
(a) Condition number. ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ●● ●●● ● ● ●● ●● ●●●● ●● ● ● ●● 0 200 400 600 800 1000 − 0.6 − 0.4 − 0.2 0.0 0.2 0.4 0.6 samples
10log( standard error )
(b) Female ego. ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ●●●●● ●● ● ●● ●●●● ●● ●●●● ●● ● ● ●● 0 200 400 600 800 1000 − 0.5 0.0 0.5 samples
10log( standard error )
(c) Female alter. ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ●● ●●●● ●●●● ●●●●●●● ● ●● ●●●● 0 200 400 600 800 1000 − 0.5 0.0 0.5 1.0 samples
10log( standard error )
(d) Reciprocity. ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ●●●●●●●●●●● ● ● ●●●●●●●● ●●●●● 0 200 400 600 800 1000 0.0 0.5 1.0 1.5 samples
10log( standard error )
(e) BMI similarity.
● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●●● ●●●●●● ● ●●●●●● ●●● ●●●●● 0 200 400 600 800 1000 0.5 1.0 1.5 2.0 samples
10log( standard error )
(f) Recip.⇥ BMI similarity.
Figure 5.2: The evolution of the condition number of the Jacobian, figure (a), and standard errors, figures (b) to (f), on N = 1000 Monte Carlo simulations. The red line indicates the standard error value corresponding to a p-value of exactly 0.05. The blue line indicates the N = 100,000 standard error estimate.
Figure 5.2 also shows the evolution of the standard errors of some parameters in the model, plotted on a logarithmic scale. The standard errors clearly have not converged yet; most of them are far from the blue lines corresponding to the N = 100,000 simulations standard error estimates. The red lines indicate the standard error values that correspond to a p-value of exactly 0.05. The female ego and female alter standard errors are still fluctuating around this line. Figure 5.3 shows the evolution of the condition number of the Jacobian and stan-dards errors on the 49,000 simulations following the first 1000 simulations. The standard errors contain fewer outliers and are no longer plotted on a logarithmic scale, which allows for better evaluation of their ‘true’ size. The standard errors corresponding to the female ego and female alter e↵ect are clearly converging to a value below the p = 0.05 line. The estimators of the BMI similarity stan-dard error ranges between 0.92 and 1.19 between N = 30,000 and N = 50,000. These estimates are all fairly large compared to the estimated BMI similarity parameter of 0.66.
● ● ●●● ●●●● ● ●● ●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●● 0 10000 20000 30000 40000 50000 5.00 5.05 5.10 5.15 samples 10log( condition n umber )
(a) Condition number. ● ● ● ●● ●● ● ● ● ● ●● ●● ●●●● ●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●● 0 10000 20000 30000 40000 50000 0.2 0.3 0.4 0.5 samples standard error (b) Female ego. ● ● ● ● ● ● ● ● ● ●● ●●●●●● ● ● ●●● ●●●●●●●●● ●●● ●●●●● ● ● ●●●● ●●●●● 0 10000 20000 30000 40000 50000 0.16 0.18 0.20 0.22 0.24 samples standard error (c) Female alter. ● ● ● ● ● ●● ● ● ●●● ●●● ●●● ●●●● ●● ●●●●●●●●●●● ●●●● ●●●●●●●●●●● 0 10000 20000 30000 40000 50000 0.4 0.6 0.8 1.0 1.2 samples standard error (d) Reciprocity. ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ●●●●● ● ●● ●● ● ● ● ● ●● ●● ● ● ●●● ● ●●●●● ●● ●● 0 10000 20000 30000 40000 50000 0.9 1.0 1.1 1.2 1.3 1.4 samples standard error
(e) BMI similarity.
● ●●● ● ●● ● ● ●●●●●●●●●●●●●●●●●●●●●●●●● ● ●●●●●●● ●●●●●●●● 0 10000 20000 30000 40000 50000 2.0 2.5 3.0 samples standard error
(f) Recip.⇥ BMI similarity.
Figure 5.3: The evolution of the condition number of the Jacobian, figure (a), and standard errors, figures (b) to (f), on N = 50,000 Monte Carlo simulations. The red line indicates the standard error value corresponding to a p-value of exactly 0.05. The blue line indicates the N = 100,000 standard error estimate.
50,000 simulations. Figure 5.4 shows the evolution of the condition number of the Jacobian and the standard errors based on 20,000 to 100,000 simulations. Between simulation N = 80,000 and 100,000, the condition number ranges between around 92,100 and 92,700. This range is considerably smaller than the ranges depicted in Figures 5.2a and 5.3a. The standard errors all seem to have converged after 100,000 simulations.
5.4.2 Bootstrap distribution
If we would have based our standard error estimates on a di↵erent draw of 1000 Monte Carlo simulations, we would have obtained di↵erent estimates in Table 5.2. These estimates could have resulted in di↵erent conclusions based on significance tests. This is generally true for estimates based on Monte Carlo methods. In this section, we explore the distribution of the standard errors using bootstrap sampling. By sampling with replacement N simulations from the first 50,000 simulations described in the previous section, we take a bootstrap sample
● ● ●● ●● ● ●●● ● ● ● ● ●● ●●●●● ●● ● ● ● ●●● ●● ●●●●● ● ● ● ● ●● ● ●●●●●●● ●●●● ● ●●●● ●● ●●● ●● ● ●●● ●● ●● ●●●● ●●
2e+04 4e+04 6e+04 8e+04 1e+05
4.965 4.970 4.975 4.980 4.985 samples 10log( condition n umber )
(a) Condition number. ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●●●●● ●●●●● ● ●●● ●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●
2e+04 4e+04 6e+04 8e+04 1e+05
0.185 0.195 0.205 0.215 samples standard error (b) Female ego. ● ● ●●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ●● ●● ●● ●●● ● ● ●●● ● ● ●● ● ●●●● ●●● ●●●●●●●●●●●●● ●●●●● ●
2e+04 4e+04 6e+04 8e+04 1e+05
0.165 0.170 0.175 samples standard error (c) Female alter. ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ●●● ● ●●●●● ● ●●● ● ● ● ● ● ●●●● ●●● ●●●●●●●●● ●● ●●●●●●●●●●●●●●●●●●●●●●●●●
2e+04 4e+04 6e+04 8e+04 1e+05
0.32 0.34 0.36 0.38 samples standard error (d) Reciprocity. ●● ● ●● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●●●● ● ● ●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●
2e+04 4e+04 6e+04 8e+04 1e+05
0.90 0.95 1.00 1.05 1.10 1.15 1.20 samples standard error
(e) BMI similarity.
● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ●●●● ● ●●● ●●●●●●●●● ●●●●●●●●●●● ●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●
2e+04 4e+04 6e+04 8e+04 1e+05
1.55 1.60 1.65 1.70 1.75 samples standard error
(f) Recip.⇥ BMI similarity.
Figure 5.4: The evolution of the condition number of the Jacobian, figure (a), and standard errors, figures (b) to (f), on N = 20,000 to 100,000 Monte Carlo simulations. The blue line indicates the N = 100,000 standard error estimate.
of standard error estimates based on N simulations. We do so 1000 times for various values of N .
Figure 5.5 shows the distributions of 1000 bootstrap samples of the female alter, reciprocity and BMI similarity standard errors estimated based on N = 1000, 5000 and 10,000 simulations. The figures show that the distributions of the standard errors have a fat tail. Increasing the number of simulations N shrinks this tail, but does not make it disappear. The three columns in Figure 5.5 represent three di↵erent situations:
1. For the female alter e↵ect, half of the estimates in the bootstrap distribu-tion for N = 1000 render the e↵ect significant and the other half does not. Given the distribution for N = 10,000, we would still want to increase the number of simulations N to get a more precise standard error estimate. 2. For the reciprocity e↵ect, 926 of the 1000 simulated standard errors
es-timated based on N = 5000 simulations would result in a significant reciprocity parameter. This strongly suggests that reciprocity plays a significant role in friendship evolution.
3. For the BMI similarity e↵ect, the standard error estimates are generally large compared to the estimated BMI similarity parameter of 0.66. This would not change by increasing the number of simulations N further.
Using the distribution of standard errors, we can pinpoint the source of the deviations in the standard errors.
log10( standard error )
frequency −1 0 1 2 3 0 100 200 300 400
(a) Female alter (452).
log10( standard error )
frequency −1 0 1 2 3 0 100 200 300 400 500 (b) Reciprocity (887).
log10( standard error )
frequency 0 1 2 3 0 100 200 300 400 (c) BMI similarity (0).
log10( standard error )
frequency −1.0 −0.5 0.0 0.5 1.0 0 100 200 300 400 500 600 700 (d) Female alter (622).
log10( standard error )
frequency −0.5 0.0 0.5 1.0 1.5 0 100 200 300 400 500 600 (e) Reciprocity (926).
log10( standard error )
frequency 0.0 0.5 1.0 1.5 2.0 0 100 200 300 400 (f) BMI similarity (0).
log10( standard error )
frequency −1.0 −0.5 0.0 0.5 1.0 1.5 0 200 400 600 800 (g) Female alter (700).
log10( standard error )
frequency −0.5 0.0 0.5 1.0 1.5 0 100 200 300 400 500 600 700 (h) Reciprocity (964).
log10( standard error )
frequency 0.0 0.5 1.0 1.5 2.0 0 100 200 300 400 500
(i) BMI similarity (0).
Figure 5.5: The bootstrap distributions of standard error estimates based on N = 1000 (first row), N = 5000 (second row) and N = 10,000 (third row) Monte Carlo simulations. The red line indicates the standard error value corresponding to a p-value of exactly 0.05. The numbers between brackets indicate the number of standard errors in the distribution rendering an e↵ect significant. Note the scale di↵erences on the horizontal axis.
● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 2.8 3.0 3.2 3.4 3.6 3.8 4.0 0 2 4 6
10log( deviation in cov(S) )
10log ( de viation in co v( θ) ) (a) Covariance ˆ⌃. ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 3.4 3.6 3.8 4.0 4.2 4.4 0 2 4 6
10log( deviation in Jacobian )
10log ( de viation in co v( θ) ) (b) Jacobian ˆD. ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ●●● ● ● ●● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● −1 0 1 2 3 0 2 4 6
10log( deviation in inverse Jacobian )
10log ( de viation in co v( θ) ) (c) Inverse Jacobian ˆD 1. ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●●● ● ● ● ● ● ● ●● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● 5.0 5.5 6.0 6.5 7.0 7.5 0 2 4 6
10log( condition number Jacobian )
10log ( de
viation in co
v(
θ) )
(d) Jacobian condition number.
Figure 5.6: Explaining deviations in cov(ˆ✓) by deviations in its constituents, the covariance of the statistics ˆ⌃ and the Jacobian ˆD.
For the 1000 bootstrap samples based on N = 1000, we compute the estimators of the covariance matrix of the statistics ˆ⌃, the Jacobian ˆD, the inverse of the Jacobian ˆD 1and the covariance of the parameter estimates cov(ˆ✓). Using the spectral norm, we determine the distances of these matrices to their counter-parts based on N = 100,000 simulations. Figures 5.6a and 5.6b show that the deviations in cov(ˆ✓) cannot be attributed to deviations in covariance ˆ⌃ or Ja-cobian ˆD. Large deviations in the inverse of the Jacobian, however, are a good predictor of deviations in cov(ˆ✓) (Figure 5.6c). The errors in standard errors are thus attributable to the deviations in the inverse of the Jacobian.
Finally, Figure 5.6d supports the theoretical findings in Section 5.3.3. Larger values of the condition number of the Jacobian correspond to larger deviations in
the parameter covariance matrix. The condition number of the Jacobian based on N = 100,000 simulations is 92,573. This is a large value, implying that small errors in ˆD may result in large errors in ˆD 1. This explains the variance in the deviations in cov(ˆ✓), visible in Figure 5.6d, for bootstrap samples with a Jacobian condition number value of around 105.
5.4.3 Exploring the dependencies
When the moment statistics are similarly sensitive to changes in the parameters ✓, the condition number of the Jacobian D is large. To explore which statistics are related, we can consider their covariance matrix ⌃.
Figure 5.7 illustrates the dependencies among the statistics, based on the co-variance matrix of the statistics computed on 100,000 Monte Carlo simulations. The graphs are based on partial correlations. In Figure 5.7b, we see that the partial correlation graph is composed of a component of BMI dynamics e↵ects, a component of BMI-related e↵ects on friendship evolution and a component containing the other e↵ects in the friendship dynamics model. If we increase the partial correlation threshold, we see that strongest dependence is among the purely structural e↵ects in the friendship dynamics model. This suggests that the problem we are studying is very general. Even in a network evolution model without co-evolving behavior or covariates, the strong dependence between the purely structural e↵ects could lead to inflated standard errors.
The dependencies among the statistics play out as expected, as is shown in Figure 5.7a. For example, the cycle “density – same home group – same home group⇥ reciprocity – reciprocity – density” constitutes an intuitive structure in a conditional independence graph. Density, which has the number of network ties as corresponding statistic, is central in the graphs.
Note that we did not choose to depict all partial correlations that significantly di↵er from 0, but only those larger than a certain threshold. A graph showing all significant partial correlations would have been much denser. In general, considering partial correlation graphs for various thresholds will give insight into the structure and strength of dependencies among statistics.
5.5 The e↵ect of a particular Monte Carlo simulation
In the previous section, we have seen that standard errors are sometimes esti-mated incorrectly. In this section, we will investigate to what extent deviations in standard error estimates are due to particular Monte Carlo simulations Yi.
Rate density recip gwespFF gwespBB gwespFFrecip inPop outAct altX egoX sameX sameX sameXRecip altX egoXsimXsimRecipX wiener feedback interceptavAlt (a)|⇢| > 0.2. Rate density recip gwespFF gwespBB gwespFFrecip inPop outAct altX egoXsameX sameX sameXRecip altX egoXsimXsimRecipX wiener feedback interceptavAlt (b)|⇢| > 0.1. Rate density recip gwespFF gwespBB gwespFFrecip inPop outAct altX egoXsameX sameX sameXRecip altX egoXsimXsimRecipX wiener feedback interceptavAlt (c)|⇢| > 0.3.
Figure 5.7: Partial correlation graphs based on N = 100,000 Monte Carlo simulations for 3 cut-o↵s of the absolute values of the partial correlations ⇢. The color of the nodes corresponds to the type of e↵ect (network dynamics: gray = purely structural, green = gender, yellow = home group, light blue = BMI; BMI dynamics: dark blue).
We investigate this in the context of the empirical example presented in Section 5.4. If ‘detrimental’ simulations exist, and if they are easily detectable, remov-ing them would o↵er a simple solution to the inflated standard error problem. In Section 5.5.1, we define a measure of the bias a particular simulation Yi in-duces in standard error estimates. Section 5.5.2 shows that none of the obvious statistics on Yiis good at detecting the Yi resulting in a large bias.
5.5.1 Defining the ‘detrimental’ simulations
One way to assess the e↵ect of a particular Monte Carlo simulation on standard error estimates based on N simulations is to re-estimate the standard errors with one simulation Yi left out. A large deviation of the newly computed estimates from the N = 100,000 estimate would indicate that the sample Yiis ‘bad’. This approach, however, would only indicate how Yi a↵ects the standard errors in a particular draw of N simulations. In the following, we instead assess the average e↵ect of samples Yi on deviations in standard error estimates.
We assess the average e↵ect of the first 50,000 Monte Carlo simulations on standard error estimates based on samples of size N = 5000. To this end, we take 1000 bootstrap samples of N = 5000 Monte Carlo simulations out of the 50,000 simulations. Based on each of the bootstrap samples, we compute standard error estimates, and the sum of the relative biases in the 21 estimates compared to the N = 100,000 standard error estimates,
21 X i=1 |s.e.5000(ˆ✓i) s.e.100,000(ˆ✓i)| s.e.100,000(ˆ✓i) . (5.25)
For each of the 50,000 Monte Carlo samples Yi, we count how often they occur in each of the bootstrap samples. The weighted average of the bias (5.25) in the bootstrap samples gives, for each of the 50,000 samples Yi, a measure of its average deviation in the standard error estimates.
Figure 5.8 shows that the distribution of the average bias among the 50,000 samples Yihas a fat tail. The average bias of 2044 of the 50,000 samples greatly exceeds 500. These samples are not included in the histogram in Figure 5.8a, but are visible in Figure 5.8b. The fat tail in the average bias distribution means that some samples are more often associated with inaccurate standard errors than others.
5.5.2 Detecting the ‘detrimental’ simulations?
The standard errors are a function of the statistics Si and scores Ji, evaluated on Monte Carlo samples Yi. Some function of the values Si and Ji therefore is a natural choice as indicator that a sample Yiis associated with a high average bias. In particular, we consider being an outlier in the multivariate distribution of the vectors (Si, Ji) as candidate indicator.
The existence of outliers in multivariate data can be checked using the Maha-lanobis distance. The sample MahaMaha-lanobis distance of a multivariate normally
average bias frequency 0 100 200 300 400 500 0 5000 10000 15000
(a) Original scale (excluding outliers).
10log( average bias )
frequency 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 0 2000 4000 6000 8000 10000 (b) Logarithmic scale.
Figure 5.8: Average bias of 50,000 Monte Carlo simulations in standard error estimates based on N = 5000 simulations.
distributed variable xifrom a sample x1, . . . , xN with mean ˆµ and covariance ˆ⌃ is given by MDi= ⇣ (xi µ)ˆ >⌃ˆ 1(xi µ)ˆ ⌘1/2 . (5.26)
For p-variate normally distributed data, the Mahalanobis distance squared is approximately chi-squared distributed with p degrees of freedom. Multivariate outliers can be observed as having a large Mahalanobis distance.
Figure 5.9a shows the density of the observed Mahalanobis distances of the vec-tors containing the statistics Si and scores Ji, and the chi-squared distribution with 42 degrees of freedom. The QQ-plot in Figure 5.9b indicates that there are some outliers in this distribution. Of the 50,000 values, 48 are larger than 80. This is more than the 18 values that would have been expected according to the theoretical distribution.
Nevertheless, Figure 5.10 shows that the average bias induced by a sample Yi is not related to the Mahalobis distance of the vector of statistics and scores (Si, Ji). Outliers in the distribution of statistics and scores do not predict a high bias in the standard errors. Considering the distributions of the statistics and scores separately yields the same conclusion.
In conclusion, our investigation into the source of a large average bias has not led to a usable diagnostic or a remedy against inflated standard errors. Other attempts, based on residuals after removing single simulations Yifrom the com-putation, have not led to usable diagnostics either. The same was true when we considered the statistics Si separately. We also tried considering the Maha-lanobis distances of the vectors containing the 21⇥ 21 = 441 entries of SiJi, but this was not possible as the variance of these vectors was not an invertible matrix.
Mahalanobis distance (MD) density 20 40 60 80 100 0.00 0.01 0.02 0.03 0.04
(a) Density plot. (b) Q-Q plot.
Figure 5.9: Distribution of the Mahalanobis distance of the statistics and scores in the N = 50,000 Monte Carlo simulations. The right figure contains a random draw of 10,000 of the 50,000 values.
Figure 5.10: Average bias by Mahalanobis distance for a random draw of 10,000 of the 50,000 values.
5.6 Alternative estimators
Because we so far have not found a good statistic indicating that a sample Yiis associated with a high average bias, in this section we consider two alternative standard error estimators. The first is motivated by the fact that deviations in the estimate of cov(✓) are caused by deviations in Jacobian ˆD, and employs a robust estimator of D. Instead of taking the mean of the summands in (5.10), this estimator uses their multivariate L1-median. We use a multivariate median here instead of an entry-wise median, because the entries of the matrix sum-mands are closely related. The multivariate L1-median of matrices A1, . . . , AN
is defined as
argminB2{A1,...,AN} N X
i=1
kvec(Ai) vec(B)k, (5.27)
wherek · k denotes the Euclidean norm. We approximate the L1-median using non-linear minimization, the best approach according to a comparison study by Fritz, Filzmoser, and Croux (2012).
The second estimator is motivated by the observation that outliers occur in the bootstrap distributions of the standard errors (see Figure 5.5). We can avoid these outliers by considering the medians of these distributions as our estimates. Therefore, the second alternative estimator is the element-wise median of b bootstrap samples of the standard errors based on N out of k given Monte Carlo simulations Yi.
Figure 5.11 compares the performance of the estimators. In this figure, the red line indicates the standard error value corresponding to a p-value of exactly 0.05. The blue line indicates the estimate based on N = 100,000 samples. The numbers between brackets are the number of standard errors in the distribution rendering an e↵ect significant.
The standard error estimates based on the L1-median Jacobian do not improve upon the regular score function estimates. The estimates based on the bootstrap median do. As expected, their range is considerably smaller than the range of the regular score function estimates, and on average they are closer to the standard error estimates based on the N = 100,000 samples.
If we increase the number of simulations k from which the bootstrap samples of size N are drawn, the range of the standard error values becomes smaller (Fig-ures 5.11 j,k,l). Notably, this range no longer includes the N = 100,000 samples estimate for the female alter e↵ect and barely includes it for the reciprocity and BMI similarity e↵ects. Considering bootstrap samples of size N = 5000 out of k = 5000 simulations, changes Figures 5.11 j,k,l only marginally (changed figures not reported here). Naturally, if we would increase N in the boot-strap estimator, the bootboot-strap median distribution would move towards the N = 100,000 estimate.
In conclusion, the alternative estimator based on the L1-median Jacobian does not outperform the regular score function estimator. Other alternative estima-tors based on a robust estimator for Jacobian D, using the entry-wise median or entry-wise trimmed mean, performed even worse. The bootstrap median estimator is more promising, as it does away with the large outliers, but yields biased results for small values of N . The shrinking range of the bootstrap es-timate for larger N and the uncertainty about for which value of N the bias is
log10( standard error ) frequency −1 0 1 2 3 0 100 200 300 400
(a) [SF] Female alter (452).
log10( standard error )
frequency −1 0 1 2 3 0 100 200 300 400 500 (b) [SF] Reciprocity (887).
log10( standard error )
frequency 0 1 2 3 0 100 200 300 400 (c) [SF] BMI similarity (0).
log10( standard error )
frequency −1.0 −0.5 0.0 0.5 1.0 1.5 2.0 0 100 200 300 400 500 (d) [L1] Female alter (307).
log10( standard error )
frequency −0.5 0.0 0.5 1.0 1.5 2.0 0 50 100 150 200 250 300 (e) [L1] Reciprocity (883).
log10( standard error )
frequency 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0 100 200 300 400 (f) [L1] BMI similarity (0).
log10( standard error )
frequency −0.8 −0.6 −0.4 −0.2 0 50 100 150 200 (g) [BS] Female alter (148).
log10( standard error )
frequency −0.4 −0.2 0.0 0.2 0 50 100 150 200 (h) [BS] Reciprocity (1000).
log10( standard error )
frequency 0.0 0.2 0.4 0.6 0 50 100 150 200
(i) [BS] BMI similarity (0).
log10( standard error )
frequency −0.7 −0.6 −0.5 −0.4 −0.3 0 50 100 150 200 250 300 (j) [BS] Female alter (209).
log10( standard error )
frequency −0.5 −0.4 −0.3 −0.2 −0.1 0.0 0 50 100 150 200 250 300 (k) [BS] Reciprocity (1000).
log10( standard error )
frequency 0.0 0.1 0.2 0.3 0 50 100 150 200 250 300 (l) [BS] BMI similarity (0).
Figure 5.11: The bootstrap distributions of standard error estimates based on the score function method with N = 1000 (SF, first row), the L1median method with N = 1000 (L1, second row), the median of 50 bootstrap samples of size N = 1000 from 1000 simulations (BS, third row), the median of 50 bootstrap samples of size N = 1000 from 5000 simulations (BS, fourth row).
