mcmcARD: Estimate network model using ARD

Description

mcmcARD estimates the network model proposed by Breza et al. (2020).

Usage

mcmcARD(
  Y,
  traitARD,
  start,
  fixv,
  consb,
  iteration = 2000L,
  sim.d = TRUE,
  sim.zeta = TRUE,
  hyperparms = NULL,
  ctrl.mcmc = list()
)

Value

A list consisting of:

n: dimension of the sample with ARD.
K: number of traits.
p: hypersphere dimension.
time: elapsed time in second.
iteration: number of MCMC steps performed.
simulations: simulations from the posterior distribution.
hyperparms: return value of hyperparameters (updated and non updated).
accept.rate: list of acceptance rates.
start: starting values.
ctrl.mcmc: return value of ctrl.mcmc.

Arguments

Y: is a matrix of ARD. The entry (i, k) is the number of i's friends having the trait k.
traitARD: is the matrix of traits for individuals with ARD. The entry (i, k) is equal to 1 if i has the trait k and 0 otherwise.
start: is a list containing starting values of z (matrix of dimension $N \times p$), v (matrix of dimension $K \times p$), d (vector of dimension $N$), b (vector of dimension $K$), eta (vector of dimension $K$) and zeta (scalar).
fixv: is a vector setting which location parameters are fixed for identifiability. These fixed positions are used to rotate the latent surface back to a common orientation at each iteration using a Procrustes transformation (see Section Identification in Details).
consb: is a vector of the subset of $\beta_k$ constrained to the total size (see Section Identification in Details).
iteration: is the number of MCMC steps to be performed.
sim.d: is logical indicating whether the degree d will be updated in the MCMC. If sim.d = FALSE, the starting value of d in the argument start is set fixed along the MCMC.
sim.zeta: is logical indicating whether the degree zeta will be updated in the MCMC. If sim.zeta = FALSE, the starting value of zeta in the argument start is set fixed along the MCMC.
hyperparms: is an 8-dimensional vector of hyperparameters (in this order) $\mu_d$, $\sigma_d$, $\mu_b$, $\sigma_b$, $\alpha_{\eta}$, $\beta_{\eta}$, $\alpha_{\zeta}$ and $\beta_{\zeta}$ (see Section Model in Details).
ctrl.mcmc: is a list of MCMC controls (see Section MCMC control in Details).

Details

The linking probability is given by

Model

$$P_{ij} \propto \exp(\nu_i + \nu_j + \zeta\mathbf{z}_i\mathbf{z}_j).$$ McCormick and Zheng (2015) write the likelihood of the model with respect to the spherical coordinate $\mathbf{z}_i$, the trait locations $\mathbf{v}_k$, the degree $d_i$, the fraction of ties in the network that are made with members of group k $b_k$, the trait intensity parameter $\eta_k$ and $\zeta$. The following prior distributions are defined. $$\mathbf{z}_i \sim Uniform ~ von ~ Mises-Fisher$$ $$\mathbf{v}_k \sim Uniform ~ von ~ Mises-Fisher$$ $$d_i \sim log-\mathcal{N}(\mu_d, \sigma_d)$$ $$b_k \sim log-\mathcal{N}(\mu_b, \sigma_b)$$ $$\eta_k \sim Gamma(\alpha_{\eta}, \beta_{\eta})$$ $$\zeta \sim Gamma(\alpha_{\zeta}, \beta_{\zeta})$$

Identification

For identification, some $\mathbf{v}_k$ and $b_k$ need to be exogenously fixed around their given starting value (see McCormick and Zheng, 2015 for more details). The parameter fixv can be used to set the desired value for $\mathbf{v}_k$ while fixb can be used to set the desired values for $b_k$.

MCMC control

During the MCMC, the jumping scales are updated following Atchade and Rosenthal (2005) in order to target the acceptance rate of each parameter to the target values. This requires to set minimal and maximal jumping scales through the parameter ctrl.mcmc. The parameter ctrl.mcmc is a list which can contain the following named components.

target: The default value is rep(0.44, 5). The target of every $\mathbf{z}_i$, $d_i$, $b_k$, $\eta_k$ and $\zeta$ is 0.44.
jumpmin: The default value is c(0,1,1e-7,1e-7,1e-7)*1e-5. The minimal jumping of every $\mathbf{z}_i$ is 0, every $d_i$ is $10^{-5}$, and every $b_k$, $\eta_k$ and $\zeta$ is $10^{-12}$.
jumpmax: The default value is c(100,1,1,1,1)*20. The maximal jumping scale is 20 except for $\mathbf{z}_i$ which is set to 2000.
print: A logical value which indicates if the MCMC progression should be printed in the console. The default value is TRUE.

Examples

Run this code

# Sample size
N       <- 30 #for a very small sample

# ARD parameters
genzeta <- 1
mu      <- -1.35
sigma   <- 0.37
K       <- 12    # number of traits
P       <- 3     # Sphere dimension


# Generate z (spherical coordinates)
genz    <- rvMF(N,rep(0,P))

# Generate nu  from a Normal distribution with parameters mu and sigma (The gregariousness)
gennu   <- rnorm(N,mu,sigma)

# compute degrees
gend <- N*exp(gennu)*exp(mu+0.5*sigma^2)*exp(logCpvMF(P,0) - logCpvMF(P,genzeta))

# Link probabilities
Probabilities <- sim.dnetwork(gennu,gend,genzeta,genz)

# Adjacency matrix
G <- sim.network(Probabilities)

# Generate vk, the trait location
genv <- rvMF(K,rep(0,P))

# set fixed some vk  distant
genv[1,] <- c(1,0,0)
genv[2,] <- c(0,1,0)
genv[3,] <- c(0,0,1)

# eta, the intensity parameter
geneta   <-abs(rnorm(K,2,1))

# Build traits matrix
densityatz       <- matrix(0,N,K)
for(k in 1:K){
  densityatz[,k] <- dvMF(genz,genv[k,]*geneta[k])
}

trait       <- matrix(0,N,K)
NK          <- floor(runif(K, 0.8, 0.95)*colSums(densityatz)/apply(densityatz, 2, max))
for (k in 1:K) {
  trait[,k]  <- rbinom(N, 1, NK[k]*densityatz[,k]/sum(densityatz[,k]))
}

# print a percentage of people having a trait
colSums(trait)*100/N

# Build ARD
ARD         <- G %*% trait

# generate b
genb        <- numeric(K)
for(k in 1:K){
  genb[k]   <- sum(G[,trait[,k]==1])/sum(G)
}

############ ARD Posterior distribution ###################
# initialization
d0     <- exp(rnorm(N)); b0 <- exp(rnorm(K)); eta0 <- rep(1,K);
zeta0  <- 05; z0 <- matrix(rvMF(N,rep(0,P)),N); v0 <- matrix(rvMF(K,rep(0,P)),K)

# We need to fix some of the vk and bk for identification (see Breza et al. (2020) for details).
vfixcolumn      <- 1:6
bfixcolumn      <- c(3, 5)
b0[bfixcolumn]  <- genb[bfixcolumn]
v0[vfixcolumn,] <- genv[vfixcolumn,]
start  <- list("z" = z0, "v" = v0, "d" = d0, "b" = b0, "eta" = eta0, "zeta" = zeta0)

# MCMC
# REMOVE COMMENTS TO RUN THE REST OF THE CODE.
# mcmcARD USES Rcpp FUNCTIONS IN PARALLEL TO BE FAST AND CRAN POLICY DOES 
# NOT ALLOW CODE IN PARALLEL. FOR THIS REASON WE PUT THE REST OF THE CODE IN 
# PARALLEL SO THAT CRAN ACCEPTS THE PACKAGE.
# out   <- mcmcARD(Y = ARD, traitARD = trait, start = start, fixv = vfixcolumn,
#                  consb = bfixcolumn, iteration = 500)
# 
# # plot simulations
# # plot d
# plot(out$simulations$d[,10], type = "l", col = "blue", ylab = "")
# abline(h = gend[10], col = "red")
# 
# # plot coordinates of individuals
# i <- 13 # individual 123
# {
#   lapply(1:3, function(x) {
#     plot(out$simulations$z[i, x,] , type = "l", ylab = "", col = "blue", ylim = c(-1, 1))
#     abline(h = genz[i, x], col = "red")
#   })
# }
# 
# # plot coordinates of traits
# k <- 8
# {
#   lapply(1:3, function(x) {
#     plot(out$simulations$v[k, x,] , type = "l", ylab = "", col = "blue", ylim = c(-1, 1))
#     abline(h = genv[k, x], col = "red")
#   })
# }

Run the code above in your browser using DataLab