## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(networkscaleup)

## ----simulation---------------------------------------------------------------
set.seed(1998)
N_i <- 50
N_k <- 5
N <- 1e5
sizes <- rbinom(N_k, N, prob = runif(N_k, min = 0.01, max = 0.15))
degrees <- round(exp(rnorm(N_i, mean = 5, sd = 1)))

ard <- matrix(NA, nrow = N_i, ncol = N_k)
for (k in 1:N_k) {
  ard[, k] <- rbinom(N_i, degrees, sizes[k] / N)
}

## Create some artificial covariates for use later
x <- matrix(sample(1:5, size = N_i * N_k, replace = T),
  nrow = N_i,
  ncol = N_k
)
z <- cbind(rbinom(N_i, 1, 0.3), rnorm(N_i), rnorm(N_i), rnorm(N_i))

## ----prep---------------------------------------------------------------------
x <- scale(x)
z <- scale(z)
z_subpop <- z[, 1:2]
z_global <- z[, 3:4]

## ----pimle--------------------------------------------------------------------
pimle.est <- killworth(ard,
  known_sizes = sizes[c(1, 2, 4)],
  known_ind = c(1, 2, 4),
  N = N, model = "PIMLE"
)


plot(degrees ~ pimle.est$degrees, xlab = "Estimated PIMLE degrees", ylab = "True Degrees")
abline(0, 1, col = "red")

round(data.frame(
  true = sizes[c(3, 5)],
  pimle = pimle.est$sizes
))

## ----mle----------------------------------------------------------------------
mle.est <- killworth(ard,
  known_sizes = sizes[c(1, 2, 4)],
  known_ind = c(1, 2, 4),
  N = N, model = "MLE"
)


plot(degrees ~ mle.est$degrees, xlab = "Estimated MLE degrees", ylab = "True Degrees")
abline(0, 1, col = "red")

round(data.frame(
  true = sizes[c(3, 5)],
  pimle = mle.est$sizes
))

## ----overdisp-----------------------------------------------------------------
overdisp_gibbs_metrop_est <- overdispersed(
  ard,
  known_sizes = sizes[c(1, 2, 4)],
  known_ind = c(1, 2, 4),
  G1_ind = 1,
  G2_ind = 2,
  B2_ind = 4,
  N = N,
  warmup = 500,
  iter = 1000,
  verbose = TRUE,
  init = "MLE"
)

overdisp_stan <- overdispersedStan(
  ard,
  known_sizes = sizes[c(1, 2, 4)],
  known_ind = c(1, 2, 4),
  G1_ind = 1,
  G2_ind = 2,
  B2_ind = 4,
  N = N,
  chains = 2,
  cores = 2,
  warmup = 250,
  iter = 500,
)


round(data.frame(
  true = sizes,
  gibbs_est = colMeans(overdisp_gibbs_metrop_est$sizes),
  stan_est = colMeans(overdisp_stan$sizes)
))

plot(degrees ~ colMeans(overdisp_stan$degrees), xlab = "Overdispersed Degree Estimates", ylab = "True Degrees")
abline(0, 1, col = "red")

## ----correlated---------------------------------------------------------------
correlated_cov_stan <- correlatedStan(
  ard,
  known_sizes = sizes[c(1, 2, 4)],
  known_ind = c(1, 2, 4),
  model = "correlated",
  scaling = "weighted",
  x = x,
  z_subpop = z_subpop,
  z_global = z_global,
  N = N,
  chains = 2,
  cores = 2,
  warmup = 250,
  iter = 500,
)

correlated_nocov_stan <- correlatedStan(
  ard,
  known_sizes = sizes[c(1, 2, 4)],
  known_ind = c(1, 2, 4),
  model = "correlated",
  scaling = "all",
  N = N,
  chains = 2,
  cores = 2,
  warmup = 250,
  iter = 500,
)


uncorrelated_cov_stan <- correlatedStan(
  ard,
  known_sizes = sizes[c(1, 2, 4)],
  known_ind = c(1, 2, 4),
  model = "uncorrelated",
  scaling = "all",
  x = x,
  z_subpop = z_subpop,
  z_global = z_global,
  N = N,
  chains = 2,
  cores = 2,
  warmup = 250,
  iter = 500,
)

uncorrelated_x_stan <- correlatedStan(
  ard,
  known_sizes = sizes[c(1, 2, 4)],
  known_ind = c(1, 2, 4),
  model = "uncorrelated",
  scaling = "all",
  x = x,
  N = N,
  chains = 2,
  cores = 2,
  warmup = 250,
  iter = 500,
)


round(data.frame(
  true = sizes,
  corr_cov_est = colMeans(correlated_cov_stan$sizes),
  corr_nocov_est = colMeans(correlated_nocov_stan$sizes),
  uncorr_cov_est = colMeans(uncorrelated_cov_stan$sizes),
  uncorr_x_est = colMeans(uncorrelated_x_stan$sizes)
))

plot(degrees ~ colMeans(correlated_cov_stan$degrees), xlab = "Correlated Covariate Degree Estimates", ylab = "True Degrees")
abline(0, 1, col = "red")

## Examine parameter estimates
colMeans(correlated_cov_stan$alpha)
colMeans(correlated_cov_stan$beta_global)
colMeans(correlated_cov_stan$beta_subpop)