varian {varian}

Variablity Analysis using a Bayesian Variability Model (VM)

Description

This function uses a linear mixed effects model that assumes the level 1 residual variance varies by Level 2 units. That is rather than assuming a homogenous residual variance, it assumes the residual standard deviations come from a Gamma distribution. In the first stage of this model, each Level 2's residual standard deviation is estimated, and in the second stage, these standard deviations are used to predict another Level 2 outcome. The interface uses an intuitive formula interface, but the underlying model is implemented in Stan, with minimally informative priors for all parameters.

The Variability Analysis in R Package

Usage

varian(y.formula, v.formula, m.formula, data, design = c("V -> Y",
  "V -> M -> Y", "V", "X -> V", "X -> V -> Y", "X -> M -> V"), useU = TRUE,
  totaliter = 2000, warmup = 1000, chains = 1, inits = NULL, modelfit,
  opts = list(SD_Tol = 0.01, pars = NULL), ...)

Arguments

y.formula

A formula describing a model for the outcome. At present, this must be a continuous, normally distributed variable.

v.formula

A formula describing a model for the variability. Note this must end with | ID, where ID is the name of the ID variable in the dataset. At present, this must be a continuous, normally distributed variable.

m.formula

An optional formula decribing a model for a mediatior variable. At present, this must be a continuous normally distributed variable.

data

A long data frame containing an both the Level 2 and Level 1 outcomes, as well as all covariates and an ID variable.

design

A character string indicating the type of model to be run. One of “V -> Y” for variability predicting an outcome, “V -> M -> Y” for mediation of variability on an outcome, “V” to take posterior samples of individual variability estimates alone.

useU

A logical value whether the latent intercept estimated in Stage 1 should also be used as a predictor. Defaults to TRUE. Note if there is a mediator as well as main outcome, the latent intercepts will be used as a predictor for both.

totaliter

The total number of iterations to be used (not including the warmup iterations), these are distributed equally across multiple independent chains.

warmup

The number of warmup iterations. Each independent chain has the same number of warmup iterations, before it starts the iterations that will be used for inference.

chains

The number of independent chains to run (default to 1).

inits

Initial values passed on to stan. If NULL, the default, initial values are estimated means, standard deviations, and coefficients from a single level linear regression.

modelfit

A compiled Stan model (e.g., from a previous run).

opts

A list giving options. Currently only SD_Tol which controls the tolerance for how small a variables standard deviation may be without stopping estimation (this ensures that duplicate variables, or variables without any variability are included as predictors).

...

Additional arguments passed to stan.

Value

A named list containing the model results, the model, the variable.names, the data, the random seeds, and the initial function .call.

Author(s)

Joshua F. Wiley <[e-mail]>

Examples

## Not run: 
  sim.data <- with(simulate_gvm(4, 60, 0, 1, 3, 2, 94367), {
    set.seed(265393)
    x2 <- MASS::mvrnorm(k, c(0, 0), matrix(c(1, .3, .3, 1), 2))
    y2 <- rnorm(k, cbind(Int = 1, x2) %*% matrix(c(3, .5, .7)) + sigma, sd = 3)
    data.frame(
      y = Data$y,
      y2 = y2[Data$ID2],
      x1 = x2[Data$ID2, 1],
      x2 = x2[Data$ID2, 2],
      ID = Data$ID2)
  })
  m <- varian(y2 ~ x1 + x2, y ~ 1 | ID, data = sim.data, design = "V -> Y",
    totaliter = 10000, warmup = 1500, thin = 10, chains = 4, verbose=TRUE)

  # check diagnostics
  vm_diagnostics(m)

  sim.data2 <- with(simulate_gvm(21, 250, 0, 1, 3, 2, 94367), {
    set.seed(265393)
    x2 <- MASS::mvrnorm(k, c(0, 0), matrix(c(1, .3, .3, 1), 2))
    y2 <- rnorm(k, cbind(Int = 1, x2) %*% matrix(c(3, .5, .7)) + sigma, sd = 3)
    data.frame(
      y = Data$y,
      y2 = y2[Data$ID2],
      x1 = x2[Data$ID2, 1],
      x2 = x2[Data$ID2, 2],
      ID = Data$ID2)
  })
  # warning: may take several minutes
  m2 <- varian(y2 ~ x1 + x2, y ~ 1 | ID, data = sim.data2, design = "V -> Y",
    totaliter = 10000, warmup = 1500, thin = 10, chains = 4, verbose=TRUE)
  # check diagnostics
  vm_diagnostics(m2)

## End(Not run)