uniNetLeroux: 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 $i$th individual in the $s$th 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)

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.

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.

Author

George Gerogiannis

Examples

Run this code

  #################################################
  #### 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
  if (FALSE) model <- uniNetLeroux(formula = formula, data = data, 
    family="binomial",  W = W,
    spatialAssignment = spatialAssignment, 
    squareSpatialNeighbourhoodMatrix = squareSpatialNeighbourhoodMatrix,
    trials = trials, numberOfSamples = 10000, 
    burnin = 10000, thin = 10, seed = 1)

Run the code above in your browser using DataLab