Learn R Programming

netcmc (version 1.0.2)

uni: A function that generates samples for a univariate fixed effects model.

Description

This function generates samples for a univariate fixed effects 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},$$ $$\boldsymbol{\beta} \sim \textrm{N}(\boldsymbol{0}, \alpha\boldsymbol{I}),$$ $$\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.

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

uni(formula, data, trials, family, numberOfSamples = 10, burnin = 0, thin = 1, seed = 1,
trueBeta = NULL, trueSigmaSquaredE = NULL, covarianceBetaPrior = 10^5, 
a3 = 0.001, b3 = 0.001)

Value

call

The matched call.

y

The response used.

X

The design matrix used.

standardizedX

The standardized design 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.

sigmaSquaredESamples

The vector of simulated samples from the posterior distribution of \(\sigma_{e}^{2}\) in the model.

acceptanceRates

The acceptance rates of parameters 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".

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 values of the \(\boldsymbol{\beta}\).

trueSigmaSquaredE

If available, the true value of \(\sigma_{e}^{2}\). Only used if \(\texttt{family}\)=``gaussian".

covarianceBetaPrior

A scalar prior \(\alpha\) for the covariance parameter of the beta prior, such that the covariance is \(\alpha\boldsymbol{I}\).

a3

The shape parameter for the Inverse-Gamma distribution \(\alpha_{3}\). Only used if \(\texttt{family}\)=``gaussian".

b3

The scale parameter for the Inverse-Gamma distribution \(\xi_{3}\). Only used if \(\texttt{family}\)=``gaussian".

Author

George Gerogiannis

Examples

Run this code
  #################################################
  #### Run the model on simulated data
  #################################################
  
  #### Generate the covariates and response data
  observations <- 100
  X <- matrix(rnorm(2 * observations), ncol = 2)
  colnames(X) <- c("x1", "x2")
  beta <- c(2, -2, 2)
  logit <- cbind(rep(1, observations), X) %*% beta
  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 <- uni(formula = formula, data = data, family="binomial", 
                        trials = trials, numberOfSamples = 10000, 
                        burnin = 10000, thin = 10, seed = 1)
                        

Run the code above in your browser using DataLab