--- title: "Getting Started with 'stansum'" output: rmarkdown::html_vignette bibliography: refs.bib vignette: > %\VignetteIndexEntry{Getting Started with 'stansum'} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} editor_options: chunk_output_type: console --- ``` r library(stansum) ``` # Introduction Aggregated Relational Data (ARD) are answers to survey questions of the form "How many people do you know who are *X*?", where *X* is a subpopulation such as "named Michael", "are medical doctors", or "died in a car accident last year". A respondent $i = 1, \dots, N$ answers such a question for each of $K$ subpopulations, and the resulting $N \times K$ table of counts $y_{ik}$ is the input of every model in this package. The main quantities of interest are the sizes of the personal networks of the respondents (their *degrees*), and -- in the Network Scale-Up Method (NSUM) -- the sizes of subpopulations that are hard to count directly. If the sizes of some of the subpopulations are known, they anchor the scale of the degrees; a hidden subpopulation can then be sized by comparing how many of its members respondents report knowing relative to the known ones. The **stansum** package wraps a collection of precompiled Stan models for ARD proposed by @zheng_etal_2006, @maltiel_etal_2015 and @baum_marsden_2023. The models are embedded in the package using **instantiate** [@r-instantiate] and sampled with **cmdstanr** [@r-cmdstanr]. Two model families are covered: - the overdispersed Poisson models of @zheng_etal_2006, in which degrees and subpopulation prevalences enter through log-linear respondent and subpopulation effects; - the models of @maltiel_etal_2015, in which degrees are explicit parameters with a lognormal prior and the counts are binomial or beta-binomial. Most models come in a *count* version, which uses the answers as they are, and in *dichotomous* (knows anyone in the subpopulation or not) and *trichotomous* (nobody / one person / two or more) versions proposed by @baum_marsden_2023 for surveys that collect coarsened answers. This vignette fits every model once. The examples are deliberately as computationally cheap as possible: they use a small subset of the data and a tiny number of MCMC iterations. Their purpose is to show the calling syntax, the structure of the fitted objects, and give a flavor of how to interpret the results -- **not** to produce meaningful estimates. # Listing and loading models All models shipped with the package are listed by `list_models()`: ``` r list_models() #> [1] "bernoulli" "CovMod_count" #> [3] "CovMod_dichotomous" "MaltielBEM_count" #> [5] "MaltielBEM_dichotomous" "MaltielBEM_trichotomous" #> [7] "MaltielCM_count" "MaltielRDM_count" #> [9] "MaltielRDM_dichotomous" "MaltielTBM_count" #> [11] "ZhengBEM_count" "ZhengBEM_count2" #> [13] "ZhengBEM_dichotomous" "ZhengBEM_trichotomous" #> [15] "ZhengGP_count" "ZhengGP_dichotomous" ``` The names are the base names of the Stan files. Each model has an R wrapper named after it in lowercase, e.g. `zheng_bem_count()` for `ZhengBEM_count`. (The `CovMod_*` models and the `bernoulli` test model are also listed; the former do not have R wrappers yet.) Use `get_model()` to load a model as a **cmdstanr** object of R6 class "CmdStanModel" without fitting it: ``` r mod <- get_model("ZhengBEM_count") class(mod) #> [1] "CmdStanModel" "R6" basename(mod$exe_file()) #> [1] "ZhengBEM_count" ``` Note that the object is created from the precompiled executable rather than a Stan file, so methods that need the source (such as `$print()` or `$variables()`) are unavailable. The Stan sources are installed alongside the executables: ``` r stan_file <- system.file("bin", "stan", "ZhengBEM_count.stan", package = "stansum") cat(readLines(stan_file), sep = "\n") #> data { #> int N; #> int K; #> array[N, K] int y; #> } #> #> parameters { #> vector[N] alpha; #> vector[K] beta; #> vector[K] inv_omega; #> real mu_beta; #> real sigma_beta; #> real sigma_alpha; #> } #> #> model { #> alpha ~ normal(0, sigma_alpha); #> beta ~ normal(mu_beta, sigma_beta); #> for (k in 1:K) { #> real omega_k_m1; #> omega_k_m1 = inv(inv(inv_omega[k]) - 1); #> for (i in 1:N) { #> if (y[i,k] >= 0) { #> real xi_i_k; #> xi_i_k = omega_k_m1 * exp(alpha[i] + beta[k]); #> y[i,k] ~ neg_binomial(xi_i_k, omega_k_m1); #> } #> } #> } #> } ``` The source is the definitive description of a model: the `data` block lists the inputs the model expects (and the R wrapper has one argument for each of them, with the same name), the `parameters` block lists what will show up in the posterior summaries, and the `model` block spells out the priors and the likelihood. The remaining sections use the R wrappers which load the model and sample from the posterior in one call. # Example data All examples use the artificial dataset `Fake_maltiel_RD` shipped with the package. It was simulated from the Random Degree model of @maltiel_etal_2015 (see below) for 1000 respondents and three subpopulations of sizes 100000, 200000 and 300000 in a population of 3 million. The degrees were drawn from a lognormal distribution with a median of about 20. In a real survey the subpopulations would be chosen so that their sizes are known from a census or a register -- these known sizes are what pin down the scale of the degree estimates. See `?Fake` for details. ``` r data("Fake_maltiel_RD", package = "stansum") head(Fake_maltiel_RD) #> p1 p2 p3 #> 1 1 3 3 #> 2 1 1 5 #> 3 2 0 0 #> 4 2 2 1 #> 5 0 2 1 #> 6 0 3 1 ``` To keep the computations cheap we use the first 50 respondents only. All models take the ARD as an integer matrix `y` with respondents in rows and subpopulations in columns, together with its dimensions `N` (number of respondents) and `K` (number of subpopulations): ``` r y <- data.matrix(Fake_maltiel_RD)[1:50, ] N <- nrow(y) K <- ncol(y) ``` The count models skip cells with negative values, so a missing answer can be encoded as, say, `-1` rather than `NA` (which Stan does not accept). Some of the models accept dichotomous (does the respondent know anyone in the subpopulation: 0/1) or trichotomous (knows nobody / one person / two or more: 0/1/2) responses instead of counts [@baum_marsden_2023]. Such data arise when the questionnaire only asks "Do you know anyone who...?" or offers a small number of answer categories. We derive them from the counts: ``` r y01 <- 1L * (y > 0) y012 <- pmin(y, 2) ``` The Maltiel et al (2015) models additionally need - `m` -- a vector of length `K` with the sizes of the subpopulations as fractions of the population size, hence numbers in (0, 1); - `L` -- a vector of length `N` with a lower bound on the degree of each respondent. The count models parametrize the degree as `L + d_raw` with `d_raw > 0`, so `L[i]` must be at least the largest number of people respondent `i` reported knowing in any subpopulation. Using exactly that maximum is the natural choice: ``` r m <- c(1, 2, 3) * 1e5 / 3e6 L <- apply(y, 1, max) ``` # Sampling settings Every wrapper passes `...` to the `$sample()` method of the "CmdStanModel" object, see `?cmdstanr::``model-method-sample``. All calls below use - `chains = 1`, `iter_warmup = 100`, `iter_sampling = 100` -- one Markov chain, 100 iterations of warmup (adaptation of the sampler; these draws are discarded) and 100 iterations retained as posterior draws. These are the cheapest settings that still run in a reasonable manner. Real applications require multiple chains (e.g. the **cmdstanr** default of 4, which can be run concurrently with `parallel_chains`) and far more iterations (hundreds to thousands of each kind) so that convergence can be assessed and the posterior summaries are precise. - `refresh = 0`, `show_messages = FALSE` and `show_exceptions = FALSE` -- suppress progress output and informational messages from Stan. Leave them at their defaults when working interactively; the messages are useful. - `seed = 666` -- for reproducibility of the random draws. With so few iterations the posterior summaries come with warnings about unreliable effective sample sizes (ESS). They are expected here and are not shown in this vignette. The arguments are spelled out in every call so that each example is self-contained; the text refers to them collectively as "the sampling settings". # A test model `test_model()` wraps a simple Bernoulli model borrowed from the Stan documentation. It is useful for checking that the toolchain works. Its inputs are `N`, the number of binary observations, and `y`, the vector of those observations. The remaining arguments are the sampling settings: ``` r y_test <- c(1, 1, 1, 1, 0, 0, 0, 0, 0, 0) fit <- test_model( N = length(y_test), y = y_test, chains = 1, iter_warmup = 100, iter_sampling = 100, refresh = 0, show_messages = FALSE, show_exceptions = FALSE, seed = 666 ) fit #> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail #> lp__ -8.77 -8.48 0.77 0.35 -10.21 -8.18 1.03 39 54 #> theta 0.47 0.51 0.14 0.15 0.24 0.69 1.00 57 78 ``` Printing a fitted object shows a summary of the posterior. The only model parameter is `theta`, the probability of success, and its posterior mean is close to the sample proportion of 0.4 (the `lp__` row is the log posterior density up to a constant, useful for diagnostics but not a parameter). # Zheng et al (2006) Barrier Effects Model See `?zheng_bem`. In this model the count $y_{ik}$ is negative binomial with mean $\exp(\alpha_i + \beta_k)$, where - `alpha[i]` is a respondent effect, the log degree of respondent $i$ up to an additive constant (higher values = more gregarious respondents), - `beta[k]` is a subpopulation effect, the log prevalence of subpopulation $k$ up to the same constant, - `inv_omega[k]` in (0, 1] is the inverse of the overdispersion of subpopulation $k$: the variance of the counts is $\omega_k$ times the mean, so `inv_omega[k]` close to 1 means Poisson-like counts, while small values mean that some respondents know many more members of subpopulation $k$ than their degree alone would suggest -- a *barrier effect* (social structure that concentrates ties to a subpopulation among some people). The hyperparameters `mu_beta` and `sigma_beta` describe the distribution of the subpopulation effects and `sigma_alpha` the spread of the respondent effects. ## Count responses The inputs are `N` and `K`, the dimensions of the ARD, and `y`, the matrix of counts. The remaining arguments are the sampling settings: ``` r fit <- zheng_bem_count( N = N, K = K, y = y, chains = 1, iter_warmup = 100, iter_sampling = 100, refresh = 0, show_messages = FALSE, show_exceptions = FALSE, seed = 666 ) fit #> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail #> lp__ -163.26 -164.74 13.02 13.47 -181.32 -139.61 1.10 7 21 #> alpha[1] 0.07 0.06 0.18 0.17 -0.21 0.36 1.00 172 78 #> alpha[2] 0.08 0.06 0.22 0.24 -0.26 0.39 0.99 159 52 #> alpha[3] -0.07 -0.07 0.16 0.16 -0.36 0.17 1.03 137 112 #> alpha[4] -0.01 0.01 0.17 0.19 -0.32 0.22 0.99 200 112 #> alpha[5] -0.03 -0.01 0.17 0.15 -0.31 0.26 1.00 200 51 #> alpha[6] 0.00 0.00 0.15 0.14 -0.27 0.26 1.02 177 60 #> alpha[7] 0.02 0.02 0.17 0.13 -0.24 0.31 1.01 200 37 #> alpha[8] 0.00 0.00 0.16 0.13 -0.29 0.28 1.00 167 61 #> alpha[9] 0.01 0.02 0.15 0.13 -0.19 0.24 1.06 146 78 #> #> # showing 10 of 60 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option) ``` The printed summary is truncated to the first 10 of the 60 rows (the `max_rows` argument of the print method controls this), so with 50 respondents we only see some of the `alpha`s. The `$summary()` method described next gives access to any subset. ## Structure of the fitted object All wrappers return an object of R6 class "CmdStanMCMC" from **cmdstanr**. Let us have a closer look at this one; all other fits below have the same structure. ``` r class(fit) #> [1] "CmdStanMCMC" "CmdStanFit" "R6" ``` Basic information about the model and the run is available through `$metadata()`, e.g. the model name, the names of all variables in the posterior, and the number of sampling iterations: ``` r fit$metadata()$model_name #> [1] "ZhengBEM_count_model" fit$metadata()$stan_variables #> [1] "lp__" "alpha" "beta" "inv_omega" "mu_beta" #> [6] "sigma_beta" "sigma_alpha" fit$metadata()$iter_sampling #> [1] 100 ``` Posterior summaries for all or selected variables are computed by `$summary()`. Vector parameters such as `alpha` are expanded to one row per element, and asking for `"alpha"` returns all of them: ``` r fit$summary() #> # A tibble: 60 × 10 #> variable mean median sd mad q5 q95 rhat ess_bulk #> #> 1 lp__ -163. -1.65e+2 13.0 13.5 -181. -140. 1.10 7.91 #> 2 alpha[1] 0.0672 5.82e-2 0.177 0.173 -0.208 0.360 1.00 173. #> 3 alpha[2] 0.0774 6.38e-2 0.217 0.242 -0.265 0.393 0.991 159. #> 4 alpha[3] -0.0741 -6.52e-2 0.159 0.157 -0.355 0.167 1.03 137. #> 5 alpha[4] -0.0103 7.59e-3 0.173 0.189 -0.323 0.216 0.994 200 #> 6 alpha[5] -0.0288 -1.32e-2 0.166 0.148 -0.305 0.264 1.00 200 #> 7 alpha[6] -0.00453 -2.83e-5 0.145 0.137 -0.274 0.260 1.02 178. #> 8 alpha[7] 0.0162 2.24e-2 0.171 0.135 -0.244 0.313 1.01 200 #> 9 alpha[8] -0.00123 -2.27e-3 0.158 0.130 -0.288 0.280 1.00 167. #> 10 alpha[9] 0.0112 2.23e-2 0.145 0.135 -0.194 0.238 1.06 147. #> # ℹ 50 more rows #> # ℹ 1 more variable: ess_tail fit$summary(variables = c("mu_beta", "sigma_beta", "sigma_alpha")) #> # A tibble: 3 × 10 #> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail #> #> 1 mu_beta 1.55 0.395 5.21 1.31 -8.64 8.92 1.88 8.38 21.6 #> 2 sigma_beta 9.78 4.97 12.5 6.42 0.298 39.5 2.08 1.47 21.6 #> 3 sigma_alpha 0.173 0.173 0.0474 0.0520 0.0982 0.255 1.02 13.6 31.5 ``` Each row gives the posterior mean and median, the standard deviation and the median absolute deviation, and the 5th and 95th percentiles (a 90% credible interval) of a variable. The last three columns are convergence diagnostics: `rhat` compares the chains (and parts of chains) with each other and should be close to 1 (values above 1.01 are a warning sign), and `ess_bulk` and `ess_tail` are the effective numbers of independent draws behind the estimate of the center and the tails of the posterior, respectively. With 100 draws from one chain these numbers are of course tiny; in a real analysis one expects several hundred at least. The posterior draws themselves are returned by `$draws()` in one of the formats supported by the **posterior** package -- by default a `draws_array` with dimensions iterations x chains x variables: ``` r draws <- fit$draws() class(draws) #> [1] "draws_array" "draws" "array" dim(draws) #> [1] 100 1 60 head(fit$draws(format = "df")) #> # A draws_df: 6 iterations, 1 chains, and 60 variables #> lp__ alpha[1] alpha[2] alpha[3] alpha[4] alpha[5] alpha[6] alpha[7] #> 1 -156 -0.094 -0.097 0.144 0.160 -0.0412 -0.015 0.0058 #> 2 -165 0.272 0.002 -0.263 -0.032 0.0059 -0.061 -0.0147 #> 3 -174 0.060 0.357 -0.012 0.043 0.0583 -0.156 -0.0461 #> 4 -171 0.186 0.132 0.074 -0.377 0.2647 -0.061 0.0259 #> 5 -167 0.168 0.380 -0.099 0.080 -0.2247 0.082 -0.0683 #> 6 -158 -0.111 -0.061 -0.154 -0.184 -0.1116 0.086 -0.1760 #> # ... with 52 more variables #> # ... hidden reserved variables {'.chain', '.iteration', '.draw'} ``` The data frame format (`format = "df"`) is the most convenient one for computing derived quantities: every row is one draw of all parameters, so any function of the parameters can be evaluated draw by draw and summarised afterwards, which propagates the posterior uncertainty correctly. Sampler diagnostics -- the number of divergent transitions, the number of iterations that hit the maximum tree depth, and the energy-based diagnostic E-BFMI -- are returned by `$diagnostic_summary()`. Divergences in particular indicate that the sampler struggled and the draws may be biased: ``` r fit$diagnostic_summary() #> $num_divergent #> [1] 0 #> #> $num_max_treedepth #> [1] 0 #> #> $ebfmi #> [1] 0.3318877 ``` Note that the draws live in CSV files written by CmdStan (see `fit$output_files()`) that are by default created in a temporary directory and are deleted when the R session ends. Use `fit$save_object()` to save the whole fitted object. See `?cmdstanr::CmdStanMCMC` for the complete list of methods. ## Interpreting the fit Because the mean of the counts depends on `alpha[i] + beta[k]` only, the model cannot tell how much of that sum is due to the respondent and how much to the subpopulation: adding a constant to all `alpha`s and subtracting it from all `beta`s leaves the likelihood unchanged. The priors center the `alpha`s at zero, which fixes the scale arbitrarily. @zheng_etal_2006 pin it down using the known subpopulation sizes: the prevalences `exp(beta[k])` of the known subpopulations should add up to their known combined share of the population. In our data all three subpopulations are known, so the shift `C` is the log of the ratio of the estimated combined share `sum(exp(beta))` to the known one `sum(m)`; subtracting it from the `beta`s makes them log prevalences and adding it to the `alpha`s makes them log degrees, `exp(alpha + C)`: ``` r beta_hat <- fit$summary("beta")$mean alpha_hat <- fit$summary("alpha")$mean C <- log(sum(exp(beta_hat))) - log(sum(m)) degree <- exp(alpha_hat + C) summary(degree) #> Min. 1st Qu. Median Mean 3rd Qu. Max. #> 19.19 20.45 21.27 21.18 21.98 23.96 ``` This is a quick point estimate based on posterior means; a proper analysis applies the same computation to every draw in `fit$draws(format = "df")` to obtain full posterior distributions of the degrees. (The data were simulated with a median degree of about 20, so even this toy run gets the order of magnitude right.) The overdispersion parameters tell whether there are barrier effects with respect to any of the subpopulations: ``` r fit$summary("inv_omega") #> # A tibble: 3 × 10 #> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail #> #> 1 inv_omega[1] 0.791 0.798 0.100 0.108 0.607 0.916 1.10 25.8 33.2 #> 2 inv_omega[2] 0.876 0.912 0.0909 0.0735 0.706 0.966 1.25 3.74 48.4 #> 3 inv_omega[3] 0.857 0.875 0.0772 0.0793 0.726 0.952 1.05 31.0 111. ``` The data were simulated without barrier effects, so `inv_omega` should be close to 1 for all three subpopulations, but with 100 draws the estimates are noisy. ## Count responses (alternative parametrization) `zheng_bem_count2()` fits the same likelihood with a different parametrization and different priors. The respondent and subpopulation effects are sampled as standardized deviations (`alpha_std`, `beta_std`) and scaled afterwards, which typically makes the sampler more efficient for hierarchical models; the overdispersion is parametrized as `eta[k] = omega[k] - 1` with a gamma prior; and the hyperparameters get weakly informative normal priors. `alpha`, `beta`, `omega` and `inv_omega` are computed as transformed parameters and generated quantities, so the summaries are directly comparable with the previous fit. The inputs are the same: `N`, `K` and the count matrix `y`, followed by the sampling settings: ``` r fit <- zheng_bem_count2( N = N, K = K, y = y, chains = 1, iter_warmup = 100, iter_sampling = 100, refresh = 0, show_messages = FALSE, show_exceptions = FALSE, seed = 666 ) fit #> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail #> lp__ -260.09 -259.93 6.18 6.76 -270.16 -250.49 1.05 27 52 #> alpha_std[1] 0.16 0.30 1.00 1.15 -1.50 1.63 1.02 152 53 #> alpha_std[2] 0.21 0.23 0.84 0.79 -1.27 1.56 1.00 90 76 #> alpha_std[3] -0.23 -0.17 0.87 0.78 -1.90 1.16 0.99 133 54 #> alpha_std[4] 0.07 0.06 1.00 1.10 -1.64 1.67 0.99 101 78 #> alpha_std[5] -0.14 -0.12 0.98 0.98 -1.71 1.33 1.00 136 113 #> alpha_std[6] -0.06 -0.09 1.16 1.35 -1.73 1.93 1.03 200 32 #> alpha_std[7] 0.14 0.15 0.84 0.75 -1.23 1.46 1.00 126 92 #> alpha_std[8] -0.09 -0.09 0.82 0.84 -1.56 1.20 1.00 99 78 #> alpha_std[9] -0.05 -0.18 1.02 0.81 -1.71 2.14 1.00 102 54 #> #> # showing 10 of 119 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option) ``` ## Dichotomous responses For dichotomous data the model only uses the probabilities, implied by the same negative binomial distribution, of reporting nobody versus at least one person in each subpopulation. The parameters and their interpretation are as for the count model, but only the "knows nobody" versus "knows someone" contrast is informative, so the degree scale is estimated much less precisely [@baum_marsden_2023]. The inputs are `N`, `K` and the 0/1 matrix `y01`, followed by the sampling settings: ``` r fit <- zheng_bem_dichotomous( N = N, K = K, y = y01, chains = 1, iter_warmup = 100, iter_sampling = 100, refresh = 0, show_messages = FALSE, show_exceptions = FALSE, seed = 666 ) #> Warning: 1 of 1 chains had an E-BFMI less than 0.3. #> See https://mc-stan.org/misc/warnings for details. fit #> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail #> lp__ -62.87 -65.27 20.23 21.34 -90.63 -24.71 1.24 3 21 #> alpha[1] 0.28 0.21 0.52 0.43 -0.46 1.05 1.05 77 78 #> alpha[2] 0.27 0.21 0.47 0.33 -0.41 1.17 1.10 7 48 #> alpha[3] -0.44 -0.37 0.45 0.46 -1.14 0.19 1.09 7 68 #> alpha[4] 0.23 0.17 0.54 0.50 -0.51 1.07 1.01 105 65 #> alpha[5] -0.06 -0.04 0.49 0.33 -0.90 0.73 1.09 137 35 #> alpha[6] -0.02 -0.02 0.37 0.29 -0.58 0.54 1.02 195 30 #> alpha[7] -0.06 -0.05 0.47 0.46 -0.79 0.67 1.10 200 31 #> alpha[8] 0.02 0.01 0.36 0.32 -0.51 0.71 1.04 200 54 #> alpha[9] -0.04 -0.01 0.45 0.40 -0.87 0.72 1.20 200 73 #> #> # showing 10 of 60 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option) ``` ## Trichotomous responses The trichotomous model additionally distinguishes knowing exactly one person from knowing two or more, which recovers some of the information lost when dichotomizing. The inputs are `N`, `K` and the 0/1/2 matrix `y012`, followed by the sampling settings: ``` r fit <- zheng_bem_trichotomous( N = N, K = K, y = y012, chains = 1, iter_warmup = 100, iter_sampling = 100, refresh = 0, show_messages = FALSE, show_exceptions = FALSE, seed = 666 ) #> Warning: 1 of 1 chains had an E-BFMI less than 0.3. #> See https://mc-stan.org/misc/warnings for details. fit #> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail #> lp__ -114.42 -123.43 27.24 16.14 -143.97 -56.96 1.17 5 16 #> alpha[1] 0.21 0.17 0.34 0.30 -0.24 0.74 1.07 15 101 #> alpha[2] 0.02 0.05 0.24 0.18 -0.40 0.38 1.02 130 76 #> alpha[3] -0.26 -0.19 0.35 0.23 -0.96 0.10 1.05 24 21 #> alpha[4] 0.06 0.01 0.27 0.27 -0.36 0.52 1.01 179 78 #> alpha[5] -0.07 -0.04 0.29 0.22 -0.54 0.42 1.05 190 51 #> alpha[6] -0.08 -0.03 0.36 0.31 -0.69 0.43 1.15 95 22 #> alpha[7] -0.04 -0.03 0.29 0.23 -0.51 0.44 1.14 195 19 #> alpha[8] 0.05 0.02 0.37 0.24 -0.47 0.48 0.99 115 44 #> alpha[9] -0.06 -0.03 0.32 0.28 -0.65 0.36 1.04 189 71 #> #> # showing 10 of 60 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option) ``` # Zheng et al (2006) Group Prevalence Model See `?zheng_gp`. This is the Poisson special case of the Barrier Effects Model: the counts $y_{ik}$ are Poisson with mean $\exp(\alpha_i + \beta_k)$ and there are no overdispersion parameters. `alpha`, `beta`, `mu_beta`, `sigma_beta` and `sigma_alpha` have the same meaning as before and the rescaling of the previous section applies unchanged. Comparing this fit with the Barrier Effects Model is one way of judging whether overdispersion matters for a given dataset. ## Count responses The inputs are `N`, `K` and the count matrix `y`, followed by the sampling settings: ``` r fit <- zheng_gp_count( N = N, K = K, y = y, chains = 1, iter_warmup = 100, iter_sampling = 100, refresh = 0, show_messages = FALSE, show_exceptions = FALSE, seed = 666 ) #> Warning: 1 of 1 chains had an E-BFMI less than 0.3. #> See https://mc-stan.org/misc/warnings for details. fit #> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail #> lp__ -58.87 -59.24 16.67 16.62 -86.15 -32.96 1.24 3 21 #> alpha[1] 0.08 0.08 0.19 0.16 -0.28 0.43 1.01 119 78 #> alpha[2] 0.10 0.07 0.25 0.20 -0.16 0.60 1.01 165 18 #> alpha[3] -0.09 -0.09 0.17 0.17 -0.36 0.16 0.99 113 49 #> alpha[4] 0.01 0.01 0.22 0.20 -0.32 0.33 1.12 200 14 #> alpha[5] -0.02 -0.02 0.19 0.21 -0.34 0.27 1.06 191 45 #> alpha[6] -0.02 -0.02 0.17 0.18 -0.28 0.25 1.05 200 60 #> alpha[7] 0.01 0.02 0.16 0.14 -0.28 0.28 1.07 200 21 #> alpha[8] -0.01 -0.01 0.16 0.15 -0.22 0.26 0.99 200 30 #> alpha[9] -0.03 -0.01 0.19 0.23 -0.32 0.26 1.01 113 30 #> #> # showing 10 of 57 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option) ``` ## Dichotomous responses The dichotomous version uses the Poisson probability of a zero count versus a positive one. The inputs are `N`, `K` and the 0/1 matrix `y01`, followed by the sampling settings: ``` r fit <- zheng_gp_dichotomous( N = N, K = K, y = y01, chains = 1, iter_warmup = 100, iter_sampling = 100, refresh = 0, show_messages = FALSE, show_exceptions = FALSE, seed = 666 ) #> Warning: 1 of 1 chains had an E-BFMI less than 0.3. #> See https://mc-stan.org/misc/warnings for details. fit #> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail #> lp__ 183.56 120.32 122.41 25.08 98.32 415.32 1.81 1 13 #> alpha[1] 0.36 0.26 0.67 0.47 -0.65 1.67 1.06 75 46 #> alpha[2] 0.36 0.36 0.48 0.44 -0.42 1.20 0.99 74 73 #> alpha[3] -0.60 -0.50 0.59 0.49 -1.57 0.12 1.03 22 71 #> alpha[4] 0.31 0.28 0.49 0.41 -0.55 1.04 1.08 106 101 #> alpha[5] -0.01 0.01 0.60 0.55 -1.11 0.89 1.07 200 78 #> alpha[6] -0.03 0.01 0.51 0.37 -0.81 0.55 1.01 126 18 #> alpha[7] -0.03 -0.08 0.64 0.68 -1.16 1.02 0.99 200 113 #> alpha[8] -0.01 -0.07 0.44 0.44 -0.59 0.62 1.01 156 54 #> alpha[9] -0.01 -0.02 0.46 0.40 -0.75 0.72 0.99 200 49 #> #> # showing 10 of 57 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option) ``` # Maltiel et al (2015) Barrier Effects Model See `?maltiel_bem`. The models of @maltiel_etal_2015 treat the degrees as parameters directly: `d[i]` is the degree of respondent $i$, and the degrees follow a lognormal distribution with parameters `mu` and `sigma` (on the log scale, with uniform priors on (3, 8) and (1/4, 2), respectively). Given the degree, the count $y_{ik}$ is the number of "successes" among $d_i$ contacts, each of whom belongs to subpopulation $k$ with probability related to its known prevalence $m_k$. In the Barrier Effects Model that probability varies between respondents around $m_k$ according to a beta distribution, which makes the counts beta-binomial. The spread of that distribution is governed by `rho[k]` in (0, 1): `rho[k]` close to 0 means no barrier effect (the counts are binomial, as in the Random Degree Model below), larger values mean stronger barrier effects with respect to subpopulation $k$. Because the degrees are explicit parameters, no rescaling is needed: `fit$summary("d")` directly gives posterior summaries of the degree of every respondent, and `exp(mu + sigma^2 / 2)` is the mean degree in the population. ## Count responses The inputs are `N`, `K` and the count matrix `y` as before, together with `m`, the fractional sizes of the subpopulations, and `L`, the per-respondent lower bounds on the degrees. The remaining arguments are the sampling settings: ``` r fit <- maltiel_bem_count( N = N, K = K, y = y, m = m, L = L, chains = 1, iter_warmup = 100, iter_sampling = 100, refresh = 0, show_messages = FALSE, show_exceptions = FALSE, seed = 666 ) fit #> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail #> lp__ -249.05 -249.30 7.26 8.35 -261.37 -238.90 1.21 4 76 #> d_raw[1] 22.79 21.61 7.44 6.54 13.84 35.70 1.09 109 48 #> d_raw[2] 21.01 20.93 5.32 5.38 12.76 28.76 0.99 200 66 #> d_raw[3] 16.76 16.43 4.42 4.41 10.92 24.48 1.01 128 71 #> d_raw[4] 20.85 19.94 5.53 4.73 12.40 29.15 1.06 200 48 #> d_raw[5] 19.04 18.30 6.79 6.29 10.39 30.65 1.01 200 92 #> d_raw[6] 18.85 18.72 4.89 4.50 12.08 27.95 1.00 174 76 #> d_raw[7] 18.90 18.79 4.92 4.56 12.19 27.28 1.01 128 78 #> d_raw[8] 19.25 18.35 5.37 5.95 11.80 28.67 1.00 186 77 #> d_raw[9] 18.43 17.51 4.75 4.33 13.47 25.59 1.03 120 76 #> #> # showing 10 of 106 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option) ``` A quick look at the degree distribution and the barrier effect parameters: ``` r fit$summary(variables = c("mu", "sigma", "rho")) #> # A tibble: 5 × 10 #> variable mean median sd mad q5 q95 rhat ess_bulk #> #> 1 mu 3.08 3.06 0.0617 0.0584 3.00 3.19 0.992 20.9 #> 2 sigma 0.291 0.278 0.0445 0.0412 0.250 0.384 1.66 1.79 #> 3 rho[1] 0.0184 0.0148 0.0130 0.0122 0.00461 0.0434 1.07 10.5 #> 4 rho[2] 0.00424 0.00167 0.00553 0.00243 0.0000137 0.0161 1.41 2.47 #> 5 rho[3] 0.00745 0.00393 0.00813 0.00481 0.000404 0.0221 0.998 15.0 #> # ℹ 1 more variable: ess_tail head(fit$summary("d")) #> # A tibble: 6 × 10 #> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail #> #> 1 d[1] 25.8 24.6 7.44 6.54 16.8 38.7 1.09 110. 48.4 #> 2 d[2] 26.0 25.9 5.32 5.38 17.8 33.8 0.993 200 67.0 #> 3 d[3] 18.8 18.4 4.42 4.41 12.9 26.5 1.01 128. 71.3 #> 4 d[4] 22.9 21.9 5.53 4.73 14.4 31.2 1.06 200 48.2 #> 5 d[5] 21.0 20.3 6.79 6.29 12.4 32.6 1.01 200 92.6 #> 6 d[6] 21.8 21.7 4.89 4.50 15.1 31.0 1.00 174. 76.8 ``` The data were simulated with `mu = 3`, which is also the lower bound of the prior on `mu`, so the posterior piles up against that bound; likewise the true `sigma` (1/16) is below the lower limit of its prior, and `rho` should be small since the data contain no barrier effects. ## Dichotomous responses The dichotomous version uses the beta-binomial probabilities of knowing nobody versus at least one person. It does not constrain the degrees from below, so `L` is not used by the Stan model, but the wrapper accepts it for consistency with the count version. The inputs are `N`, `K`, the 0/1 matrix `y01`, `m` and `L`, followed by the sampling settings: ``` r fit <- maltiel_bem_dichotomous( N = N, K = K, y = y01, m = m, L = L, chains = 1, iter_warmup = 100, iter_sampling = 100, refresh = 0, show_messages = FALSE, show_exceptions = FALSE, seed = 666 ) #> Warning: 1 of 1 chains had an E-BFMI less than 0.3. #> See https://mc-stan.org/misc/warnings for details. fit #> variable mean median sd mad q5 q95 rhat ess_bulk #> lp__ -157.50 -155.50 21.44 25.39 -187.77 -120.78 1.84 1 #> d[1] 2389.41 144.58 6436.69 165.79 32.40 10601.21 1.82 1 #> d[2] 2074.20 135.77 4486.31 137.73 27.14 11874.01 1.77 1 #> d[3] 263.37 31.71 784.62 28.39 7.10 1303.60 1.38 2 #> d[4] 2166.86 176.88 8759.35 214.32 23.82 5225.35 1.78 1 #> d[5] 695.97 127.00 1460.25 146.43 20.16 2860.69 1.74 1 #> d[6] 1268.28 82.30 4986.20 87.35 21.76 3608.46 1.76 1 #> d[7] 1186.33 133.84 2845.13 164.42 13.57 6174.97 1.62 1 #> d[8] 1688.52 109.80 7120.70 108.02 17.14 4844.39 1.51 2 #> d[9] 2044.52 113.38 10057.82 127.52 14.18 8392.27 1.61 1 #> ess_tail #> 32 #> 30 #> 33 #> 48 #> 31 #> 73 #> 32 #> 38 #> 60 #> 60 #> #> # showing 10 of 56 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option) ``` ## Trichotomous responses The trichotomous version distinguishes nobody / one person / two or more and, like the dichotomous one, ignores `L`. The inputs are `N`, `K`, the 0/1/2 matrix `y012`, `m` and `L`, followed by the sampling settings: ``` r fit <- maltiel_bem_trichotomous( N = N, K = K, y = y012, m = m, L = L, chains = 1, iter_warmup = 100, iter_sampling = 100, refresh = 0, show_messages = FALSE, show_exceptions = FALSE, seed = 666 ) fit #> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail #> lp__ -185.96 -185.50 10.84 12.22 -201.35 -168.53 1.06 11 30 #> d[1] 36.37 33.13 17.10 14.39 15.78 60.32 1.01 124 48 #> d[2] 27.64 26.23 9.60 8.44 15.83 44.43 1.00 106 53 #> d[3] 23.29 22.66 9.12 9.17 9.66 37.57 1.10 110 56 #> d[4] 32.93 31.24 15.06 10.26 16.38 64.42 1.05 102 42 #> d[5] 26.85 24.01 14.49 10.94 13.04 45.81 1.03 110 15 #> d[6] 25.99 25.73 8.93 8.09 13.37 40.03 1.04 88 59 #> d[7] 27.63 23.89 13.82 13.85 11.35 52.35 1.02 120 24 #> d[8] 33.14 29.48 13.36 9.61 18.42 52.12 1.01 120 78 #> d[9] 27.21 26.44 11.67 9.00 12.90 45.50 1.04 18 66 #> #> # showing 10 of 56 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option) ``` # Maltiel et al (2015) Random Degree Model See `?maltiel_rdm`. This is the simplest of the Maltiel et al (2015) models: the count $y_{ik}$ is binomial with $d_i$ trials and success probability $m_k$, so the only parameters are the degrees `d` and the lognormal hyperparameters `mu` and `sigma`. It is also the model the example data were simulated from, so the fit should recover `mu` close to 3. ## Count responses The inputs are `N`, `K`, the count matrix `y`, the fractional subpopulation sizes `m` and the degree lower bounds `L`, followed by the sampling settings: ``` r fit <- maltiel_rdm_count( N = N, K = K, y = y, m = m, L = L, chains = 1, iter_warmup = 100, iter_sampling = 100, refresh = 0, show_messages = FALSE, show_exceptions = FALSE, seed = 666 ) fit #> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail #> lp__ -227.79 -227.65 5.64 4.92 -238.16 -219.42 1.12 9 30 #> d_raw[1] 23.30 23.23 6.08 5.07 14.31 33.76 1.07 120 52 #> d_raw[2] 21.43 20.72 5.62 5.84 13.43 30.78 1.03 158 112 #> d_raw[3] 17.00 16.51 4.65 4.07 9.39 25.29 1.05 200 71 #> d_raw[4] 20.64 20.05 5.18 5.85 13.51 30.04 1.00 200 61 #> d_raw[5] 18.73 18.59 5.46 4.40 9.95 29.78 1.00 200 101 #> d_raw[6] 18.93 19.17 5.85 5.29 10.50 29.91 0.99 200 90 #> d_raw[7] 19.20 18.89 6.21 5.46 9.68 27.46 1.02 200 65 #> d_raw[8] 19.57 19.06 4.60 5.05 13.00 27.36 1.00 200 116 #> d_raw[9] 18.29 17.38 4.71 4.52 11.80 26.78 1.03 200 78 #> #> # showing 10 of 103 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option) ``` The summary starts with `d_raw`, the increments over the lower bounds `L` that the sampler actually works with; the degrees `d = L + d_raw` are further down the list. The hyperparameters of the degree distribution: ``` r fit$summary(variables = c("mu", "sigma")) #> # A tibble: 2 × 10 #> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail #> #> 1 mu 3.08 3.08 0.0499 0.0403 3.00 3.18 1.07 21.3 21.0 #> 2 sigma 0.285 0.276 0.0352 0.0259 0.253 0.332 1.13 7.11 30.6 ``` ## Dichotomous responses The dichotomous version (note the wrapper name, `maltiel_rd_dichotomous()`) uses the binomial probabilities of a zero versus a positive count and does not use `L`. The inputs are `N`, `K`, the 0/1 matrix `y01`, `m` and `L`, followed by the sampling settings: ``` r fit <- maltiel_rd_dichotomous( N = N, K = K, y = y01, m = m, L = L, chains = 1, iter_warmup = 100, iter_sampling = 100, refresh = 0, show_messages = FALSE, show_exceptions = FALSE, seed = 666 ) #> Warning: 1 of 1 chains had an E-BFMI less than 0.3. #> See https://mc-stan.org/misc/warnings for details. fit #> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail #> lp__ -104.21 -103.59 12.97 16.88 -125.26 -84.08 1.35 2 17 #> d[1] 31.01 27.75 12.59 12.31 16.71 58.50 1.07 60 112 #> d[2] 33.83 26.62 25.39 11.54 14.11 56.80 1.18 4 21 #> d[3] 17.26 16.97 7.09 6.44 6.13 31.08 1.03 66 112 #> d[4] 33.28 27.24 22.41 11.56 14.12 67.84 1.05 28 25 #> d[5] 23.92 21.93 9.21 5.76 12.41 44.57 0.99 173 68 #> d[6] 24.00 21.95 8.52 8.25 13.63 40.97 0.99 166 101 #> d[7] 23.30 21.99 9.05 7.70 9.76 38.48 0.99 200 71 #> d[8] 23.12 21.60 8.38 7.82 11.85 38.50 1.07 81 47 #> d[9] 23.69 21.88 9.21 6.75 13.07 40.17 1.00 155 76 #> #> # showing 10 of 53 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option) ``` # Maltiel et al (2015) Transmission Bias Model See `?maltiel_tbm`. Respondents may not know that some of their contacts belong to a subpopulation -- membership in stigmatized or private groups is not always "transmitted" across ties. The Transmission Bias Model extends the Random Degree Model by replacing the prevalence $m_k$ with $\tau_k m_k$, where `tau[k]` in (0, 1) is the fraction of the members of subpopulation $k$ whose membership is visible to their contacts. Each `tau[k]` gets a beta prior with shape parameters `eta[k]` and `v[k]` supplied by the user. Since the expected counts depend on the product $d_i \tau_k m_k$, the data alone cannot separate a low `tau[k]` from a small `m[k]`, nor either of them from the degrees: inflating the degrees and deflating the biases leaves the fitted counts unchanged. These priors are therefore where knowledge about transmission (e.g. from other studies) enters the model: `eta = 1, v = 1` is a uniform prior, whereas `eta = 9, v = 1` encodes a belief that about 90% of the members are visible. ## Count responses The inputs are `N`, `K`, the count matrix `y`, the fractional subpopulation sizes `m` and the degree lower bounds `L` as in the Random Degree Model, plus the vectors `eta` and `v` of length `K` with the beta prior shape parameters for each `tau[k]` -- here uniform priors for all three subpopulations. The remaining arguments are the sampling settings: ``` r fit <- maltiel_tbm_count( N = N, K = K, y = y, m = m, L = L, eta = rep(1, K), v = rep(1, K), chains = 1, iter_warmup = 100, iter_sampling = 100, refresh = 0, show_messages = FALSE, show_exceptions = FALSE, seed = 666 ) fit #> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail #> lp__ -234.24 -234.32 7.34 8.51 -246.02 -222.72 0.99 25 72 #> d_raw[1] 31.04 29.51 11.51 10.73 14.60 47.63 1.05 41 36 #> d_raw[2] 30.73 28.09 12.73 10.90 16.24 56.12 1.00 26 13 #> d_raw[3] 23.41 20.05 11.20 7.10 11.00 44.74 1.00 22 31 #> d_raw[4] 28.13 25.26 10.82 8.34 16.33 53.45 1.02 20 28 #> d_raw[5] 24.59 22.83 9.45 7.56 13.84 43.87 1.00 23 31 #> d_raw[6] 24.61 23.26 8.57 7.05 13.11 39.33 1.01 24 42 #> d_raw[7] 26.19 23.16 11.21 8.28 14.44 46.25 1.00 23 31 #> d_raw[8] 26.77 23.84 12.38 7.02 13.83 49.02 1.00 29 52 #> d_raw[9] 26.10 23.23 12.39 9.34 12.17 51.04 1.00 38 30 #> #> # showing 10 of 106 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option) ``` The transmission parameters are summarised as usual: ``` r fit$summary("tau") #> # A tibble: 3 × 10 #> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail #> #> 1 tau[1] 0.696 0.696 0.158 0.179 0.429 0.941 1.04 24.2 12.8 #> 2 tau[2] 0.825 0.876 0.142 0.0904 0.520 0.966 1.02 12.9 12.7 #> 3 tau[3] 0.750 0.782 0.174 0.179 0.401 0.966 1.05 11.0 12.9 ``` Without informative priors `tau` and the degrees `d` are only weakly separated -- a lower visibility of a subpopulation can be compensated by larger degrees -- so the posterior of `tau` remains wide. The data were simulated with full transmission, i.e. `tau = 1`. # Maltiel et al (2015) Combined Model See `?maltiel_cm`. This is the model of section 2.4 of @maltiel_etal_2015, carrying both the barrier effects of the Barrier Effects Model and the transmission bias of the Transmission Bias Model: $$ y_{ik} \sim \mathrm{Binomial}(d_i, \tau_k q_{ik}), \qquad d_i \sim \mathrm{LogNormal}(\mu, \sigma^2), \qquad q_{ik} \sim \mathrm{Beta}(m_k, \rho_k) $$ where $\mathrm{Beta}(m_k, \rho_k)$ is parametrized by its mean $m_k$ and overdispersion $\rho_k$, as in the Barrier Effects Model. Multiplying $q_{ik}$ by $\tau_k$ breaks the beta-binomial conjugacy that `maltiel_bem_count()` relies on, so the latent $q_{ik}$ cannot be integrated out and are sampled instead. The model therefore has `N * K` more parameters than the others and is the slowest in the package. ## Count responses The inputs are those of the Transmission Bias Model -- `N`, `K`, the count matrix `y`, the fractional subpopulation sizes `m`, the degree lower bounds `L` and the beta prior shapes `eta` and `v` for the transmission biases -- followed by the sampling settings. Because the `tau`/degree confounding described above applies here too, and because this model has the barrier effects to fit as well, we use an informative `Beta(9, 1)` prior (mean 0.9) rather than the flat one: ``` r fit <- maltiel_cm_count( N = N, K = K, y = y, m = m, L = L, eta = rep(9, K), v = rep(1, K), chains = 1, iter_warmup = 100, iter_sampling = 100, refresh = 0, show_messages = FALSE, show_exceptions = FALSE, seed = 666 ) #> Warning: 1 of 1 chains had an E-BFMI less than 0.3. #> See https://mc-stan.org/misc/warnings for details. fit #> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail #> lp__ -327.34 -337.89 55.55 73.54 -402.21 -243.36 1.99 1 21 #> d_raw[1] 25.54 24.51 8.25 6.51 15.57 39.52 1.05 138 112 #> d_raw[2] 22.68 22.24 6.50 6.69 13.18 33.38 1.00 96 78 #> d_raw[3] 17.40 16.75 5.85 5.25 9.77 27.30 1.01 66 31 #> d_raw[4] 22.13 21.23 6.59 6.41 13.27 35.20 1.01 95 24 #> d_raw[5] 20.49 19.13 6.48 5.58 12.52 31.65 1.02 69 37 #> d_raw[6] 20.55 19.47 5.28 4.65 12.91 29.07 1.00 96 76 #> d_raw[7] 19.14 18.13 6.56 5.16 9.84 31.39 1.00 98 76 #> d_raw[8] 21.85 20.92 5.72 6.20 14.28 32.38 1.00 140 53 #> d_raw[9] 19.82 18.71 5.49 5.90 13.54 29.75 1.01 115 74 #> #> # showing 10 of 259 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option) ``` The summary is dominated by the `N * K` latent `q` values. The two sets of parameters that describe the subpopulations are more informative: ``` r fit$summary(variables = c("tau", "rho")) #> # A tibble: 6 × 10 #> variable mean median sd mad q5 q95 rhat ess_bulk #> #> 1 tau[1] 0.883 0.907 0.0848 0.0812 0.718 0.982 1.01 68.7 #> 2 tau[2] 0.940 0.952 0.0510 0.0443 0.848 0.998 1.03 21.6 #> 3 tau[3] 0.893 0.906 0.0705 0.0686 0.770 0.984 1.04 35.6 #> 4 rho[1] 0.0192 0.0130 0.0160 0.00949 0.00278 0.0511 1.48 2.18 #> 5 rho[2] 0.00769 0.00703 0.00606 0.00653 0.000553 0.0177 1.22 3.37 #> 6 rho[3] 0.0141 0.0104 0.0109 0.00874 0.00271 0.0364 1.59 1.86 #> # ℹ 1 more variable: ess_tail ``` `tau` stays close to the prior, as the discussion above leads one to expect, and `rho` measures the barrier effects on top of it. With so few iterations and this many parameters the sampler has not explored the posterior properly; a real fit of this model needs considerably longer runs than the other models in this vignette. # Where to go from here The examples above show the calling syntax and the shape of the results. In a real analysis: - run several chains for many more iterations, check `rhat`, the effective sample sizes and `$diagnostic_summary()` before trusting any summary; - use `$draws()` with the **posterior** package to compute derived quantities (degrees, subpopulation sizes, mean degree) for every draw and summarize them, and packages such as **bayesplot** to visualize the posterior; - save the fitted object with `$save_object()` because the underlying CSV files are temporary. The help pages `?zheng_bem`, `?zheng_gp`, `?maltiel_bem`, `?maltiel_rdm`, `?maltiel_tbm` and `?maltiel_cm` document the wrappers and point to the original papers. # References