R/mcmcARD.R

In PartialNetwork: Estimating Peer Effects Using Partial Network Data

#' @title Estimate network model using ARD
#' @description  \code{mcmcARD} estimates the network model proposed by Breza et al. (2020).
#' @param Y is a matrix of ARD. The entry (i, k) is the number of i's friends having the trait k.
#' @param 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.
#' @param start is a list containing starting values of `z` (matrix of dimension \eqn{N \times p}), `v` (matrix of dimension \eqn{K \times p}),
#'  `d` (vector of dimension \eqn{N}), `b` (vector of dimension \eqn{K}), `eta` (vector of dimension \eqn{K}) and `zeta` (scalar).
#' @param 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).
#' @param consb is a vector of the subset of \eqn{\beta_k}{bk} constrained to the total size (see Section Identification in Details).
#' @param iteration is the number of MCMC steps to be performed.
#' @param 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.
#' @param 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.
#' @param hyperparms is an 8-dimensional vector of hyperparameters (in this order) \eqn{\mu_d}{mud},  \eqn{\sigma_d}{sigmad},
#' \eqn{\mu_b}{mub}, \eqn{\sigma_b}{sigmab}, \eqn{\alpha_{\eta}}{alphaeta}, \eqn{\beta_{\eta}}{betaeta}, 
#' \eqn{\alpha_{\zeta}}{alphazeta} and \eqn{\beta_{\zeta}}{betazeta} (see Section Model in Details).
#' @param ctrl.mcmc is a list of MCMC controls (see Section MCMC control in Details).
#' 
#' @details The linking probability is given by
#' ## Model
#' \deqn{P_{ij} \propto \exp(\nu_i + \nu_j + \zeta\mathbf{z}_i\mathbf{z}_j).}{Pij is proportional to exp(nui + nuj + zeta * zi * zj).}
#' McCormick and Zheng (2015) write the likelihood of the model with respect to the spherical coordinate \eqn{\mathbf{z}_i}{zi},
#' the trait locations \eqn{\mathbf{v}_k}{vk}, the degree \eqn{d_i}{di}, the fraction of ties in the network that are
#' made with members of group k \eqn{b_k}{bk}, the trait intensity parameter \eqn{\eta_k}{etak} and \eqn{\zeta}{zeta}. 
#' The following
#' prior distributions are defined.
#' \deqn{\mathbf{z}_i \sim Uniform ~ von ~ Mises-Fisher}{zi ~ Uniform von Mises Fisher}
#' \deqn{\mathbf{v}_k \sim Uniform ~ von ~ Mises-Fisher}{vk ~ Uniform von Mises Fisher}
#' \deqn{d_i \sim log-\mathcal{N}(\mu_d, \sigma_d)}{di ~ log-Normal(mud, sigmad)}
#' \deqn{b_k \sim log-\mathcal{N}(\mu_b, \sigma_b)}{bk ~ log-Normal(mub, sigmab)}
#' \deqn{\eta_k \sim Gamma(\alpha_{\eta}, \beta_{\eta})}{etak ~ Gamma(alphaeta, betaeta)}
#' \deqn{\zeta \sim Gamma(\alpha_{\zeta}, \beta_{\zeta})}{zeta ~ Gamma(alphazeta, betazeta)} 
#' ## Identification
#' For identification, some \eqn{\mathbf{v}_k}{vk} and \eqn{b_k}{bk} 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 \eqn{\mathbf{v}_k}{vk} while `fixb` can be used to set the desired values for \eqn{b_k}{bk}.\cr
#' ## 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.
#' \itemize{
#' \item{`target`}: The default value is \code{rep(0.44, 5)}. 
#' The target of every \eqn{\mathbf{z}_i}{zi}, \eqn{d_i}{di}, \eqn{b_k}{bk}, \eqn{\eta_k}{etak} and \eqn{\zeta}{zeta} is  0.44.
#' \item{`jumpmin`}: The default value is \code{c(0,1,1e-7,1e-7,1e-7)*1e-5}. 
#' The minimal jumping of every \eqn{\mathbf{z}_i}{zi} is 0, every \eqn{d_i}{di} is \eqn{10^{-5}}{1e-5}, and every \eqn{b_k}{bk}, \eqn{\eta_k}{etak} and \eqn{\zeta}{zeta} is \eqn{10^{-12}}{1e-12}.
#' \item{`jumpmax`}: The default value is \code{c(100,1,1,1,1)*20}. The maximal jumping scale is 20 except for \eqn{\mathbf{z}_i}{zi} which is set to 2000.
#' \item{`print`}: A logical value which indicates if the MCMC progression should be printed in the console. The default value is `TRUE`.
#' }
#' @return A list consisting of:
#'     \item{n}{dimension of the sample with ARD.}
#'     \item{K}{number of traits.}
#'     \item{p}{hypersphere dimension.}
#'     \item{time}{elapsed time in second.}
#'     \item{iteration}{number of MCMC steps performed.}
#'     \item{simulations}{simulations from the posterior distribution.}
#'     \item{hyperparms}{return value of hyperparameters (updated and non updated).}
#'     \item{accept.rate}{list of acceptance rates.}
#'     \item{start}{starting values.}
#'     \item{ctrl.mcmc}{return value of `ctrl.mcmc`.}
#' @examples 
#' # 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")
#' #   })
#' # }
#' @export
mcmcARD        <- function(Y, traitARD, start, fixv, consb, iteration = 2000L, sim.d = TRUE, sim.zeta = TRUE, hyperparms = NULL, ctrl.mcmc = list()) {
  t1           <- Sys.time()
  
  # hyperparameters
  if (is.null(hyperparms)) {
    hyperparms <- c(0,1,0,1,5,0.5,1,1)
  } else {
    if(length(hyperparms) != 8) {
      stop("hyperparms should be a 8-dimensional vector")
    }
  }
  
  # MCMC control
  target       <- ctrl.mcmc$target
  jumpmin      <- ctrl.mcmc$jumpmin
  jumpmax      <- ctrl.mcmc$jumpmax
  print        <- ctrl.mcmc$print
  c            <- ctrl.mcmc$c
  
  if (is.null(target)) {
    target     <- c(1,1,1,1,1)*0.44
  } else {
    if(length(target) != 5) {
      stop("target in ctrl.mcmc should be a 5-dimensional vector")
    }
  }
  if (is.null(jumpmin)) {
    jumpmin    <- c(0,1,1e-7,1e-7,1e-7)*1e-5
  } else {
    if(length(jumpmin) != 5) {
      stop("jumpmin in ctrl.mcmc should be a 5-dimensional vector")
    }
  }
  if (is.null(jumpmax)) {
    jumpmax    <- c(100,1,1,1,1)*20
  } else {
    if(length(jumpmax) != 5) {
      stop("jumpmax in ctrl.mcmc should be a 5-dimensional vector")
    }
  }
  if (is.null(c)) {
    c          <- 0.6
  } else {
    if(length(c) != 1) {
      stop("c in ctrl.mcmc should be a scalar")
    }
  }
  if (is.null(print)) {
    print     <- TRUE
  } else {
    if(length(print) != 1) {
      stop("print in ctrl.mcmc should be a scalar")
    }
  }
  namjum         <- c("z", "d", "b", "eta", "zeta")
  names(target)  <- namjum
  names(jumpmin) <- namjum
  names(jumpmax) <- namjum
  
  
  ctrl.mcmc    <- list(
    target     = target,
    jumpmin    = jumpmin,
    jumpmax    = jumpmax,
    print      = print,
    c          = c
  )
  
  z0           <- start$z
  v0           <- start$v
  d0           <- start$d
  b0           <- start$b
  eta0         <- start$eta
  zeta0        <- start$zeta
  stopifnot(start$d > 0)
  stopifnot(start$b > 0)
  stopifnot(start$eta > 0)
  stopifnot(start$zeta > 0)
  
  if (is.null(z0) | is.null(v0) | is.null(d0) | is.null(b0) | is.null(eta0) | is.null(zeta0)) {
    stop("Elements in start should be named as z, v, b, d, eta and zeta")
  }
  
  n            <- nrow(Y)
  K            <- ncol(Y)
  p            <- ncol(z0)
  
  out          <- updateGP(Y, traitARD, z0, v0, d0, b0, eta0, zeta0, fixv, consb, iteration, !sim.d, !sim.zeta,
                           hyperparms, target, jumpmin, jumpmax,  c, print)
  
  out$accept.rate$z      <- c(out$accept.rate$z)
  out$accept.rate$d      <- c(out$accept.rate$d)
  out$accept.rate$b      <- c(out$accept.rate$b)
  out$accept.rate$eta    <- c(out$accept.rate$eta)
  
  if(!sim.d) {
    out$accept.rate$d    <- "Fixed"
  }
  if(!sim.zeta) {
    out$accept.rate$zeta <- "Fixed"
  }
  
  zaccept     <- out$accept.rate$z
  daccept     <- out$accept.rate$d
  baccept     <- out$accept.rate$b
  etaaccept   <- out$accept.rate$eta
  zetaaccept  <- out$accept.rate$zeta
  
  colnames(out$hyperparms$updated) <- c("mud", "sigmad", "mub", "sigmab")
  
  t2          <- Sys.time()
  timer       <- as.numeric(difftime(t2, t1, units = "secs"))
  
  out         <- c(list("n" = n,
                        "K" = K,
                        "p" = p,
                        "time" = timer,
                        "iteration" = iteration),
                   out,
                   list("start" = start, "ctrl.mcmc" = ctrl.mcmc))
  
  class(out)  <- "estim.ARD"
  if(print) {
    cat("\n")
    cat("The program successfully executed \n")
    cat("\n")
    cat("********SUMMARY******** \n")
    cat("n              : ", n, "\n")
    cat("K              : ", K, "\n")
    cat("Dimension      : ", p, "\n")
    cat("Iteration      : ", iteration, "\n")
    
    
    # Print the processing time
    nhours     <- floor(timer/3600)
    nminutes   <- floor((timer-3600*nhours)/60)%%60
    nseconds   <- timer-3600*nhours-60*nminutes
    cat("Elapsed time   : ", nhours, " HH ", nminutes, " mm ", round(nseconds), " ss \n \n")
    cat("Average acceptance rate \n")
    cat("                      z: ", mean(zaccept), "\n")
    cat("                      d: ", ifelse(sim.d, mean(daccept), daccept), "\n")
    cat("                      b: ", mean(baccept), "\n")
    cat("                    eta: ", mean(etaaccept), "\n")
    cat("                   zeta: ", zetaaccept, "\n")
  }
  
  out
}