R: A function that generates samples for a univariate network...

uniNetLeroux {netcmc}

R Documentation

A function that generates samples for a univariate network Leroux model.

Description

This function generates samples for a univariate network Leroux model, which is given by

Y_{i_s}|\mu_{i_s} \sim f(y_{i_s}| \mu_{i_s}, \sigma_{e}^{2}) ~~~ i=1,\ldots, N_{s},~s=1,\ldots,S ,

g(\mu_{i_s}) = \boldsymbol{x}^\top_{i_s} \boldsymbol{\beta} + \phi_{s} + \sum_{j\in \textrm{net}(i_s)}w_{i_sj}u_{j} + w^{*}_{i_s}u^{*},

\boldsymbol{\beta} \sim \textrm{N}(\boldsymbol{0}, \alpha\boldsymbol{I}),

\phi_{s} | \boldsymbol{\phi}_{-s} \sim \textrm{N}\bigg(\frac{ \rho \sum_{l = 1}^{S} a_{sl} \phi_{l} }{ \rho \sum_{l = 1}^{S} a_{sl} + 1 - \rho }, \frac{ \tau^{2} }{ \rho \sum_{l = 1}^{S} a_{sl} + 1 - \rho } \bigg),

u_{j} \sim \textrm{N}(0, \sigma_{u}^{2}),

u^{*} \sim \textrm{N}(0, \sigma_{u}^{2}),

\tau^{2} \sim \textrm{Inverse-Gamma}(\alpha_{1}, \xi_{1}),

\rho \sim \textrm{Uniform}(0, 1),

\sigma_{u}^{2} \sim \textrm{Inverse-Gamma}(\alpha_{2}, \xi_{2}),

\sigma_{e}^{2} \sim \textrm{Inverse-Gamma}(\alpha_{3}, \xi_{3}).

The covariates for the ith individual in the sth spatial unit or other grouping are included in a p \times 1 vector \boldsymbol{x}_{i_s}. The corresponding p \times 1 vector of fixed effect parameters are denoted by \boldsymbol{\beta}, which has an assumed multivariate Gaussian prior with mean \boldsymbol{0} and diagonal covariance matrix \alpha\boldsymbol{I} that can be chosen by the user. A conjugate Inverse-Gamma prior is specified for \sigma_{e}^{2}, and the corresponding hyperparamaterers (\alpha_{3}, \xi_{3}) can be chosen by the user.

Spatial correlation in these areal unit level random effects is most often modelled by a conditional autoregressive (CAR) prior distribution. Using this model spatial correlation is induced into the random effects via a non-negative spatial adjacency matrix \boldsymbol{A} = (a_{sl})_{S \times S}, which defines how spatially close the S areal units are to each other. The elements of \boldsymbol{A}_{S \times S} can be binary or non-binary, and the most common specification is that a_{sl} = 1 if a pair of areal units (\mathcal{G}_{s}, \mathcal{G}_{l}) share a common border or are considered neighbours by some other measure, and a_{sl} = 0 otherwise. Note, a_{ss} = 0 for all s. \boldsymbol{\phi}_{-s}=(\phi_1,\ldots,\phi_{s-1}, \phi_{s+1},\ldots,\phi_{S}). Here \tau^{2} is a measure of the variance relating to the spatial random effects \boldsymbol{\phi}, while \rho controls the level of spatial autocorrelation, with values close to one and zero representing strong autocorrelation and independence respectively. A non-conjugate uniform prior on the unit interval is specified for the single level of spatial autocorrelation \rho. In contrast, a conjugate Inverse-Gamma prior is specified for the random effects variance \tau^{2}, and corresponding hyperparamaterers (\alpha_{1}, \xi_{1}) can be chosen by the user.

The J \times 1 vector of alter random effects are denoted by \boldsymbol{u} = (u_{1}, \ldots, u_{J})_{J \times 1} and modelled as independently Gaussian with mean zero and a constant variance, and due to the row standardised nature of \boldsymbol{W}, \sum_{j \in \textrm{net}(i_s)} w_{i_sj} u_{j} represents the average (mean) effect that the peers of individual i in spatial unit or group s have on that individual. w^{*}_{i_s}u^{*} is an isolation effect, which is an effect for individuals who don't nominate any friends. This is achieved by setting w^{*}_{i_s}=1 if individual i_s nominates no peers and w^{*}_{i_s}=0 otherwise, and if w^{*}_{i_s}=1 then clearly \sum_{j\in \textrm{net}(i_{s})}w_{i_{s}j}u_{jr}=0 as \textrm{net}(i_{s}) is the empty set. A conjugate Inverse-Gamma prior is specified for the random effects variance \sigma_{u}^{2}, and the corresponding hyperparamaterers (\alpha_{2}, \xi_{2}) can be chosen by the user.

The exact specification of each of the likelihoods (binomial, Gaussian, and Poisson) are given below:

\textrm{Binomial:} ~ Y_{i_s} \sim \textrm{Binomial}(n_{i_s}, \theta_{i_s}) ~ \textrm{and} ~ g(\mu_{i_s}) = \textrm{ln}(\theta_{i_s} / (1 - \theta_{i_s})),

\textrm{Gaussian:} ~ Y_{i_s} \sim \textrm{N}(\mu_{i_s}, \sigma_{e}^{2}) ~ \textrm{and} ~ g(\mu_{i_s}) = \mu_{i_s},

\textrm{Poisson:} ~ Y_{i_s} \sim \textrm{Poisson}(\mu_{i_s}) ~ \textrm{and} ~ g(\mu_{i_s}) = \textrm{ln}(\mu_{i_s}).

Usage

uniNetLeroux(formula, data, trials, family,
squareSpatialNeighbourhoodMatrix, spatialAssignment, W, numberOfSamples = 10, 
burnin = 0, thin = 1, seed = 1, trueBeta = NULL, 
trueSpatialRandomEffects = NULL, trueURandomEffects = NULL, 
trueSpatialTauSquared = NULL, trueSpatialRho = NULL, trueSigmaSquaredU = NULL,
trueSigmaSquaredE = NULL, covarianceBetaPrior = 10^5, a1 = 0.001, b1 = 0.001, 
a2 = 0.001, b2 = 0.001, a3 = 0.001, b3 = 0.001, 
centerSpatialRandomEffects = TRUE, centerURandomEffects = TRUE)

Arguments

`formula`	A formula for the covariate part of the model using a similar syntax to that used in the lm() function.
`data`	An optional data.frame containing the variables in the formula.
`trials`	A vector the same length as the response containing the total number of trials `n_{i_s}`. Only used if `\texttt{family}`=“binomial".
`family`	The data likelihood model that must be “gaussian", “poisson" or “binomial".
`squareSpatialNeighbourhoodMatrix`	An `S \times S` symmetric and non-negative neighbourhood matrix `\boldsymbol{A} = (a_{sl})_{S \times S}`.
`W`	A matrix `\boldsymbol{W}` that encodes the social network structure and whose rows sum to 1.
`spatialAssignment`	The binary matrix of individual's assignment to spatial area used in the model fitting process.
`numberOfSamples`	The number of samples to generate pre-thin.
`burnin`	The number of MCMC samples to discard as the burn-in period.
`thin`	The value by which to thin `\texttt{numberOfSamples}`.
`seed`	A seed for the MCMC algorithm.
`trueBeta`	If available, the true value of `\boldsymbol{\beta}`.
`trueSpatialRandomEffects`	If available, the true value of `\boldsymbol{\phi}`.
`trueURandomEffects`	If available, the true value of `\boldsymbol{u}`.
`trueSpatialTauSquared`	If available, the true value of `\tau^{2}`.
`trueSpatialRho`	If available, the true value of`\rho`.
`trueSigmaSquaredU`	If available, the true value of `\sigma_{u}^{2}`.
`trueSigmaSquaredE`	If available, the true value of `\sigma_{e}^{2}`.
`covarianceBetaPrior`	A scalar prior `\alpha` for the covariance parameter of the beta prior, such that the covariance is `\alpha\boldsymbol{I}`.
`a1`	The shape parameter for the Inverse-Gamma distribution relating to the spatial random effects `\alpha_{1}`.
`b1`	The scale parameter for the Inverse-Gamma distribution relating to the spatial random effects `\xi_{1}`.
`a2`	The shape parameter for the Inverse-Gamma distribution relating to the network random effects `\alpha_{2}`.
`b2`	The scale parameter for the Inverse-Gamma distribution relating to the network random effects `\xi_{2}`.
`a3`	The shape parameter for the Inverse-Gamma distribution relating to the error terms `\alpha_{3}`. Only used if `\texttt{family}`=“gaussian".
`b3`	The scale parameter for the Inverse-Gamma distribution relating to the error terms `\xi_{3}`. Only used if `\texttt{family}`=“gaussian".
`centerSpatialRandomEffects`	A choice to center the spatial random effects after each iteration of the MCMC sampler.
`centerURandomEffects`	A choice to center the network random effects after each iteration of the MCMC sampler.

Value

`call`	The matched call.
`y`	The response used.
`X`	The design matrix used.
`standardizedX`	The standardized design matrix used.
`squareSpatialNeighbourhoodMatrix`	The spatial neighbourhood matrix used.
`spatialAssignment`	The spatial assignment matrix used.
`W`	The network matrix used.
`samples`	The matrix of simulated samples from the posterior distribution of each parameter in the model (excluding random effects).
`betaSamples`	The matrix of simulated samples from the posterior distribution of `\boldsymbol{\beta}` parameters in the model.
`spatialTauSquaredSamples`	The vector of simulated samples from the posterior distribution of `\tau^{2}` in the model.
`spatialRhoSamples`	The vector of simulated samples from the posterior distribution of `\rho` in the model.
`sigmaSquaredUSamples`	The vector of simulated samples from the posterior distribution of `\sigma_{u}^{2}` in the model.
`sigmaSquaredESamples`	The vector of simulated samples from the posterior distribution of `\sigma_{e}^{2}` in the model.
`spatialRandomEffectsSamples`	The matrix of simulated samples from the posterior distribution of spatial/grouping random effects `\boldsymbol{\phi}` in the model.
`uRandomEffectsSamples`	The matrix of simulated samples from the posterior distribution of network random effects `\boldsymbol{u}` in the model.
`acceptanceRates`	The acceptance rates of parameters in the model (excluding random effects) from the MCMC sampling scheme .
`spatialRandomEffectsAcceptanceRate`	The acceptance rates of spatial/grouping random effects in the model from the MCMC sampling scheme.
`uRandomEffectsAcceptanceRate`	The acceptance rates of network random effects in the model from the MCMC sampling scheme.
`timeTaken`	The time taken for the model to run.
`burnin`	The number of MCMC samples to discard as the burn-in period.
`thin`	The value by which to thin `\texttt{numberOfSamples}`.
`DBar`	DBar for the model.
`posteriorDeviance`	The posterior deviance for the model.
`posteriorLogLikelihood`	The posterior log likelihood for the model.
`pd`	The number of effective parameters in the model.
`DIC`	The DIC for the model.

Author(s)

George Gerogiannis

Examples

  #################################################
  #### Run the model on simulated data
  #################################################
  #### Load other libraries required
  library(MCMCpack)
  
  #### Set up a network
  observations <- 200
  numberOfMultipleClassifications <- 50
  W <- matrix(rbinom(observations * numberOfMultipleClassifications, 1, 0.05), 
              ncol = numberOfMultipleClassifications)
  numberOfActorsWithNoPeers <- sum(apply(W, 1, function(x) { sum(x) == 0 }))
  peers <- sample(1:numberOfMultipleClassifications, numberOfActorsWithNoPeers,
  TRUE)
  actorsWithNoPeers <- which(apply(W, 1, function(x) { sum(x) == 0 }))
  for(i in 1:numberOfActorsWithNoPeers) {
    W[actorsWithNoPeers[i], peers[i]] <- 1
  }
  W <- t(apply(W, 1, function(x) { x / sum(x) }))
  
  #### Set up a spatial structure
  numberOfSpatialAreas <- 100
  factor = sample(1:numberOfSpatialAreas, observations, TRUE)
  spatialAssignment = matrix(NA, ncol = numberOfSpatialAreas, 
                             nrow = observations)
  for(i in 1:length(factor)){
    for(j in 1:numberOfSpatialAreas){
      if(factor[i] == j){
        spatialAssignment[i, j] = 1
      } else {
        spatialAssignment[i, j] = 0
      }
    }
  }
  
  gridAxis = sqrt(numberOfSpatialAreas)
  easting = 1:gridAxis
  northing = 1:gridAxis
  grid = expand.grid(easting, northing)
  numberOfRowsInGrid = nrow(grid)
  distance = as.matrix(dist(grid))
  squareSpatialNeighbourhoodMatrix = array(0, c(numberOfRowsInGrid, 
                                                numberOfRowsInGrid))
  squareSpatialNeighbourhoodMatrix[distance==1] = 1

  #### Generate the covariates and response data
  X <- matrix(rnorm(2 * observations), ncol = 2)
  colnames(X) <- c("x1", "x2")
  beta <- c(2, -2, 2)
  
  spatialRho <- 0.5
  spatialTauSquared <- 2
  spatialPrecisionMatrix = spatialRho * 
    (diag(apply(squareSpatialNeighbourhoodMatrix, 1, sum)) -
     squareSpatialNeighbourhoodMatrix) + (1 - spatialRho) * 
     diag(rep(1, numberOfSpatialAreas))
  spatialCovarianceMatrix = solve(spatialPrecisionMatrix)
  spatialPhi = mvrnorm(n = 1, mu = rep(0, numberOfSpatialAreas), 
                       Sigma = (spatialTauSquared * spatialCovarianceMatrix))
  
  sigmaSquaredU <- 2
  uRandomEffects <- rnorm(numberOfMultipleClassifications, mean = 0, 
                          sd = sqrt(sigmaSquaredU))
  
  logit <- cbind(rep(1, observations), X) %*% beta + 
    spatialAssignment %*% spatialPhi + W %*% uRandomEffects
  prob <- exp(logit) / (1 + exp(logit))
  trials <- rep(50, observations)
  Y <- rbinom(n = observations, size = trials, prob = prob)
  data <- data.frame(cbind(Y, X))
  
  #### Run the model
  formula <- Y ~ x1 + x2
  ## Not run: model <- uniNetLeroux(formula = formula, data = data, 
    family="binomial",  W = W,
    spatialAssignment = spatialAssignment, 
    squareSpatialNeighbourhoodMatrix = squareSpatialNeighbourhoodMatrix,
    trials = trials, numberOfSamples = 10000, 
    burnin = 10000, thin = 10, seed = 1)
## End(Not run)

[Package netcmc version 1.0.2 Index]