. * nwergm tutorial 3: MCMC diagnostics and goodness of fit . * Requires the core nwcommands and the package nwergm. . . nwwebuse lazega, nwclear . encode practice, generate(prac) . . * a fixed-size sample with a short interval between draws (the fits of nwergm up to 1.2): a chain that mixes badly . nwergm lazega_fr, edges mutual nodematch(prac) nodeicov(seniority) gwesp(0.5) seed(1) mcmcinterval(50) mcmcburnin(3000) > mcmleeffectivesize(0) Exponential-family random graph model Network: = lazega_fr Nodes: = 71 Ties: = 575 Directed: = Yes Estimation: = MCMLE MCMC sample size: = 3000 ------------------------------------------------------------------------------------ lazega_fr | Coefficient Std. err. z P>|z| [95% conf. interval] -------------------+---------------------------------------------------------------- edges | -3.322119 .1437037 -23.12 0.000 -3.603773 -3.040465 mutual | .810562 .1634309 4.96 0.000 .4902434 1.130881 nodematch_prac | .6795895 .090628 7.50 0.000 .5019619 .8572172 nodeicov_seniority | -.0043864 .0020888 -2.10 0.036 -.0084803 -.0002924 gwesp_0.5 | .5023746 .0653666 7.69 0.000 .3742584 .6304907 ------------------------------------------------------------------------------------ MCMLE converged after 1 iteration(s). Evaluating log-likelihood at the estimate (bridge sampling with 16 bridges) ... Null Deviance: 6889.9 on 4970 degrees of freedom Residual Deviance: 3339.2 on 4965 degrees of freedom Log-likelihood: -1669.6228 (df = 5) AIC: 3349.2 BIC: 3381.8 (MC std. err. = .31714) . estat mcmcdiag MCMC diagnostics for the final simulation (nwergm, method(mcmle)) Sample size: = 3000 Burn-in: = 3000 Thinning interval: = 50 Acceptance rate: = 0.6694 ------------------------------------------------------------ Statistic | Mean SD Autocorr ESS ----------------+------------------------------------------- edges | 578.0303 45.5309 0.9921 13.7 mutual | 74.7633 11.5793 0.9850 18.8 nodematch_prac | 407.8230 34.4362 0.9920 13.7 nodeicov_seniority|19478.4000 1781.8600 0.9913 15.1 gwesp_0.5 | 561.3006 86.7234 0.9913 14.6 ------------------------------------------------------------ Geweke diagnostic (first 10% vs. last 50% of the retained sample): --------------------------------------- Statistic | z p-value ----------------+---------------------- edges | 4.066 0.0000 * mutual | 3.279 0.0010 * nodematch_prac | 3.428 0.0006 * nodeicov_seniority| 5.118 0.0000 * gwesp_0.5 | 3.957 0.0001 * --------------------------------------- Note: |z| > 1.96 (flagged *) rejects convergence of that parameter's own chain at the 5% level (Geweke 1992) - a large sy > stematic difference between the retained sample's own early and late portions, not proof either way for a borderline z. Heidelberger-Welch stationarity + halfwidth test (pvalue=.05, eps=.1): -------------------------------------------------------------- Statistic | Stationary Start Halfwidth HW test ----------------+--------------------------------------------- edges | passed 1 24.1402 passed mutual | passed 1 5.2340 passed nodematch_prac | passed 1 18.2623 passed nodeicov_seniority| passed 301 658.1714 passed gwesp_0.5 | passed 1 44.4853 passed -------------------------------------------------------------- Note: Stationary = passed a Cramer-von-Mises test on the chain's own cumulative-sum path (Heidelberger & Welch 1983), dis > carding an increasing initial fraction (up to 50%) until it does - Start is the first retained iteration once it passes > , FAILED means no discard fraction achieved stationarity. HW test additionally requires the retained portion's own 95% > CI halfwidth to be within eps() (default 10%) of its own mean - a separate, stricter precision check, only meaningful o > nce stationarity itself has passed. MCMLE's own convergence test was satisfied after 1 iteration(s). Note: a satisfied convergence test is a necessary check, not proof of a converged fit - inspect the autocorrelation/ESS a > bove, and consider raising mcmleeffectivesize() (the target effective sample size of the adaptive sample; mcmleiteratio > ns() if the iterations ran out) if acceptance rate or ESS looks poor. For a fit with a fixed-size sample (mcmleeffectiv > esize(0)) raise mcmcsamplesize() or mcmcinterval() instead. Sample statistics summary (centered at the observed statistics): 1. Empirical mean and standard deviation for each variable, plus standard error of the mean: Mean SD Naive SE Time-series SE edges 3.0303 45.531 .83128 12.321 mutual .76333 11.579 .21141 2.6727 nodematch_prac -4.177 34.436 .62872 9.3206 nodeicov_senio~y 99.402 1781.9 32.532 458.57 gwesp_0.5 5.6815 86.723 1.5833 22.704 2. Quantiles for each variable: 2.5% 25% 50% 75% 97.5% edges -90 -23.25 4 30 101 mutual -22 -7 0 9 24 nodematch_prac -76 -24 -2 16 68 nodeicov_senio~y -3245.3 -1050.5 75 1189 4136 gwesp_0.5 -172 -43.098 6.4661 57.39 190.71 Are sample statistics significantly different from observed? edges mutual nodematch_prac nodeicov_sen~y gwesp_0.5 (Omni) diff. 3.0303 .76333 -4.177 99.402 5.6815 NA test stat. .24596 .28561 -.44815 .21677 .25024 17.149 P-val. .80571 .77518 .65405 .82839 .8024 .017808 Sample statistics cross-correlations: edges mutual nodematch_prac nodeicov_sen~y gwesp_0.5 edges 1.0000000 0.8271945 0.9363273 0.9605527 0.9815277 mutual 0.8271945 1.0000000 0.8181614 0.7958244 0.8318259 nodematch_prac 0.9363273 0.8181614 1.0000000 0.9043872 0.9372578 nodeicov_senio~y 0.9605527 0.7958244 0.9043872 1.0000000 0.9443242 gwesp_0.5 0.9815277 0.8318259 0.9372578 0.9443242 1.0000000 Sample statistics auto-correlation: Chain 1 edges mutual nodematch_prac nodeicov_sen~y gwesp_0.5 Lag 0 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 Lag 50 0.9909336 0.9847219 0.9909380 0.9899812 0.9903170 Lag 100 0.9820540 0.9706894 0.9822828 0.9801249 0.9812110 Lag 150 0.9732507 0.9573892 0.9734945 0.9705890 0.9724443 Lag 200 0.9647549 0.9457720 0.9647575 0.9612061 0.9642691 Lag 250 0.9561144 0.9346162 0.9559702 0.9513154 0.9559531 Sample statistics burn-in diagnostic (Geweke): Chain 1 Fraction in 1st window = 0.1 Fraction in 2nd window = 0.5 z-score Individual P-value edges 3.5560406 0.0003765 mutual 3.2698161 0.0010762 nodematch_prac 3.1562684 0.0015980 nodeicov_senio~y 4.6906069 0.0000027 gwesp_0.5 3.6125794 0.0003032 Joint P-value (lower = worse): 0.0180385 Note: MCMC diagnostics shown here are from the last round of simulation, prior to computation of final parameter estimate > s. Because the final estimates are refinements of those used for this simulation run, these diagnostics may understate > model performance. To directly assess the performance of the final model on in-model statistics, please use estat gof. . . * the defaults of nwergm 1.3 (R's interval and burn-in, adaptive sample size) . nwergm lazega_fr, edges mutual nodematch(prac) nodeicov(seniority) gwesp(0.5) seed(1) Exponential-family random graph model Network: = lazega_fr Nodes: = 71 Ties: = 575 Directed: = Yes Estimation: = MCMLE MCMC sample size: = 486 ------------------------------------------------------------------------------------ lazega_fr | Coefficient Std. err. z P>|z| [95% conf. interval] -------------------+---------------------------------------------------------------- edges | -3.32184 .131386 -25.28 0.000 -3.579352 -3.064328 mutual | .7969492 .1553399 5.13 0.000 .4924887 1.10141 nodematch_prac | .7214706 .080024 9.02 0.000 .5646264 .8783149 nodeicov_seniority | -.0046162 .0018014 -2.56 0.010 -.0081468 -.0010856 gwesp_0.5 | .4935235 .061269 8.06 0.000 .3734384 .6136085 ------------------------------------------------------------------------------------ MCMLE converged after 2 iteration(s). Evaluating log-likelihood at the estimate (bridge sampling with 16 bridges) ... Null Deviance: 6889.9 on 4970 degrees of freedom Residual Deviance: 3334.9 on 4965 degrees of freedom Log-likelihood: -1667.4731 (df = 5) AIC: 3344.9 BIC: 3377.5 (MC std. err. = .15171) . estat mcmcdiag MCMC diagnostics for the final simulation (nwergm, method(mcmle)) Sample size: = 486 Burn-in: = 16384 Thinning interval: = 2048 Requested interval: = 1024 (raised by the adaptive sampler) Acceptance rate: = 0.6692 ------------------------------------------------------------ Statistic | Mean SD Autocorr ESS ----------------+------------------------------------------- edges | 571.1543 39.0767 0.6726 95.6 mutual | 73.6770 10.4820 0.5538 120.0 nodematch_prac | 406.3889 30.0663 0.6635 98.8 nodeicov_seniority|19282.2500 1453.6430 0.6317 73.9 gwesp_0.5 | 546.8287 74.5329 0.6491 103.9 ------------------------------------------------------------ Geweke diagnostic (first 10% vs. last 50% of the retained sample): --------------------------------------- Statistic | z p-value ----------------+---------------------- edges | -1.610 0.1073 mutual | -1.972 0.0486 * nodematch_prac | -1.128 0.2594 nodeicov_seniority| -1.297 0.1947 gwesp_0.5 | -1.501 0.1335 --------------------------------------- Note: |z| > 1.96 (flagged *) rejects convergence of that parameter's own chain at the 5% level (Geweke 1992) - a large sy > stematic difference between the retained sample's own early and late portions, not proof either way for a borderline z. Heidelberger-Welch stationarity + halfwidth test (pvalue=.05, eps=.1): -------------------------------------------------------------- Statistic | Stationary Start Halfwidth HW test ----------------+--------------------------------------------- edges | passed 1 7.8316 passed mutual | passed 1 1.8753 passed nodematch_prac | passed 1 5.9273 passed nodeicov_seniority| passed 1 331.4554 passed gwesp_0.5 | passed 1 14.3329 passed -------------------------------------------------------------- Note: Stationary = passed a Cramer-von-Mises test on the chain's own cumulative-sum path (Heidelberger & Welch 1983), dis > carding an increasing initial fraction (up to 50%) until it does - Start is the first retained iteration once it passes > , FAILED means no discard fraction achieved stationarity. HW test additionally requires the retained portion's own 95% > CI halfwidth to be within eps() (default 10%) of its own mean - a separate, stricter precision check, only meaningful o > nce stationarity itself has passed. MCMLE's own convergence test was satisfied after 2 iteration(s). Note: a satisfied convergence test is a necessary check, not proof of a converged fit - inspect the autocorrelation/ESS a > bove, and consider raising mcmleeffectivesize() (the target effective sample size of the adaptive sample; mcmleiteratio > ns() if the iterations ran out) if acceptance rate or ESS looks poor. For a fit with a fixed-size sample (mcmleeffectiv > esize(0)) raise mcmcsamplesize() or mcmcinterval() instead. Sample statistics summary (centered at the observed statistics): 1. Empirical mean and standard deviation for each variable, plus standard error of the mean: Mean SD Naive SE Time-series SE edges -3.8457 39.077 1.7726 4.004 mutual -.32305 10.482 .47547 .95976 nodematch_prac -5.6111 30.066 1.3638 3.0304 nodeicov_senio~y -96.751 1453.6 65.939 169.99 gwesp_0.5 -8.7904 74.533 3.3809 7.3278 2. Quantiles for each variable: 2.5% 25% 50% 75% 97.5% edges -74.75 -31 -4.5 23 73 mutual -19.875 -8 -1 7 22 nodematch_prac -63 -27 -6.5 15 51 nodeicov_senio~y -2849.9 -1102 -83.5 961.5 2683.9 gwesp_0.5 -144.69 -61.393 -9.6405 40.06 147.53 Are sample statistics significantly different from observed? edges mutual nodematch_prac nodeicov_sen~y gwesp_0.5 (Omni) diff. -3.8457 -.32305 -5.6111 -96.751 -8.7904 NA test stat. -.96047 -.33659 -1.8516 -.56917 -1.1996 20.115 P-val. .33682 .73643 .064078 .56924 .23029 .0019311 Sample statistics cross-correlations: edges mutual nodematch_prac nodeicov_sen~y gwesp_0.5 edges 1.0000000 0.7530858 0.8939733 0.9175342 0.9702475 mutual 0.7530858 1.0000000 0.7393374 0.6798716 0.7810688 nodematch_prac 0.8939733 0.7393374 1.0000000 0.8465985 0.8966733 nodeicov_senio~y 0.9175342 0.6798716 0.8465985 1.0000000 0.8750858 gwesp_0.5 0.9702475 0.7810688 0.8966733 0.8750858 1.0000000 Sample statistics auto-correlation: Chain 1 edges mutual nodematch_prac nodeicov_sen~y gwesp_0.5 Lag 0 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 Lag 2048 0.6716965 0.5535584 0.6625514 0.6307216 0.6483841 Lag 4096 0.4637820 0.3595955 0.4536637 0.4201381 0.4494919 Lag 6144 0.3403473 0.2263601 0.3195065 0.3090722 0.3114812 Lag 8192 0.2762225 0.1376349 0.1992015 0.2888740 0.2483208 Lag 10240 0.1956788 0.0805804 0.1292762 0.2237654 0.1643815 Sample statistics burn-in diagnostic (Geweke): Chain 1 Fraction in 1st window = 0.1 Fraction in 2nd window = 0.5 z-score Individual P-value edges -1.2303640 0.2185608 mutual -1.9098929 0.0561470 nodematch_prac -1.1016837 0.2705992 nodeicov_senio~y -1.1809769 0.2376119 gwesp_0.5 -1.4362798 0.1509227 Joint P-value (lower = worse): 0.7788128 Note: MCMC diagnostics shown here are from the last round of simulation, prior to computation of final parameter estimate > s. Because the final estimates are refinements of those used for this simulation run, these diagnostics may understate > model performance. To directly assess the performance of the final model on in-model statistics, please use estat gof. . . * goodness of fit . estat gof, nsim(50) seed(2) Goodness-of-fit for model statistics obs min mean max MC p-value edges 575 503 586.20 685 0.80 mutual 74 58 74.20 96 1.00 nodematch_prac 412 366 421.06 493 0.84 nodeicov_senio~y 19379 15657 19581.72 22829 0.80 gwesp_0.5 555.61913 431.86804 580.31 786.826 0.76 Goodness-of-fit for minimum geodesic distance obs min mean max MC p-value 1 575 503 586.20 685 0.80 2 2316 2086 2553.72 2968 0.16 3 1782 1201 1684.34 1968 0.56 4 154 18 115.78 335 0.48 5 3 0 1.98 10 0.44 Inf 140 0 27.98 140 0.04 Goodness-of-fit for out-degree obs min mean max MC p-value 0 2 0 0.24 1 0.00 1 1 0 0.88 3 1.00 2 2 0 1.70 5 1.00 3 7 0 3.26 7 0.08 4 5 0 4.18 9 0.76 5 5 1 5.76 13 1.00 6 6 2 6.90 13 0.96 7 6 1 7.24 14 0.92 8 9 4 8.22 13 0.84 9 3 2 7.36 12 0.04 10 7 3 7.24 12 1.00 11 4 2 5.72 10 0.68 12 3 1 4.60 12 0.68 13 3 0 2.82 7 1.00 14 1 0 1.88 7 0.84 15 1 0 1.08 3 1.00 16 2 0 0.86 3 0.48 17 1 0 0.40 3 0.56 18 1 0 0.42 2 0.76 19 0 0 0.14 1 1.00 20 1 0 0.06 1 0.12 21 0 0 0.02 1 1.00 22 1 0 0.02 1 0.04 Goodness-of-fit for in-degree obs min mean max MC p-value 0 0 0 0.16 1 1.00 1 0 0 0.84 4 0.80 2 3 0 1.86 7 0.52 3 2 0 3.24 9 0.64 4 8 0 4.48 9 0.12 5 6 1 5.72 12 1.00 6 5 3 6.76 13 0.68 7 8 3 8.54 15 1.00 8 9 2 8.16 14 1.00 9 5 2 6.88 12 0.48 10 7 3 5.96 12 0.72 11 6 1 5.22 10 1.00 12 4 0 4.64 9 0.88 13 2 0 2.86 9 0.80 14 4 0 2.18 5 0.36 15 1 0 1.34 5 1.00 16 0 0 1.02 3 0.60 17 1 0 0.44 3 0.72 18 0 0 0.26 1 1.00 19 0 0 0.20 1 1.00 20 0 0 0.08 1 1.00 21 0 0 0.08 1 1.00 22 0 0 0.04 1 1.00 23 0 0 0.04 1 1.00 Goodness-of-fit for edgewise shared partner obs min mean max MC p-value 0 138 87 125.60 160 0.56 1 190 163 204.92 241 0.24 2 131 110 148.82 199 0.48 3 69 42 69.54 114 0.88 4 26 10 26.56 53 1.00 5 11 0 7.62 23 0.44 6 7 0 2.42 11 0.08 7 2 0 0.56 3 0.28 8 1 0 0.14 2 0.20 9 0 0 0.02 1 1.00 . . * a model that cannot be fitted: triangles in the Kapferer tailor shop network . nwwebuse kapferer, nwclear . capture noisily nwergm kapferer_t2, edges triangle seed(1) mcmcsamplesize(1000) Warning: MCMLE iter 2 - statistic(s) at column(s) 1 2 not varying this iteration - this may indicate an extreme point or > dead-end configuration. Unconstrained MCMC sampling did not mix at all. Optimization cannot continue. . display "return code: " _rc return code: 498 . . * the geometrically weighted version . nwergm kapferer_t2, edges gwesp(0.5) seed(1) Exponential-family random graph model Network: = kapferer_t2 Nodes: = 43 Ties: = 190 Directed: = No Estimation: = MCMLE MCMC sample size: = 463 ------------------------------------------------------------------------------ kapferer_t2 | Coefficient Std. err. z P>|z| [95% conf. interval] -------------+---------------------------------------------------------------- edges | -4.270901 .4311131 -9.91 0.000 -5.115867 -3.425935 gwesp_0.5 | 1.505898 .2434484 6.19 0.000 1.028747 1.983048 ------------------------------------------------------------------------------ MCMLE converged after 4 iteration(s). Evaluating log-likelihood at the estimate (bridge sampling with 16 bridges) ... Null Deviance: 1251.8 on 903 degrees of freedom Residual Deviance: 857.4 on 901 degrees of freedom Log-likelihood: -428.6868 (df = 2) AIC: 861.4 BIC: 871.0 (MC std. err. = .22782) . estat gof, nsim(50) seed(2) Goodness-of-fit for model statistics obs min mean max MC p-value edges 190 146 186.54 227 0.88 gwesp_0.5 275.77262 195.37931 269.97 342.59428 0.84 Goodness-of-fit for minimum geodesic distance obs min mean max MC p-value 1 190 146 186.54 227 0.88 2 535 397 526.12 614 0.84 3 163 62 139.84 233 0.60 4 15 0 7.38 48 0.28 5 0 0 0.22 6 1.00 Inf 0 0 42.90 200 0.88 Goodness-of-fit for degree obs min mean max MC p-value 0 0 0 1.00 5 0.88 1 2 0 0.54 3 0.24 2 3 0 1.28 6 0.32 3 2 0 1.70 6 1.00 4 0 0 2.36 6 0.20 5 6 0 3.08 6 0.08 6 1 0 3.34 6 0.20 7 4 1 3.54 8 0.92 8 2 2 4.08 8 0.32 9 3 0 4.26 9 0.60 10 8 0 3.78 9 0.08 11 4 1 3.42 7 0.84 12 1 0 3.14 8 0.32 13 0 0 2.22 8 0.20 14 0 0 1.80 6 0.32 15 3 0 1.16 4 0.20 16 1 0 0.96 3 1.00 17 0 0 0.60 3 1.00 18 1 0 0.30 2 0.48 19 1 0 0.22 2 0.36 20 0 0 0.16 2 1.00 21 0 0 0.02 1 1.00 22 0 0 0.02 1 1.00 24 0 0 0.02 1 1.00 25 1 0 0.00 0 0.00 Goodness-of-fit for edgewise shared partner obs min mean max MC p-value 0 6 0 2.56 7 0.08 1 23 7 22.78 36 1.00 2 33 31 46.90 66 0.08 3 24 30 49.12 69 0.00 4 32 21 33.60 54 0.96 5 27 0 18.52 38 0.28 6 19 1 8.24 21 0.04 7 13 0 3.22 12 0.00 8 5 0 1.04 9 0.04 9 2 0 0.38 3 0.16 10 2 0 0.16 2 0.04 11 0 0 0.02 1 1.00 12 3 0 0.00 0 0.00 13 1 0 0.00 0 0.00 .