mcmcARD: Estimate network model using ARD

In ahoundetoungan/PartialNetwork: Estimating Peer Effects Using Partial Network Data

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()
)

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.

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.

Examples

  # Sample size
  N       <- 500
  
  # 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
  out   <- mcmcARD(Y = ARD, traitARD = trait, start = start, fixv = vfixcolumn,
                   consb = bfixcolumn, iteration = 5000)
  
  # plot simulations
  # plot d
  plot(out$simulations$d[,100], type = "l", col = "blue", ylab = "")
  abline(h = gend[100], col = "red")
  
  # plot coordinates of individuals
  i <- 123 # 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")
    })
  }

ahoundetoungan/PartialNetwork documentation built on March 15, 2024, 4:35 p.m.