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

uniNetRand {netcmc}

R Documentation

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

Description

This function generates samples for a univariate network group 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} + v_{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}),

v_{s} \sim \textrm{N}(0, \tau^{2}),

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}),

\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.

The S \times 1 vector of random effects for the groups are collectively denoted by \boldsymbol{v} = (v_{1}, \ldots, v_{S})_{S \times 1}, and each element is assigned an independent zero-mean Gaussian prior distribution with a constant variance \tau^{2}. A conjugate Inverse-Gamma prior is specified for \tau^{2}. The 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

uniNetRand(formula, data, trials, family, groupAssignment, W, numberOfSamples = 10, 
burnin = 0, thin = 1, seed = 1, trueBeta = NULL, 
trueGroupRandomEffects = NULL, trueURandomEffects = NULL, 
trueTauSquared = 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, 
centerGroupRandomEffects = 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".
`W`	A matrix `\boldsymbol{W}` that encodes the social network structure and whose rows sum to 1.
`groupAssignment`	The binary matrix of individual's assignment to groups 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}`.
`trueGroupRandomEffects`	If available, the true value of `\boldsymbol{v}`.
`trueURandomEffects`	If available, the true value of `\boldsymbol{u}`.
`trueTauSquared`	If available, the true value `\tau^{2}`.
`trueSigmaSquaredU`	If available, the true value `\sigma_{u}^{2}`.
`trueSigmaSquaredE`	If available, the true value `\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 group random effects `\alpha_{1}`.
`b1`	The shape parameter for the Inverse-Gamma distribution relating to the group 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".
`centerGroupRandomEffects`	A choice to center the group 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.
`groupAssignment`	The group 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.
`tauSquaredSamples`	The vector of simulated samples from the posterior distribution of `\tau^{2}` 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.
`groupRandomEffectsSamples`	The matrix of simulated samples from the posterior distribution of spatial/grouping random effects `\boldsymbol{v}` 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 .
`groupRandomEffectsAcceptanceRate`	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 single level classification
  numberOfSingleClassifications <- 20
  factor = sample(1:numberOfSingleClassifications, observations, TRUE)
  V = matrix(NA, ncol = numberOfSingleClassifications, nrow = observations)
  for(i in 1:length(factor)){
    for(j in 1:numberOfSingleClassifications){
      if(factor[i] == j){
        V[i, j] = 1
      } else {
        V[i, j] = 0
      }
    }
  }
  
  #### Generate the covariates and response data
  X <- matrix(rnorm(2 * observations), ncol = 2)
  colnames(X) <- c("x1", "x2")
  beta <- c(1, -0.5, 0.5)
  tauSquared <- 0.5
  vRandomEffects <- rnorm(numberOfSingleClassifications, mean = 0, 
                          sd = sqrt(tauSquared))
  sigmaSquaredU <- 1
  uRandomEffects <- rnorm(numberOfMultipleClassifications, mean = 0, 
                          sd = sqrt(sigmaSquaredU))

  logTheta <- cbind(rep(1, observations), X) %*% beta + V %*% vRandomEffects 
            + W %*% uRandomEffects
  Y <- rpois(n = observations, lambda = exp(logTheta))
  data <- data.frame(cbind(Y, X))
  
  #### Run the model
  formula <- Y ~ x1 + x2
  ## Not run: model <- uniNetRand(formula = formula, data = data, family="poisson", 
                               W = W, groupAssignment = V, 
                               numberOfSamples = 10000, burnin = 10000,
                               thin = 10, seed = 1)
## End(Not run)

[Package netcmc version 1.0.2 Index]