In donaldRwilliams/hypMuVar: Bayesian Hypothesis Testing of Mean-Variance Relations

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  dev = "png",
  dpi = 500,
  fig.align = "center"
)

hypMuVar: Bayesian Hypothesis
Testing of Mean-Variance Relations

Overview

This is a dedicated R package for testing mean-variance relations. Currently the methodology is for Bayesian mixed-effects location scale models [MELSM, mel$\cdot$zem; @Hedeker2008], where mean and (within-person) variance relations are tested in the distribution of random effects. This is accomplished with covariance selection, in which the random effects correlations are tested with a mixture of prior distributions. In future releases, this package will allow for testing other mean-variance relations, including linear models (no random effects) and fixed effects in the MELSM.

Who actually studies mean-variance relations ?

Across the sciences, there is large literature investigating mean-variance relations [e.g., @nakagawa2015meta; @xiao2015process; @moreno2003relationship]. The typical pattern is that means are positively correlated with variances (or standard deviations)^[Interestingly, I recently came across a paper investigating the mean-variance relation in guppies [@mitchell2016towards]. Apparently, in that field, the MELSM has been termed the doubly hierarchical model (see here dhglm).]. As one example, in ecology the relationship between the mean and variance of species abundance is known as Taylor's law. That Wikipedia page in particular includes interesting references that date back to the 1920's (e.g., Jerzy Neyman studied the mean-variance relation in 1926). Furthermore, a related line of research has set out to identify processes that violate the "law" [e.g., @Pratte2010; @Schwarz2012].

In my field (psychology), there is also a rich (perhaps the richest!) tradition of studying mean-variance relations. The most prominent example comes from cognitive psychology, where a lawful, linear relation between reaction time means and standard deviations was proposed in @wagenmakers2007linear. Accordingly, we demonstrated the utility of our methodology in cognitive inhibition tasks [see @williams2019beneath].

Is hypMuVar necessary ?

hypMuVar is not meant to be a general package for mixed-effects modeling (see here brms). There is certainly a need for hypMuVar, however, because the testing strategy is not possible with current software. I anticipate this package will be most useful for researchers specifically interested in testing (possibly confirmatory) hypotheses about the mean-variance relation or the (within-person) variance structure in general (in future releases). This extends inference beyond the customary goal (of testing mean differences) and thus opens the door for answering novel research questions.

Installation

Note you will also have to intall jags (see here download page).

# install.packages("devtools")
devtools::install_github("donaldRwilliams/hypMuVar")

library(hypMuVar)
library(ggplot2)

Random intercepts for the location and scale

This model effectively estimates the reaction time means and variances (on the log scale) for each person. These are hierarchical estimates, which provides shrinkage that can improve accuracy. Note a customary, location only, mixed effects model assumes that each person shares a common variance. This is implemented with

fit_int <- melsm(fixed_location = rt ~ 1, 
                 random_location = ~ 1| ID, 
                 fixed_scale = sigma ~ 1, 
                 random_scale = ~ 1 | ID, k = 3, 
                 mixture = "SSVS", adapt = 5000, 
                 iter = 10000, 
                 rho_test = "muvar",
                 data = stroop)

Note mixture = "SSVS" corresponds to stochastic search variable selection [@George1993]. This is a "spike" and slab approach that utilizes a mixture of prior distributions [see here for reviews @OHara2009; @Malsiner-Walli2011]. In the case of SSVS, there are typically two normal distributions, with the "spike" concentrated around zero. One innovation of hypMuVar is to extend that approach to several mixture components. k = 3 specifies a prior comprised of three distributions that place restrictions to positive values (k = 2), negative values (k = 3), or a null region (k = 1). Setting k = 2 is also an option, as well as the Dirac spike and slab formulation [ @Mitchell1988; @kuo_mallick].

load("C:/Users/willidon/Box Sync/r_packages/rmd_files/fit_int.Rdata")

options(width = 250)
summary(fit_int)

As seen in the summary output, component number 2 has an inclusion probability of 1. This results in an infinite Bayes factor in favor of a positive relation. Furthermore, it is possible to test all random effects correlations by setting rho_test = "all". This would employ the chosen mixture approach for all off-diagonal elements in the random effects covariance matrix.

The individual-specific effects can also be visualized in a variety of ways. For example,

# coefficients
coefs <- coef(fit_int, cred = 0.90)
plot(coefs)[[2]] +
   theme_bw() +
   theme(panel.grid = element_blank(), 
         legend.position = "none") +
   scale_x_discrete(labels = c(1, rep("", 119), 121),
                   expand = c(0.015, 0.015)) 

plots the reaction time standard deviations for each person. These are on the log scale. Note that, in a location only mixed effects model, it is assumed that each person is adequately represented by the dotted line. In these data, there is clearly heterogeneous within-person variance.

Further, it is also possible to visualize the relations between the individual-specific effects. Testing these mean-relations in the distribution of random effects is the primary objective in @williams2019beneath.

cor_plot(fit_int)[[1]]

Random slopes for the location and scale:

fit <- melsm(fixed_location = rt ~  congruency,
             random_location = ~  congruency | ID,
             fixed_scale = sigma ~  congruency,
             random_scale = ~ congruency | ID,  
             mixture = "SSVS", k = 3,
             adapt = 5000, 
             iter = 10000, 
             data = stroop)

load("C:/Users/willidon/Box Sync/r_packages/rmd_files/fit.Rdata")

options(width = 250)
summary(fit)

Marginal Bayes Factors

k2 (positive component) vs. the compliment

bf <- marginal_bf(fit, H1 = "k2")

bf

k2 (positive component) vs. k1(null component)

# note testing against the compliment is preferable
bf <- marginal_bf(fit, H1 = "k2", H2 = "k1")

bf

Posterior Model Probabilities

prob <- model_prob(fit)

prob

Confirmatory Hypothesis Testing

hyp <- list(h1 = c("rho_02 = k2", 
                   "rho_03 = k2", 
                   "rho_12 = k2", 
                   "rho_13 = k2"))

confirm(fit, hyp = hyp)

References

donaldRwilliams/hypMuVar documentation built on Jan. 10, 2020, 9:45 a.m.