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

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 <josh@elkhartgroup.com>

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)

Variability Measures

Description

Variability Measures

by_id - Internal function to allow a simple statistic (e.g., SD) to be calculated individually by an ID variable and returned either as per ID (i.e., wide form) or for every observation of an ID (i.e., long form).

sd_id - Calculates the standard deviation of observations by ID.

rmssd - Calculates the root mean square of successive differences (RMSSD). Note that missing values are removed.

rmssd_id - Calculates the RMSSD by ID.

rolling_diff - Calculates the average rolling difference of the data. Within each window, the difference between the maximum and minimum value is computed and these are averaged across all windows. The equation is:

\frac{\sum_{t = 1}^{N - k} max(x_{t}, \ldots, x_{t + k}) - min(x_{t}, \ldots, x_{t + k})}{N - k}

rolling_diff_id - Calculates the average rolling difference by ID

Usage

by_id(x, ID, fun, long = TRUE, ...)

sd_id(x, ID, long = TRUE)

rmssd(x)

rmssd_id(x, ID, long = TRUE)

rolling_diff(x, window = 4)

rolling_diff_id(x, ID, long = TRUE, window = 4)

Arguments

Value

by_id - A vector the same length as x if long=TRUE, or the length of unique IDs if long=FALSE.

sd_id - A vector of the standard deviations by ID

rmssd - The RMSSD for the data.

rmssd_id - A vector of the RMSSDs by ID

rolling_diff - The average of the rolling differences between maximum and minimum.

rolling_diff_id - A vector of the average rolling differences by ID

Note

These are a set of functions designed to calculate various measures of variability either on a single data vector, or calculate them by an ID.

Author(s)

Joshua F. Wiley <josh@elkhartgroup.com>

Examples

sd_id(mtcars$mpg, mtcars$cyl, long=TRUE)
sd_id(mtcars$mpg, mtcars$cyl, long=FALSE)
rmssd(1:4)
rmssd(c(1, 3, 2, 4))
rmssd_id(mtcars$mpg, mtcars$cyl)
rmssd_id(mtcars$mpg, mtcars$cyl, long=FALSE)
rolling_diff(1:7, window = 4)
rolling_diff(c(1, 4, 3, 4, 5))
rolling_diff_id(mtcars$mpg, mtcars$cyl, window = 3)

Calculates an empirical p-value based on the data

Description

This function takes a vector of statistics and calculates the empirical p-value, that is, how many fall on the other side of zero. It calculates a two-tailed p-value.

Usage

empirical_pvalue(x, na.rm = TRUE)

Arguments

Value

a named vector with the number of values falling at or below zero, above zero, and the empirical p-value.

Author(s)

Joshua F. Wiley <josh@elkhartgroup.com>

Examples

empirical_pvalue(rnorm(100))

Estimate the parameters for a Gamma distribution

Description

This is a simple function to estimate what the parameters for a Gamma distribution would be from a data vector. It is used internally to generate start values.

Usage

gamma_params(x)

Arguments

Value

a list of the shape (alpha) and rate (beta) parameters and the mean and variance

Author(s)

Joshua F. Wiley <josh@elkhartgroup.com>

Wrapper for the stan function to parallelize chains

Description

This funcntion takes Stan model code, compiles the Stan model, and then runs multiple chains in parallel.

Usage

parallel_stan(model_code, standata, totaliter, warmup, thin = 1, chains, cl,
  cores, seeds, modelfit, verbose = FALSE, pars = NA, sample_file = NA,
  diagnostic_file = NA, init = "random", ...)

Arguments

Value

a named list with three elements, the results, compiled Stan model, and the random seeds

Author(s)

Joshua F. Wiley <josh@elkhartgroup.com>

Examples

# Make me!

Calculates summaries for a parameter

Description

This function takes a vector of statistics and calculates several summaries: mean, median, 95 the empirical p-value, that is, how many fall on the other side of zero.

Usage

param_summary(x, digits = 2, pretty = FALSE, ..., na.rm = TRUE)

Arguments

Value

.

Author(s)

Joshua F. Wiley <josh@elkhartgroup.com>

Examples

param_summary(rnorm(100))
param_summary(rnorm(100), pretty = TRUE)

nice formatting for p-values

Description

nice formatting for p-values

Usage

pval_smartformat(p, d = 3, sd = 5)

Arguments

Author(s)

Joshua F. Wiley <josh@elkhartgroup.com>

Examples

varian:::pval_smartformat(c(1, .15346, .085463, .05673, .04837, .015353462,
  .0089, .00164, .0006589, .0000000053326), 3, 5)

Estimates the parameters of a Gamma distribution from SDs

Description

This function calcualtes the parameters of a Gamma distribution from the residuals from an individuals' own mean. That is, the distribution of (standard) deviations from individuals' own mean are calculated and then an estimate of the parameters of a Gamma distribution are calculated.

Usage

res_gamma(x, ID)

Arguments

Value

a list of the shape (alpha) and rate (beta) parameters and the mean and variance

Author(s)

Joshua F. Wiley <josh@elkhartgroup.com>

Examples

set.seed(1234)
y <- rgamma(100, 3, 2)
x <- rnorm(100 * 10, mean = 0, sd = rep(y, each = 10))
ID <- rep(1:100, each = 10)
res_gamma(x, ID)

Simulate a Gamma Variability Model

Description

This function facilitates simulation of a Gamma Variability Model and allows the number of units and repeated measures to be varied as well as the degree of variability.

Usage

simulate_gvm(n, k, mu, mu.sigma, sigma.shape, sigma.rate, seed = 5346)

Arguments

Value

a list of the data, IDs, and the parameters used for the simulation

Author(s)

Joshua F. Wiley <josh@elkhartgroup.com>

Examples

raw.sim <- simulate_gvm(12, 140, 0, 1, 4, .1, 94367)
sim.data <- with(raw.sim, {
  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)
})

Calculate Initial Values for Stan VM Model

Description

Internal function used to get rough starting values for a variability model in Stan. Uses inidivudal standard deviations, means, and linear regressions.

Usage

stan_inits(stan.data, design = c("V -> Y", "V -> M -> Y", "V", "X -> V",
  "X -> V -> Y", "X -> M -> V"), useU, ...)

Arguments

Value

A named list containing the initial values for Stan.

Author(s)

Joshua F. Wiley <josh@elkhartgroup.com>

Examples

# make me!

Plot diagnostics from a VM model

Description

This function plots a variety of diagnostics from a Variability Model. These include a histogram of the Rhat values (so-called percent scale reduction factors). An Rhat value of 1 indicates that no reduction in the variability of the estimates is possible from running the chain longer. Values below 1.10 or 1.05 are typically considered indicative of convergence, with higher values indicating the model did not converge and should be changed or run longer. A histogram of the effective sample size indicates for every parameter estimated how many effective posterior samples are available for inference. Low values may indicate high autocorrelation in the samples and may be a sign of failure to converge. The maximum possible will be the total iterations available. Histograms of the posterior medians for the latent variability and intercept estimates are also shown.

Usage

vm_diagnostics(object, plot = TRUE, ...)

Arguments

Value

A graphical object

Author(s)

Joshua F. Wiley <josh@elkhartgroup.com>

Examples

# Make Me!

Create a Stan class VM object

Description

Internal function to create and compile a Stan model.

Usage

vm_stan(design = c("V -> Y", "V -> M -> Y", "V", "X -> V", "X -> V -> Y",
  "X -> M -> V"), useU = TRUE, ...)

Arguments

Value

A compiled Stan model.

Author(s)

Joshua F. Wiley <josh@elkhartgroup.com>

Examples

# Make Me!
## Not run: 
  test1 <- vm_stan("V -> Y", useU=TRUE)
  test2 <- vm_stan("V -> Y", useU=FALSE)
  test3 <- vm_stan("V -> M -> Y", useU=TRUE)
  test4 <- vm_stan("V -> M -> Y", useU=FALSE)
  test5 <- vm_stan("V")

## End(Not run)

Plot the posterior distributions of the focal parameters from a VM model

Description

This function plots the univariate and bivariate (if applicable) distributions of the focal (alpha) parameters from a Variability Model where the variability is used as a predictor in a second-stage model. The latent variability estimates are referred to as “Sigma” and, if used, the latent intercepts are referred to as “U”.

Usage

vmp_plot(alpha, useU = TRUE, plot = TRUE, digits = 3, ...)

Arguments

Value

A list containing the Combined and the Individual plot objects.

Author(s)

Joshua F. Wiley <josh@elkhartgroup.com>

Examples

# Using made up data because the real models take a long time to run
set.seed(1234) # make reproducible
vmp_plot(matrix(rnorm(1000), ncol = 2))