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

In netcmc: Spatio-Network Generalised Linear Mixed Models for Areal Unit and Network Data

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}|μ_{i_s} \sim f(y_{i_s}| μ_{i_s}, σ_{e}^{2}) ~~~ i=1,…, N_{s},~s=1,…,S ,

g(μ_{i_s}) = \boldsymbol{x}^\top_{i_s} \boldsymbol{β} + φ_{s} + ∑_{j\in \textrm{net}(i_s)}w_{i_sj}u_{j} + w^{*}_{i_s}u^{*},

\boldsymbol{β} \sim \textrm{N}(\boldsymbol{0}, α\boldsymbol{I}),

φ_{s} | \boldsymbol{φ}_{-s} \sim \textrm{N}\bigg(\frac{ ρ ∑_{l = 1}^{S} a_{sl} φ_{l} }{ ρ ∑_{l = 1}^{S} a_{sl} + 1 - ρ }, \frac{ τ^{2} }{ ρ ∑_{l = 1}^{S} a_{sl} + 1 - ρ } \bigg),

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

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

τ^{2} \sim \textrm{Inverse-Gamma}(α_{1}, ξ_{1}),

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

σ_{u}^{2} \sim \textrm{Inverse-Gamma}(α_{2}, ξ_{2}),

σ_{e}^{2} \sim \textrm{Inverse-Gamma}(α_{3}, ξ_{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{β}, which has an assumed multivariate Gaussian prior with mean \boldsymbol{0} and diagonal covariance matrix α\boldsymbol{I} that can be chosen by the user. A conjugate Inverse-Gamma prior is specified for σ_{e}^{2}, and the corresponding hyperparamaterers (α_{3}, ξ_{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{φ}_{-s}=(φ_1,…,φ_{s-1}, φ_{s+1},…,φ_{S}). Here τ^{2} is a measure of the variance relating to the spatial random effects \boldsymbol{φ}, while ρ 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 ρ. In contrast, a conjugate Inverse-Gamma prior is specified for the random effects variance τ^{2}, and corresponding hyperparamaterers (α_{1}, ξ_{1}) can be chosen by the user.

The J \times 1 vector of alter random effects are denoted by \boldsymbol{u} = (u_{1}, …, 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}, ∑_{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 ∑_{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 σ_{u}^{2}, and the corresponding hyperparamaterers (α_{2}, ξ_{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}, θ_{i_s}) ~ \textrm{and} ~ g(μ_{i_s}) = \textrm{ln}(θ_{i_s} / (1 - θ_{i_s})),

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

\textrm{Poisson:} ~ Y_{i_s} \sim \textrm{Poisson}(μ_{i_s}) ~ \textrm{and} ~ g(μ_{i_s}) = \textrm{ln}(μ_{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{β}.
trueSpatialRandomEffects	If available, the true value of \boldsymbol{φ}.
trueURandomEffects	If available, the true value of \boldsymbol{u}.
trueSpatialTauSquared	If available, the true value of τ^{2}.
trueSpatialRho	If available, the true value ofρ.
trueSigmaSquaredU	If available, the true value of σ_{u}^{2}.
trueSigmaSquaredE	If available, the true value of σ_{e}^{2}.
covarianceBetaPrior	A scalar prior α for the covariance parameter of the beta prior, such that the covariance is α\boldsymbol{I}.
a1	The shape parameter for the Inverse-Gamma distribution relating to the spatial random effects α_{1}.
b1	The scale parameter for the Inverse-Gamma distribution relating to the spatial random effects ξ_{1}.
a2	The shape parameter for the Inverse-Gamma distribution relating to the network random effects α_{2}.
b2	The scale parameter for the Inverse-Gamma distribution relating to the network random effects ξ_{2}.
a3	The shape parameter for the Inverse-Gamma distribution relating to the error terms α_{3}. Only used if \texttt{family}=“gaussian".
b3	The scale parameter for the Inverse-Gamma distribution relating to the error terms ξ_{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{β} parameters in the model.
spatialTauSquaredSamples	The vector of simulated samples from the posterior distribution of τ^{2} in the model.
spatialRhoSamples	The vector of simulated samples from the posterior distribution of ρ in the model.
sigmaSquaredUSamples	The vector of simulated samples from the posterior distribution of σ_{u}^{2} in the model.
sigmaSquaredESamples	The vector of simulated samples from the posterior distribution of σ_{e}^{2} in the model.
spatialRandomEffectsSamples	The matrix of simulated samples from the posterior distribution of spatial/grouping random effects \boldsymbol{φ} 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)

netcmc documentation built on Nov. 10, 2022, 5:11 p.m.