This tutorial will focus on the use of Bayesian estimation to fit simple linear regression models …
Software, Statistics, Stan
This tutorial will focus on the use of Bayesian estimation to fit simple linear regression models. BUGS (Bayesian inference Using Gibbs Sampling) is an algorithm and supporting language (resembling R) dedicated to performing the Gibbs sampling implementation of Markov Chain Monte Carlo (MCMC) method. Dialects of the BUGS language are implemented within three main projects:
OpenBUGS - written in component pascal.
JAGS - (Just Another Gibbs Sampler) - written in C++.
Stan - a dedicated Bayesian modelling framework written in C++ and implementing Hamiltonian MCMC samplers.
Whilst the above programs can be used stand-alone, they do offer the rich data pre-processing and graphical capabilities of R, and thus, they are best accessed from within R itself. As such there are multiple packages dedicated to interfacing with the various software implementations:
R2OpenBUGS - interfaces with OpenBUGS
R2jags - interfaces with JAGS
rstan - interfaces with Stan
This tutorial will demonstrate how to fit models in Stan (Gelman, Lee, and Guo (2015)) using the package rstan (Stan Development Team (2018)) as interface, which also requires to load some other packages.
Previous tutorials have concentrated on designs for either continuous (Regression) or categorical (ANOVA) predictor variables. Analysis of covariance (ANCOVA) models are essentially ANOVA models that incorporate one or more continuous and categorical variables (covariates). Although the relationship between a response variable and a covariate may itself be of substantial clinical interest, typically covariate(s) are incorporated to reduce the amount of unexplained variability in the model and thereby increase the power of any treatment effects.
In ANCOVA, a reduction in unexplained variability is achieved by adjusting the response (to each treatment) according to slight differences in the covariate means as well as accounting for any underlying trends between the response and covariate(s). To do so, the extent to which the within treatment group small differences in covariate means between groups and treatment groups are essentially compared via differences in their \(y\)-intercepts. The total variation is thereafter partitioned into explained (using the deviations between the overall trend and trends approximated for each of the treatment groups) and unexplained components (using the deviations between the observations and the approximated within group trends). In this way, ANCOVA can be visualized as a regular ANOVA in which the group and overall means are replaced by group and overall trendlines. Importantly, it should be apparent that ANCOVA is only appropriate when each of the within group trends have the same slope and are thus parallel to one another and the overall trend. Furthermore, ANCOVA is not appropriate when the resulting adjustments must be extrapolated from a linear relationship outside the measured range of the covariate.
As an example, an experiment might be set up to investigate the energetic impacts of sexual vs parthenogenetic (egg development without fertilization) reproduction on leaf insect food consumption. To do so, researchers could measure the daily food intake of individual adult female leaf insects from female only (parthenogenetic) and mixed (sexual) populations. Unfortunately, the available individual leaf insects varied substantially in body size which was expected to increase the variability of daily food intake of treatment groups. Consequently, the researchers also measured the body mass of the individuals as a covariate, thereby providing a means by which daily food consumption could be standardized for body mass. ANCOVA attempts to reduce unexplained variability by standardising the response to the treatment by the effects of the specific covariate condition. Thus ANCOVA provides a means of exercising some statistical control over the variability when it is either not possible or not desirable to exercise experimental control (such as blocking or using otherwise homogeneous observations).
Factor A: the main treatment effect
The adjusted population group means are all equal. The mean of population \(1\) adjusted for the covariate is equal to that of population \(2\) adjusted for the covariate and so on, and thus all population means adjusted for the covariate are equal to an overall adjusted mean. If the effect of the \(i\)-th group is the difference between the \(i\)-th group adjusted mean and the overall adjusted mean (\(\alpha_i(adj)=\mu_i(adj)−\mu(adj)\)) then the \(H_0\) can alternatively be written as:
The effect of each group equals zero. If one or more of the \(\alpha_i(adj)\) are different from zero (the response mean for this treatment differs from the overall response mean), the null hypothesis is not true, indicating that the treatment does affect the response variable.
Factor B: the covariate effect
The pooled population slope equals zero. Note, that this null hypothesis is rarely of much interest. It is precisely because of this nuisance relationship that ANCOVA designs are applied.
One or more covariates can be incorporated into single factor, nested, factorial and partly nested designs in order to reduce the unexplained variation. Fundamentally, the covariate(s) are purely used to adjust the response values prior to the regular analysis. The difficulty is in determining the appropriate adjustments. Following is a list of the appropriate linear models and adjusted response calculations for a range of ANCOVA designs. Note that these linear models do not include interactions involving the covariates as these are assumed to be zero. The inclusion of these interaction terms is a useful means of testing the homogeneity of slopes assumption.
Single categorical and single covariate
Linear model: \(y_{ij}=\mu + \alpha_i + \beta(x_{ij}-\bar{x}) + \epsilon_{ij}\)
Adjustments: \(y_{ij(adj)}=y_{ij} - b(x_{ij} - \bar{x})\)
Single categorical and two covariates
Linear model: \(y_{ij}=\mu + \alpha_i + \beta_{YX}(x_{ij}-\bar{x}) + \beta_{YZ}(z_{ij}-\bar{z}) + \epsilon_{ij}\)
Adjustments: \(y_{ij(adj)}=y_{ij} - b_{YX}(x_{ij} - \bar{x}) - b_{YZ}(z_{ij} - \bar{z})\)
Factorial designs
Linear model: \(y_{ij}=\mu + \alpha_i + \gamma_j + (\alpha\gamma)_{ij}+ \beta(x_{ijk}-\bar{x}) + \epsilon_{ijk}\)
Adjustments: \(y_{ijk(adj)}=y_{ijk} - b(x_{ijk} - \bar{x})\)
Nested designs
Linear model: \(y_{ijk}=\mu + \alpha_i + \gamma_{j(i)} + \beta(x_{ijk}-\bar{x}) + \epsilon_{ijk}\)
Adjustments: \(y_{ijk(adj)}=y_{ijk} - b(x_{ijk} - \bar{x})\)
Partly nested designs
Linear model: \(y_{ijkl}=\mu + \alpha_i + \gamma_{j(i)} + \delta_k + (\alpha\delta)_{ik} + (\gamma\delta)_{j(i)k} + \beta(x_{ijk}-\bar{x}) + \epsilon_{ijkl}\)
Adjustments: \(y_{ijk(adj)}=y_{ijkl} - b_{between}(x_{i} - \bar{x}) - b_{within}(x_{ijk} - \bar{x}_i)\)
In ANCOVA, the total variability of the response variable is sequentially partitioned into components explained by each of the model terms, starting with the covariate and is therefore equivalent to performing a regular analysis of variance on the response variables that have been adjusted for the covariate. The appropriate unexplained residuals and therefore the appropriate F-ratios for each factor differ according to the different null hypotheses associated with different linear models as well as combinations of fixed and random factors in the model (see the following tables). Note that since the covariate levels measured are typically different for each group, ANCOVA designs are inherently non-orthogonal (unbalanced). Consequently, sequential (Type I sums of squares) should not be used. For very simple Ancova designs that incorporate a single categorical and single covariate, Type I sums of squares can be used provided the covariate appears in the linear model first (and thus is partitioned out last) as we are typically not interested in estimating this effect.
ancova_tableNA df MS F-ratio (A&B fixed) F-ratio (B fixed)
NA Factor A "a-1" "MS A" "(MS A)/(MS res)" "(MS A)/(MS res)"
NA Factor B "1" "MS B" "(MS B)/(MS res)" "(MS B)/(MS res)"
NA Factor AB "a-1" "MS AB" "(MS AB)/(MS res)" "(MS AB)/(MS res)"
NA Residual "(n-2)a" "MS res" "" ""The corresponding R syntax is given below.
anova(lm(DV ~ B * A, dataset))
# OR
anova(aov(DV ~ B * A, dataset))
# OR (make sure not using treatment contrasts)
Anova(lm(DV ~ B * A, dataset), type = "III")As ANCOVA designs are essentially regular ANOVA designs that are first adjusted (centered) for the covariate(s), ANCOVA designs inherit all of the underlying assumptions of the appropriate ANOVA design. Specifically, hypothesis tests assume that:
The appropriate residuals are normally distributed. Boxplots using the appropriate scale of replication (reflecting the appropriate residuals/F-ratio denominator, see the above tables) should be used to explore normality. Scale transformations are often useful.
The appropriate residuals are equally varied. Boxplots and plots of means against variance (using the appropriate scale of replication) should be used to explore the spread of values. Residual plots should reveal no patterns. Scale transformations are often useful.
The appropriate residuals are independent of one another.
The relationship between the response variable and the covariate should be linear. Linearity can be explored using scatterplots and residual plots should reveal no patterns.
For repeated measures and other designs in which treatment levels within blocks can not be be randomly ordered, the variance/covariance matrix is assumed to display sphericity.
For designs that utilise blocking, it is assumed that there are no block by within block interactions.
Homogeneity of Slopes
In addition to the above assumptions, ANCOVA designs also assume that slopes of relationships between the response variable and the covariate(s) are the same for each treatment level (group). That is, all the trends are parallel. If the individual slopes deviate substantially from each other (and thus the overall slope), then adjustments made to each of the observations are nonsensical. This situation is analogous to an interaction between two or more factors. In ANCOVA, interactions involving the covariate suggest that the nature of the relationship between the response and the covariate differs between the levels of the categorical treatment. More importantly, they also indicate that whether or not there is an effect of the treatment depends on what range of the covariate you are focussed on. Clearly then, it is not possible to make conclusions about the main effects of treatments in the presence of such interactions. The assumption of homogeneity of slopes can be examined via interaction plots or more formally, by testing hypotheses about the interactions between categorical variables and the covariate(s). There are three broad approaches for dealing with ANCOVA designs with heterogeneous slopes and selection depends on the primary focus of the study.
When the primary objective of the analysis is to investigate the effects of categorical treatments, it is possible to adopt an approach similar to that taken when exploring interactions in multiple regression. The effect of treatments can be examined at specific values of the covariate (such as the mean and \(\pm\) one standard deviation). This approach is really only useful at revealing broad shifts in patterns over the range of the covariate and if the selected values of the covariate do not have some inherent clinical meaning (selected arbitrarily), then the outcomes can be of only limited clinical interest.
Alternatively, the Johnson-Neyman technique (or Wilxon modification thereof) procedure indicates the ranges of the covariate over which the individual regression lines of pairs of treatment groups overlap or cross. Although less powerful than the previous approach, the Wilcox(J-N) procedure has the advantage of revealing the important range (ranges for which the groups are different and not different) of the covariate rather than being constrained by specific levels selected.
Use contrast treatments to split up the interaction term into its constituent contrasts for each level of the treatment. Essentially this compares each of the treatment level slopes to the slope from the “control” group and is useful if the primary focus is on the relationships between the response and the covariate.
Similar covariate ranges
Adjustments made to the response means in an attempt to statistically account for differences in the covariate involve predicting mean response values along displaced linear relationships between the overall response and covariate variables. The degree of trend displacement for any given group is essentially calculated by multiplying the overall regression slope by the degree of difference between the overall covariate mean and the mean of the covariate for that group. However, when the ranges of the covariate within each of the groups differ substantially from one another, these adjustments are effectively extrapolations and therefore of unknown reliability. If a simple ANOVA of the covariate modelled against the categorical factor indicates that the covariate means differ significantly between groups, it may be necessary to either remove extreme observations or reconsider the analysis.
Robust ANCOVA
ANCOVA based on rank transformed data can be useful for accommodating data with numerous problematic outliers. Nevertheless, problems about the difficulties of detecting interactions from rank transformed data, obviously have implications for inferential tests of homogeneity of slopes. Randomisation tests that maintain response0covariate pairs and repeatedly randomise these observations amongst the levels of the treatments can also be useful, particularly when there is doubt over the independence of observations. Both planned and unplanned comparisons follow those of other ANOVA chapters without any real additional complications. Notably, recent implementations of the Tukey’s test (within R) accommodate unbalanced designs and thus negate the need for some of the more complicated and specialised techniques that have been highlighted in past texts.
Consider an experimental design aimed at exploring the effects of a categorical variable with three levels (Group A, Group B and Group C) on a response. From previous studies, we know that the response is influenced by another variable (covariate). Unfortunately, it was not possible to ensure that all sampling units were the same degree of the covariate. Therefore, in an attempt to account for this anticipated extra source of variability, we measured the level of the covariate for each sampling unit. Actually, in allocating treatments to the various treatment groups, we tried to ensure a similar mean and range of the covariate within each group.
set.seed(123)
n <- 10
p <- 3
A.eff <- c(40, -15, -20)
beta <- -0.45
sigma <- 4
B <- rnorm(n * p, 0, 15)
A <- gl(p, n, lab = paste("Group", LETTERS[1:3]))
mm <- model.matrix(~A + B)
data <- data.frame(A = A, B = B, Y = as.numeric(c(A.eff, beta) %*% t(mm)) + rnorm(n * p, 0, 4))
data$B <- data$B + 20
head(data)NA A B Y
NA 1 Group A 11.59287 45.48907
NA 2 Group A 16.54734 40.37341
NA 3 Group A 43.38062 33.05922
NA 4 Group A 21.05763 43.03660
NA 5 Group A 21.93932 42.41363
NA 6 Group A 45.72597 31.17787library(car)
scatterplot(Y ~ B | A, data = data)
boxplot(Y ~ A, data)
# OR via ggplot
library(ggplot2)
ggplot(data, aes(y = Y, x = B, group = A)) + geom_point() + geom_smooth(method = "lm")
ggplot(data, aes(y = Y, x = A)) + geom_boxplot()
Conclusions
There is no evidence of obvious non-normality. The assumption of linearity seems reasonable. The variability of the three groups seems approximately equal. The slopes (\(Y\) vs B trends) appear broadly similar for each treatment group.
We can explore inferential evidence of unequal slopes by examining estimated effects of the interaction between the categorical variable and the covariate. Note, pay no attention to the main effects - only the interaction. Even though I intend to illustrate Bayesian analyses here, for such a simple model, it is considerably simpler to use traditional OLS for testing for the presence of an interaction.
anova(lm(Y ~ B * A, data = data))NA Analysis of Variance Table
NA
NA Response: Y
NA Df Sum Sq Mean Sq F value Pr(>F)
NA B 1 989.99 989.99 92.6782 1.027e-09 ***
NA A 2 2320.05 1160.02 108.5956 9.423e-13 ***
NA B:A 2 51.36 25.68 2.4041 0.1118
NA Residuals 24 256.37 10.68
NA ---
NA Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There is very little evidence to suggest that the assumption of equal slopes will be inappropriate.
The observed response (\(y_i\)) are assumed to be drawn from a normal distribution with a given mean (\(\mu\)) and standard deviation (\(\sigma\)). The expected values are themselves determined by the linear predictor (\(\boldsymbol X \boldsymbol \beta\)). In this case, \(\boldsymbol \beta\) represents the vector of \(\beta\)’s - the intercept associated with the first group, the (effects) differences between this intercept and the intercepts for each other group as well as the slope associated with the continuous covariate. \(\boldsymbol X\) is the model matrix. MCMC sampling requires priors on all parameters. We will employ weakly informative priors. Specifying ‘uninformative’ priors is always a bit of a balancing act. If the priors are too vague (wide) the MCMC sampler can wander off into nonscence areas of likelihood rather than concentrate around areas of highest likelihood (desired when wanting the outcomes to be largely driven by the data). On the other hand, if the priors are too strong, they may have an influence on the parameters. In such a simple model, this balance is very forgiving - it is for more complex models that prior choice becomes more important. For this simple model, we will go with zero-centered Gaussian (normal) priors with relatively large standard deviations (\(100\)) for both the intercept and the treatment effect and a wide half-cauchy (\(\text{scale}=5\)) for the standard deviation.
\[ y_i \sim N(\mu_i,\sigma), \]
where \(\mu_i=\beta_0 +\boldsymbol \beta \boldsymbol X\). The assumed priors are: \(\beta \sim N(0,100)\) and \(\sigma \sim \text{Cauchy}(0,5)\). Note, exploratory data analysis suggests that while the intercept (intercept of Group A) and categorical predictor effects (differences between intercepts of each of the Group and Group A’s intercept) could be drawn from a similar distribution (with mean in the \(10\)’s and variances in the \(100\)’s), the slope (effect associated with Group A linear relationship) is likely to be an order of magnitude less. We might therefore be tempted to provide different priors for the intercept, categorical effects and slope effect. For a simple model such as this, it is unlikely to be necessary. However, for more complex models, where prior specification becomes more critical, separate priors would probably be necessary. We proceed to code the model into Stan.
modelString = "
data {
int<lower=1> n;
int<lower=1> nX;
vector [n] y;
matrix [n,nX] X;
}
parameters {
vector[nX] beta;
real<lower=0> sigma;
}
transformed parameters {
vector[n] mu;
mu = X*beta;
}
model {
//Likelihood
y~normal(mu,sigma);
//Priors
beta ~ normal(0,100);
sigma~cauchy(0,5);
}
generated quantities {
vector[n] log_lik;
for (i in 1:n) {
log_lik[i] = normal_lpdf(y[i] | mu[i], sigma);
}
}
"
## write the model to a stan file
writeLines(modelString, con = "ancovaModel.stan")Arrange the data as a list (as required by Stan). As input, Stan will need to be supplied with: the response variable, the predictor variable, the total number of observed items. This all needs to be contained within a list object. We will create two data lists, one for each of the hypotheses.
Xmat <- model.matrix(~A + B, data)
data.list <- with(data, list(y = Y, X = Xmat, nX = ncol(Xmat), n = nrow(data)))Define the nodes (parameters and derivatives) to monitor and chain parameters.
params <- c("beta", "sigma", "log_lik")
nChains = 2
burnInSteps = 500
thinSteps = 1
numSavedSteps = 2000 #across all chains
nIter = ceiling(burnInSteps + (numSavedSteps * thinSteps)/nChains)
nIterNA [1] 1500Now compile and run the Stan code via the rstan interface. Note that the first time stan is run after the rstan package is loaded, it is often necessary to run any kind of randomization function just to initiate the .Random.seed variable.
library(rstan)During the warmup stage, the No-U-Turn sampler (NUTS) attempts to determine the optimum stepsize - the stepsize that achieves the target acceptance rate (\(0.8\) or \(80\)% by default) without divergence (occurs when the stepsize is too large relative to the curvature of the log posterior and results in approximations that are likely to diverge and be biased) - and without hitting the maximum treedepth (\(10\)). At each iteration of the NUTS algorithm, the number of leapfrog steps doubles (as it increases the treedepth) and only terminates when either the NUTS criterion are satisfied or the tree depth reaches the maximum (\(10\) by default).
data.rstan <- stan(data = data.list, file = "ancovaModel.stan", chains = nChains, pars = params, iter = nIter, warmup = burnInSteps, thin = thinSteps)NA
NA SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
NA Chain 1:
NA Chain 1: Gradient evaluation took 2.6e-05 seconds
NA Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.26 seconds.
NA Chain 1: Adjust your expectations accordingly!
NA Chain 1:
NA Chain 1:
NA Chain 1: Iteration: 1 / 1500 [ 0%] (Warmup)
NA Chain 1: Iteration: 150 / 1500 [ 10%] (Warmup)
NA Chain 1: Iteration: 300 / 1500 [ 20%] (Warmup)
NA Chain 1: Iteration: 450 / 1500 [ 30%] (Warmup)
NA Chain 1: Iteration: 501 / 1500 [ 33%] (Sampling)
NA Chain 1: Iteration: 650 / 1500 [ 43%] (Sampling)
NA Chain 1: Iteration: 800 / 1500 [ 53%] (Sampling)
NA Chain 1: Iteration: 950 / 1500 [ 63%] (Sampling)
NA Chain 1: Iteration: 1100 / 1500 [ 73%] (Sampling)
NA Chain 1: Iteration: 1250 / 1500 [ 83%] (Sampling)
NA Chain 1: Iteration: 1400 / 1500 [ 93%] (Sampling)
NA Chain 1: Iteration: 1500 / 1500 [100%] (Sampling)
NA Chain 1:
NA Chain 1: Elapsed Time: 0.045 seconds (Warm-up)
NA Chain 1: 0.035 seconds (Sampling)
NA Chain 1: 0.08 seconds (Total)
NA Chain 1:
NA
NA SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
NA Chain 2:
NA Chain 2: Gradient evaluation took 4e-06 seconds
NA Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
NA Chain 2: Adjust your expectations accordingly!
NA Chain 2:
NA Chain 2:
NA Chain 2: Iteration: 1 / 1500 [ 0%] (Warmup)
NA Chain 2: Iteration: 150 / 1500 [ 10%] (Warmup)
NA Chain 2: Iteration: 300 / 1500 [ 20%] (Warmup)
NA Chain 2: Iteration: 450 / 1500 [ 30%] (Warmup)
NA Chain 2: Iteration: 501 / 1500 [ 33%] (Sampling)
NA Chain 2: Iteration: 650 / 1500 [ 43%] (Sampling)
NA Chain 2: Iteration: 800 / 1500 [ 53%] (Sampling)
NA Chain 2: Iteration: 950 / 1500 [ 63%] (Sampling)
NA Chain 2: Iteration: 1100 / 1500 [ 73%] (Sampling)
NA Chain 2: Iteration: 1250 / 1500 [ 83%] (Sampling)
NA Chain 2: Iteration: 1400 / 1500 [ 93%] (Sampling)
NA Chain 2: Iteration: 1500 / 1500 [100%] (Sampling)
NA Chain 2:
NA Chain 2: Elapsed Time: 0.039 seconds (Warm-up)
NA Chain 2: 0.034 seconds (Sampling)
NA Chain 2: 0.073 seconds (Total)
NA Chain 2:print(data.rstan, par = c("beta", "sigma"))NA Inference for Stan model: anon_model.
NA 2 chains, each with iter=1500; warmup=500; thin=1;
NA post-warmup draws per chain=1000, total post-warmup draws=2000.
NA
NA mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
NA beta[1] 50.93 0.05 1.54 47.90 49.90 50.93 51.95 53.92 964 1
NA beta[2] -16.18 0.05 1.61 -19.24 -17.25 -16.20 -15.10 -13.00 1241 1
NA beta[3] -20.56 0.05 1.63 -23.79 -21.66 -20.51 -19.47 -17.36 1134 1
NA beta[4] -0.48 0.00 0.05 -0.57 -0.51 -0.48 -0.45 -0.40 1194 1
NA sigma 3.54 0.01 0.51 2.71 3.19 3.49 3.82 4.75 1265 1
NA
NA Samples were drawn using NUTS(diag_e) at Mon Jul 22 11:49:42 2024.
NA For each parameter, n_eff is a crude measure of effective sample size,
NA and Rhat is the potential scale reduction factor on split chains (at
NA convergence, Rhat=1).In addition to the regular model diagnostic checks (such as residual plots), for Bayesian analyses, it is necessary to explore the characteristics of the MCMC chains and the sampler in general. Recall that the purpose of MCMC sampling is to replicate the posterior distribution of the model likelihood and priors by drawing a known number of samples from this posterior (thereby formulating a probability distribution). This is only reliable if the MCMC samples accurately reflect the posterior. Unfortunately, since we only know the posterior in the most trivial of circumstances, it is necessary to rely on indirect measures of how accurately the MCMC samples are likely to reflect the likelihood. I will briefly outline the most important diagnostics.
Traceplots for each parameter illustrate the MCMC sample values after each successive iteration along the chain. Bad chain mixing (characterised by any sort of pattern) suggests that the MCMC sampling chains may not have completely traversed all features of the posterior distribution and that more iterations are required to ensure the distribution has been accurately represented.
Autocorrelation plot for each parameter illustrate the degree of correlation between MCMC samples separated by different lags. For example, a lag of \(0\) represents the degree of correlation between each MCMC sample and itself (obviously this will be a correlation of \(1\)). A lag of \(1\) represents the degree of correlation between each MCMC sample and the next sample along the chain and so on. In order to be able to generate unbiased estimates of parameters, the MCMC samples should be independent (uncorrelated).
Potential scale reduction factor (Rhat) statistic for each parameter provides a measure of sampling efficiency/effectiveness. Ideally, all values should be less than \(1.05\). If there are values of \(1.05\) or greater it suggests that the sampler was not very efficient or effective. Not only does this mean that the sampler was potentially slower than it could have been but, more importantly, it could indicate that the sampler spent time sampling in a region of the likelihood that is less informative. Such a situation can arise from either a misspecified model or overly vague priors that permit sampling in otherwise nonscence parameter space.
Prior to examining the summaries, we should have explored the convergence diagnostics. We use the package mcmcplots to obtain density and trace plots.
library(mcmcplots)
s = as.array(data.rstan)
wch = grep("beta", dimnames(s)$parameters)
s = s[, , wch]
mcmc <- do.call(mcmc.list, plyr:::alply(s[, , -(length(s[1, 1, ]))], 2, as.mcmc))
denplot(mcmc, parms = c("beta"))
traplot(mcmc, parms = c("beta"))
These plots show no evidence that the chains have not reasonably traversed the entire multidimensional parameter space.
#Raftery diagnostic
raftery.diag(mcmc)NA $`1`
NA
NA Quantile (q) = 0.025
NA Accuracy (r) = +/- 0.005
NA Probability (s) = 0.95
NA
NA You need a sample size of at least 3746 with these values of q, r and s
NA
NA $`2`
NA
NA Quantile (q) = 0.025
NA Accuracy (r) = +/- 0.005
NA Probability (s) = 0.95
NA
NA You need a sample size of at least 3746 with these values of q, r and sThe Raftery diagnostics for each chain estimate that we would require no more than \(5000\) samples to reach the specified level of confidence in convergence. As we have \(10500\) samples, we can be confidence that convergence has occurred.
#Autocorrelation diagnostic
stan_ac(data.rstan, pars = c("beta"))
A lag of 10 appears to be sufficient to avoid autocorrelation (poor mixing).
stan_rhat(data.rstan, pars = c("beta"))
stan_ess(data.rstan, pars = c("beta"))
Rhat and effective sample size. In this instance, most of the parameters have reasonably high effective samples and thus there is likely to be a good range of values from which to estimate parameter properties.
Model validation involves exploring the model diagnostics and fit to ensure that the model is broadly appropriate for the data. As such, exploration of the residuals should be routine. Ideally, a good model should also be able to predict the data used to fit the model. Residuals are not computed directly within rstan However, we can calculate them manually form the posteriors.
library(dplyr)
mcmc = as.data.frame(data.rstan) %>% dplyr:::select(contains("beta"),
sigma) %>% as.matrix
# generate a model matrix
newdata = data
Xmat = model.matrix(~A + B, newdata)
## get median parameter estimates
coefs = apply(mcmc[, 1:4], 2, median)
fit = as.vector(coefs %*% t(Xmat))
resid = data$Y - fit
ggplot() + geom_point(data = NULL, aes(y = resid, x = fit)) + theme_classic()
Residuals against predictors
library(tidyr)
mcmc = as.data.frame(data.rstan) %>% dplyr:::select(contains("beta"),
sigma) %>% as.matrix
# generate a model matrix
newdata = newdata
Xmat = model.matrix(~A + B, newdata)
## get median parameter estimates
coefs = apply(mcmc[, 1:4], 2, median)
fit = as.vector(coefs %*% t(Xmat))
resid = data$Y - fit
newdata = newdata %>% cbind(fit, resid)
ggplot(newdata) + geom_point(aes(y = resid, x = A)) + theme_classic()
ggplot(newdata) + geom_point(aes(y = resid, x = B)) + theme_classic()
And now for studentised residuals
mcmc = as.data.frame(data.rstan) %>% dplyr:::select(contains("beta"),
sigma) %>% as.matrix
# generate a model matrix
newdata = data
Xmat = model.matrix(~A + B, newdata)
## get median parameter estimates
coefs = apply(mcmc[, 1:4], 2, median)
fit = as.vector(coefs %*% t(Xmat))
resid = data$Y - fit
sresid = resid/sd(resid)
ggplot() + geom_point(data = NULL, aes(y = sresid, x = fit)) + theme_classic()
For this simple model, the studentized residuals yield the same pattern as the raw residuals (or the Pearson residuals for that matter). Lets see how well data simulated from the model reflects the raw data.
mcmc = as.data.frame(data.rstan) %>% dplyr:::select(contains("beta"),
sigma) %>% as.matrix
# generate a model matrix
Xmat = model.matrix(~A + B, data)
## get median parameter estimates
coefs = mcmc[, 1:4]
fit = coefs %*% t(Xmat)
## draw samples from this model
yRep = sapply(1:nrow(mcmc), function(i) rnorm(nrow(data), fit[i,
], mcmc[i, "sigma"]))
newdata = data.frame(A = data$A, B = data$B, yRep) %>% gather(key = Sample,
value = Value, -A, -B)
ggplot(newdata) + geom_violin(aes(y = Value, x = A, fill = "Model"),
alpha = 0.5) + geom_violin(data = data, aes(y = Y, x = A,
fill = "Obs"), alpha = 0.5) + geom_point(data = data, aes(y = Y,
x = A), position = position_jitter(width = 0.1, height = 0),
color = "black") + theme_classic()
ggplot(newdata) + geom_violin(aes(y = Value, x = B, fill = "Model",
group = B, color = A), alpha = 0.5) + geom_point(data = data,
aes(y = Y, x = B, group = B, color = A)) + theme_classic()
The predicted trends do encapsulate the actual data, suggesting that the model is a reasonable representation of the underlying processes. Note, these are prediction intervals rather than confidence intervals as we are seeking intervals within which we can predict individual observations rather than means. We can also explore the posteriors of each parameter.
library(bayesplot)
mcmc_intervals(as.matrix(data.rstan), regex_pars = "beta|sigma")
mcmc_areas(as.matrix(data.rstan), regex_pars = "beta|sigma")
Although all parameters in a Bayesian analysis are considered random and are considered a distribution, rarely would it be useful to present tables of all the samples from each distribution. On the other hand, plots of the posterior distributions have some use. Nevertheless, most workers prefer to present simple statistical summaries of the posteriors. Popular choices include the median (or mean) and \(95\)% credibility intervals.
mcmcpvalue <- function(samp) {
## elementary version that creates an empirical p-value for the
## hypothesis that the columns of samp have mean zero versus a general
## multivariate distribution with elliptical contours.
## differences from the mean standardized by the observed
## variance-covariance factor
## Note, I put in the bit for single terms
if (length(dim(samp)) == 0) {
std <- backsolve(chol(var(samp)), cbind(0, t(samp)) - mean(samp),
transpose = TRUE)
sqdist <- colSums(std * std)
sum(sqdist[-1] > sqdist[1])/length(samp)
} else {
std <- backsolve(chol(var(samp)), cbind(0, t(samp)) - colMeans(samp),
transpose = TRUE)
sqdist <- colSums(std * std)
sum(sqdist[-1] > sqdist[1])/nrow(samp)
}
}First, we look at the results from the additive model.
print(data.rstan, pars = c("beta", "sigma"))NA Inference for Stan model: anon_model.
NA 2 chains, each with iter=1500; warmup=500; thin=1;
NA post-warmup draws per chain=1000, total post-warmup draws=2000.
NA
NA mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
NA beta[1] 50.93 0.05 1.54 47.90 49.90 50.93 51.95 53.92 964 1
NA beta[2] -16.18 0.05 1.61 -19.24 -17.25 -16.20 -15.10 -13.00 1241 1
NA beta[3] -20.56 0.05 1.63 -23.79 -21.66 -20.51 -19.47 -17.36 1134 1
NA beta[4] -0.48 0.00 0.05 -0.57 -0.51 -0.48 -0.45 -0.40 1194 1
NA sigma 3.54 0.01 0.51 2.71 3.19 3.49 3.82 4.75 1265 1
NA
NA Samples were drawn using NUTS(diag_e) at Mon Jul 22 11:49:42 2024.
NA For each parameter, n_eff is a crude measure of effective sample size,
NA and Rhat is the potential scale reduction factor on split chains (at
NA convergence, Rhat=1).# OR
library(broom)
library(broom.mixed)
tidyMCMC(data.rstan, conf.int = TRUE, conf.method = "HPDinterval", pars = c("beta", "sigma"))NA # A tibble: 5 × 5
NA term estimate std.error conf.low conf.high
NA <chr> <dbl> <dbl> <dbl> <dbl>
NA 1 beta[1] 50.9 1.54 47.8 53.9
NA 2 beta[2] -16.2 1.61 -19.3 -13.0
NA 3 beta[3] -20.6 1.63 -24.0 -17.5
NA 4 beta[4] -0.483 0.0465 -0.575 -0.397
NA 5 sigma 3.54 0.506 2.65 4.56Conclusions
The intercept of the first group (Group A) is \(51\).
The mean of the second group (Group B) is \(-16.3\) units greater than (A).
The mean of the third group (Group C) is \(-20.7\) units greater than (A).
A one unit increase in B in Group A is associated with a \(-0.484\) units increase in \(Y\).
The \(95\)% confidence interval for the effects of Group B, Group C and the partial slope associated with B do not overlapp with 0 implying a significant difference between group A and groups B, C and a significant negative relationship with B. While workers attempt to become comfortable with a new statistical framework, it is only natural that they like to evaluate and comprehend new structures and output alongside more familiar concepts. One way to facilitate this is via Bayesian p-values that are somewhat analogous to the frequentist p-values for investigating the hypothesis that a parameter is equal to zero.
## since values are less than zero
mcmcpvalue(as.matrix(data.rstan)[, "beta[2]"]) # effect of (B-A = 0)NA [1] 0mcmcpvalue(as.matrix(data.rstan)[, "beta[3]"]) # effect of (C-A = 0)NA [1] 0mcmcpvalue(as.matrix(data.rstan)[, "beta[4]"]) # effect of (slope = 0)NA [1] 0mcmcpvalue(as.matrix(data.rstan)[, 2:4]) # effect of (model)NA [1] 0There is evidence that the reponse differs between the groups. There is evidence suggesting that the response of group D differs from that of group A. In a Bayesian context, we can compare models using the leave-one-out cross-validation statistics. Leave-one-out (LOO) cross-validation explores how well a series of models can predict withheld values Vehtari, Gelman, and Gabry (2017). The LOO Information Criterion (LOOIC) is analogous to the AIC except that the LOOIC takes priors into consideration, does not assume that the posterior distribution is drawn from a multivariate normal and integrates over parameter uncertainty so as to yield a distribution of looic rather than just a point estimate. The LOOIC does however assume that all observations are equally influential (it does not matter which observations are left out). This assumption can be examined via the Pareto \(k\) estimate (values greater than \(0.5\) or more conservatively \(0.75\) are considered overly influential). We can compute LOOIC if we store the loglikelihood from our Stan model, which can then be extracted to compute the information criterion using the package loo.
library(loo)
(full = loo(extract_log_lik(data.rstan)))NA
NA Computed from 2000 by 30 log-likelihood matrix.
NA
NA Estimate SE
NA elpd_loo -83.0 4.3
NA p_loo 4.8 1.3
NA looic 166.0 8.7
NA ------
NA MCSE of elpd_loo is NA.
NA MCSE and ESS estimates assume independent draws (r_eff=1).
NA
NA Pareto k diagnostic values:
NA Count Pct. Min. ESS
NA (-Inf, 0.7] (good) 29 96.7% 300
NA (0.7, 1] (bad) 1 3.3% <NA>
NA (1, Inf) (very bad) 0 0.0% <NA>
NA See help('pareto-k-diagnostic') for details.# now fit a model without main factor
modelString2 = "
data {
int<lower=1> n;
int<lower=1> nX;
vector [n] y;
matrix [n,nX] X;
}
parameters {
vector[nX] beta;
real<lower=0> sigma;
}
transformed parameters {
vector[n] mu;
mu = X*beta;
}
model {
// Likelihood
y~normal(mu,sigma);
// Priors
beta ~ normal(0,1000);
sigma~cauchy(0,5);
}
generated quantities {
vector[n] log_lik;
for (i in 1:n) {
log_lik[i] = normal_lpdf(y[i] | mu[i], sigma);
}
}
"
## write the model to a stan file
writeLines(modelString2, con = "ancovaModel2.stan")
Xmat <- model.matrix(~1, data)
data.list <- with(data, list(y = Y, X = Xmat, n = nrow(data), nX = ncol(Xmat)))
data.rstan.red <- stan(data = data.list, file = "ancovaModel2.stan", chains = nChains,
iter = nIter, warmup = burnInSteps, thin = thinSteps)NA
NA SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
NA Chain 1:
NA Chain 1: Gradient evaluation took 2.6e-05 seconds
NA Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.26 seconds.
NA Chain 1: Adjust your expectations accordingly!
NA Chain 1:
NA Chain 1:
NA Chain 1: Iteration: 1 / 1500 [ 0%] (Warmup)
NA Chain 1: Iteration: 150 / 1500 [ 10%] (Warmup)
NA Chain 1: Iteration: 300 / 1500 [ 20%] (Warmup)
NA Chain 1: Iteration: 450 / 1500 [ 30%] (Warmup)
NA Chain 1: Iteration: 501 / 1500 [ 33%] (Sampling)
NA Chain 1: Iteration: 650 / 1500 [ 43%] (Sampling)
NA Chain 1: Iteration: 800 / 1500 [ 53%] (Sampling)
NA Chain 1: Iteration: 950 / 1500 [ 63%] (Sampling)
NA Chain 1: Iteration: 1100 / 1500 [ 73%] (Sampling)
NA Chain 1: Iteration: 1250 / 1500 [ 83%] (Sampling)
NA Chain 1: Iteration: 1400 / 1500 [ 93%] (Sampling)
NA Chain 1: Iteration: 1500 / 1500 [100%] (Sampling)
NA Chain 1:
NA Chain 1: Elapsed Time: 0.009 seconds (Warm-up)
NA Chain 1: 0.014 seconds (Sampling)
NA Chain 1: 0.023 seconds (Total)
NA Chain 1:
NA
NA SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
NA Chain 2:
NA Chain 2: Gradient evaluation took 4e-06 seconds
NA Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
NA Chain 2: Adjust your expectations accordingly!
NA Chain 2:
NA Chain 2:
NA Chain 2: Iteration: 1 / 1500 [ 0%] (Warmup)
NA Chain 2: Iteration: 150 / 1500 [ 10%] (Warmup)
NA Chain 2: Iteration: 300 / 1500 [ 20%] (Warmup)
NA Chain 2: Iteration: 450 / 1500 [ 30%] (Warmup)
NA Chain 2: Iteration: 501 / 1500 [ 33%] (Sampling)
NA Chain 2: Iteration: 650 / 1500 [ 43%] (Sampling)
NA Chain 2: Iteration: 800 / 1500 [ 53%] (Sampling)
NA Chain 2: Iteration: 950 / 1500 [ 63%] (Sampling)
NA Chain 2: Iteration: 1100 / 1500 [ 73%] (Sampling)
NA Chain 2: Iteration: 1250 / 1500 [ 83%] (Sampling)
NA Chain 2: Iteration: 1400 / 1500 [ 93%] (Sampling)
NA Chain 2: Iteration: 1500 / 1500 [100%] (Sampling)
NA Chain 2:
NA Chain 2: Elapsed Time: 0.01 seconds (Warm-up)
NA Chain 2: 0.015 seconds (Sampling)
NA Chain 2: 0.025 seconds (Total)
NA Chain 2:(reduced = loo(extract_log_lik(data.rstan.red)))NA
NA Computed from 2000 by 30 log-likelihood matrix.
NA
NA Estimate SE
NA elpd_loo -116.4 3.1
NA p_loo 1.7 0.4
NA looic 232.9 6.2
NA ------
NA MCSE of elpd_loo is 0.0.
NA MCSE and ESS estimates assume independent draws (r_eff=1).
NA
NA All Pareto k estimates are good (k < 0.7).
NA See help('pareto-k-diagnostic') for details.par(mfrow = 1:2, mar = c(5, 3.8, 1, 0) + 0.1, las = 3)
plot(full, label_points = TRUE)
plot(reduced, label_points = TRUE)
The expected out-of-sample predictive accuracy is substantially lower for the model that includes \(x\). This might be used to suggest that the inferential evidence for a general effect of \(x\) on \(y\).
A nice graphic is often a great accompaniment to a statistical analysis. Although there are no fixed assumptions associated with graphing (in contrast to statistical analyses), we often want the graphical summaries to reflect the associated statistical analyses. After all, the sample is just one perspective on the population(s). What we are more interested in is being able to estimate and depict likely population parameters/trends. Thus, whilst we could easily provide a plot displaying the raw data along with simple measures of location and spread, arguably, we should use estimates that reflect the fitted model. In this case, it would be appropriate to plot the credibility interval associated with each group.
mcmc = as.matrix(data.rstan)
## Calculate the fitted values
newdata = expand.grid(A = levels(data$A), B = seq(min(data$B), max(data$B),
len = 100))
Xmat = model.matrix(~A + B, newdata)
coefs = mcmc[, c("beta[1]", "beta[2]", "beta[3]", "beta[4]")]
fit = coefs %*% t(Xmat)
newdata = newdata %>% cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval"))
ggplot(newdata, aes(y = estimate, x = B, fill = A)) + geom_ribbon(aes(ymin = conf.low,
ymax = conf.high), alpha = 0.2) + geom_line() + scale_y_continuous("Y") +
scale_x_continuous("B") + theme_classic()
As this is simple single factor ANOVA, we can simple add the raw data to this figure. For more complex designs with additional predictors, it is necessary to plot partial residuals.
## Calculate partial residuals fitted values
fdata = rdata = data
fMat = rMat = model.matrix(~A + B, fdata)
fit = as.vector(apply(coefs, 2, median) %*% t(fMat))
resid = as.vector(data$Y - apply(coefs, 2, median) %*% t(rMat))
rdata = rdata %>% mutate(partial.resid = resid + fit)
ggplot(newdata, aes(y = estimate, x = B, fill = A)) + geom_point(data = rdata,
aes(y = partial.resid, x = B, color = A)) + geom_ribbon(aes(ymin = conf.low,
ymax = conf.high), alpha = 0.2) + geom_line() + scale_y_continuous("Y") +
scale_x_continuous("B") + theme_classic()
In frequentist statistics, when we have more than two groups, we are typically not only interested in whether there is evidence for an overall “effect” of a factor - we are also interested in how various groups compare to one another. To explore these trends, we either compare each group to each other in a pairwise manner (controlling for family-wise Type I error rates) or we explore an independent subset of the possible comparisons. Although these alternate approaches can adequately address a specific research agenda, often they impose severe limitations and compromises on the scope and breadth of questions that can be asked of your data. The reason for these limitations is that in a frequentist framework, any single hypothesis carries with it a (nominally) \(5\)% chance of a false rejection (since it is based on long-run frequency). Thus, performing multiple tests are likely to compound this error rate. The point is, that each comparison is compared to its own probability distribution (and each carries a \(5\)% error rate). By contrast, in Bayesian statistics, all comparisons (contrasts) are drawn from the one (hopefully stable and convergent) posterior distribution and this posterior is invariant to the type and number of comparisons drawn. Hence, the theory clearly indicates that having generated our posterior distribution, we can then query this distribution in any way that we wish thereby allowing us to explore all of our research questions simultaneously.
Bayesian “contrasts” can be performed either:
within the Bayesian sampling model or
construct them from the returned MCMC samples (they are drawn from the posteriors)
Only the latter will be demonstrated as it provides a consistent approach across all routines. In order to allow direct comparison to the frequentist equivalents, I will explore the same set of planned and Tukey’s test comparisons described here. For the “planned comparison” we defined two contrasts: 1) group B vs group C; and 2) group A vs the average of groups B and C. Of course each of these could be explored at multiple values of B, however, since we fit an additive model (which assumes that the slopes are homogeneous), the contrasts will be constant throughout the domain of B.
Lets start by comparing each group to each other group in a pairwise manner. Arguably the most elegant way to do this is to generate a Tukey’s contrast matrix. This is a model matrix specific to comparing each group to each other group. Again, since the lines are parallel, it does not really matter what level of B we estimate these efffects at - so lets use the mean B.
mcmc = data.rstan
coefs <- as.matrix(mcmc)[, 1:4]
newdata <- data.frame(A = levels(data$A), B = mean(data$B))
# A Tukeys contrast matrix
library(multcomp)
tuk.mat <- contrMat(n = table(newdata$A), type = "Tukey")
Xmat <- model.matrix(~A + B, data = newdata)
pairwise.mat <- tuk.mat %*% Xmat
pairwise.matNA (Intercept) AGroup B AGroup C B
NA Group B - Group A 0 1 0 0
NA Group C - Group A 0 0 1 0
NA Group C - Group B 0 -1 1 0mcmc_areas(coefs %*% t(pairwise.mat))
(comps = tidyMCMC(coefs %*% t(pairwise.mat), conf.int = TRUE, conf.method = "HPDinterval"))NA # A tibble: 3 × 5
NA term estimate std.error conf.low conf.high
NA <chr> <dbl> <dbl> <dbl> <dbl>
NA 1 Group B - Group A -16.2 1.61 -19.3 -13.0
NA 2 Group C - Group A -20.6 1.63 -24.0 -17.5
NA 3 Group C - Group B -4.38 1.57 -7.43 -1.36ggplot(comps, aes(y = estimate, x = term)) + geom_pointrange(aes(ymin = conf.low,
ymax = conf.high)) + geom_hline(yintercept = 0, linetype = "dashed") +
scale_y_continuous("Effect size") + scale_x_discrete("") + coord_flip() +
theme_classic()
With a couple of modifications, we could also express this as percentage changes. A percentage change represents the change (difference between groups) divided by one of the groups (determined by which group you want to express the percentage change to). Hence, we generate an additional mcmc matrix that represents the cell means for the divisor group (group we want to express change relative to). Since the tuk.mat defines comparisons as \(-1\) and \(1\) pairs, if we simply replace all the \(-1\) with \(0\), the eventual matrix multiplication will result in estimates of the divisor cell means instead of the difference. We can then divide the original mcmc matrix above with this new divisor mcmc matrix to yield a mcmc matrix of percentage change.
# Modify the tuk.mat to replace -1 with 0. This will allow us to get a
# mcmc matrix of ..
tuk.mat[tuk.mat == -1] = 0
comp.mat <- tuk.mat %*% Xmat
comp.matNA (Intercept) AGroup B AGroup C B
NA Group B - Group A 1 1 0 19.29344
NA Group C - Group A 1 0 1 19.29344
NA Group C - Group B 1 0 1 19.29344comp.mcmc = 100 * (coefs %*% t(pairwise.mat))/coefs %*% t(comp.mat)
(comps = tidyMCMC(comp.mcmc, conf.int = TRUE, conf.method = "HPDinterval"))NA # A tibble: 3 × 5
NA term estimate std.error conf.low conf.high
NA <chr> <dbl> <dbl> <dbl> <dbl>
NA 1 Group B - Group A -63.9 8.58 -79.4 -46.1
NA 2 Group C - Group A -98.2 12.1 -122. -76.4
NA 3 Group C - Group B -21.1 8.37 -38.7 -6.23ggplot(comps, aes(y = estimate, x = term)) + geom_pointrange(aes(ymin = conf.low,
ymax = conf.high)) + geom_hline(yintercept = 0, linetype = "dashed") +
scale_y_continuous("Effect size (%)") + scale_x_discrete("") + coord_flip() +
theme_classic()
And now for the specific planned comparisons (Group B vs Group C as well as Group A vs the average of Groups B and C). This is achieved by generating our own contrast matrix (defining the contributions of each group to each contrast).
c.mat = rbind(c(0, 1, -1), c(1/2, -1/3, -1/3))
c.matNA [,1] [,2] [,3]
NA [1,] 0.0 1.0000000 -1.0000000
NA [2,] 0.5 -0.3333333 -0.3333333mcmc = data.rstan
coefs <- as.matrix(mcmc)[, 1:4]
newdata <- data.frame(A = levels(data$A), B = mean(data$B))
Xmat <- model.matrix(~A + B, data = newdata)
c.mat = c.mat %*% Xmat
c.matNA (Intercept) AGroup B AGroup C B
NA [1,] 0.0000000 1.0000000 -1.0000000 0.000000
NA [2,] -0.1666667 -0.3333333 -0.3333333 -3.215574(comps = tidyMCMC(as.mcmc(coefs %*% t(c.mat)), conf.int = TRUE, conf.method = "HPDinterval"))NA # A tibble: 2 × 5
NA term estimate std.error conf.low conf.high
NA <chr> <dbl> <dbl> <dbl> <dbl>
NA 1 var1 4.38 1.57 1.36 7.43
NA 2 var2 5.31 0.795 3.74 6.88Variance components, the amount of added variance attributed to each influence, are traditionally estimated for so called random effects. These are the effects for which the levels employed in the design are randomly selected to represent a broader range of possible levels. For such effects, effect sizes (differences between each level and a reference level) are of little value. Instead, the “importance” of the variables are measured in units of variance components. On the other hand, regular variance components for fixed factors (those whose measured levels represent the only levels of interest) are not logical - since variance components estimate variance as if the levels are randomly selected from a larger population. Nevertheless, in order to compare and contrast the scale of variability of both fixed and random factors, it is necessary to measure both on the same scale (sample or population based variance).
Finite-population variance components assume that the levels of all factors (fixed and random) in the design are all the possible levels available (Gelman et al. (2005)). In other words, they are assumed to represent finite populations of levels. Sample (rather than population) statistics are then used to calculate these finite-population variances (or standard deviations). Since standard deviation (and variance) are bound at zero, standard deviation posteriors are typically non-normal. Consequently, medians and HPD intervals are more robust estimates.
mcmc = as.matrix(data.rstan)
head(mcmc)NA parameters
NA iterations beta[1] beta[2] beta[3] beta[4] sigma log_lik[1]
NA [1,] 48.68797 -14.90045 -18.42159 -0.4094795 3.031655 -2.158433
NA [2,] 48.43319 -15.01286 -18.67982 -0.4294815 2.921049 -2.233506
NA [3,] 50.43634 -16.82976 -19.72486 -0.4513570 3.537004 -2.185470
NA [4,] 48.64655 -13.89231 -18.70568 -0.4382509 2.484179 -2.128525
NA [5,] 48.15618 -14.92421 -17.71449 -0.4328598 2.243914 -2.276010
NA [6,] 49.17068 -14.51494 -16.74284 -0.4175619 2.329755 -1.888469
NA parameters
NA iterations log_lik[2] log_lik[3] log_lik[4] log_lik[5] log_lik[6] log_lik[7]
NA [1,] -2.156859 -2.275958 -2.508336 -2.427388 -2.108191 -2.140407
NA [2,] -2.044102 -2.612585 -2.770405 -2.669474 -2.323688 -2.282355
NA [3,] -2.451184 -2.376191 -2.359273 -2.323440 -2.258361 -2.208821
NA [4,] -1.913384 -2.778915 -2.889791 -2.755609 -2.364344 -2.240115
NA [5,] -1.765343 -3.072510 -3.312357 -3.126649 -2.513844 -2.397504
NA [6,] -2.092972 -2.134149 -2.415902 -2.297064 -1.876290 -1.891225
NA parameters
NA iterations log_lik[8] log_lik[9] log_lik[10] log_lik[11] log_lik[12]
NA [1,] -2.028076 -2.120660 -2.194535 -2.951822 -2.175850
NA [2,] -1.996090 -2.033820 -2.079288 -2.513460 -2.025906
NA [3,] -2.295370 -2.462257 -2.527712 -2.399819 -2.188735
NA [4,] -1.829593 -1.907385 -1.971099 -3.114748 -2.116560
NA [5,] -1.760407 -1.756790 -1.808722 -2.434182 -1.751882
NA [6,] -1.783484 -2.033781 -2.180475 -3.781417 -2.256792
NA parameters
NA iterations log_lik[13] log_lik[14] log_lik[15] log_lik[16] log_lik[17]
NA [1,] -3.923757 -5.518989 -3.165092 -4.103392 -2.371347
NA [2,] -3.465041 -6.540221 -3.558522 -3.382449 -2.140047
NA [3,] -3.040044 -5.490873 -3.280067 -2.833088 -2.237948
NA [4,] -4.864954 -6.591713 -3.087007 -4.551258 -2.414516
NA [5,] -3.957672 -9.900892 -4.622784 -3.761895 -1.898696
NA [6,] -5.730245 -6.697238 -3.093358 -5.858561 -2.683569
NA parameters
NA iterations log_lik[18] log_lik[19] log_lik[20] log_lik[21] log_lik[22]
NA [1,] -2.115685 -2.266963 -2.050612 -2.037393 -2.258567
NA [2,] -2.060764 -2.544361 -1.990902 -1.992742 -2.076796
NA [3,] -2.270484 -2.687233 -2.182458 -2.201401 -2.312137
NA [4,] -2.339028 -2.155251 -1.953336 -1.828934 -1.955807
NA [5,] -1.814116 -2.851421 -1.731678 -1.751657 -2.064323
NA [6,] -2.216344 -1.964871 -1.947308 -2.360778 -3.301713
NA parameters
NA iterations log_lik[23] log_lik[24] log_lik[25] log_lik[26] log_lik[27]
NA [1,] -2.171641 -2.830499 -2.395792 -3.395862 -6.222735
NA [2,] -2.051781 -3.196050 -2.196774 -3.713123 -5.390383
NA [3,] -2.322982 -2.751841 -2.452424 -2.937520 -4.652791
NA [4,] -1.939132 -3.416791 -2.142931 -4.009022 -6.408910
NA [5,] -2.011512 -3.220021 -2.361215 -3.934801 -8.395927
NA [6,] -3.059341 -2.047610 -3.778508 -2.492069 -12.270572
NA parameters
NA iterations log_lik[28] log_lik[29] log_lik[30] lp__
NA [1,] -2.029591 -2.071100 -2.240269 -52.19822
NA [2,] -2.065377 -1.996808 -2.018447 -51.76482
NA [3,] -2.199483 -2.240936 -2.207482 -50.07295
NA [4,] -1.933361 -1.847537 -1.851903 -55.18946
NA [5,] -1.753449 -1.825760 -1.880609 -59.66453
NA [6,] -2.067420 -2.608989 -3.110990 -65.84951# A
wch = grep("beta.[2-3]", colnames(mcmc))
# Get the rowwise standard deviations between effects parameters
sd.A = apply(mcmc[, wch], 1, sd)
# B
wch = grep("beta.[4]", colnames(mcmc))
sd.B = sd(data$B) * abs(mcmc[, wch])
# generate a model matrix
newdata = data
Xmat = model.matrix(~A + B, newdata)
## get median parameter estimates
wch = grep("beta", colnames(mcmc))
coefs = mcmc[, wch]
fit = coefs %*% t(Xmat)
resid = sweep(fit, 2, data$Y, "-")
sd.resid = apply(resid, 1, sd)
sd.all = cbind(sd.A, sd.B, sd.resid)
(fpsd = tidyMCMC(sd.all, conf.int = TRUE, conf.method = "HPDinterval"))NA # A tibble: 3 × 5
NA term estimate std.error conf.low conf.high
NA <chr> <dbl> <dbl> <dbl> <dbl>
NA 1 sd.A 3.10 1.10 0.963 5.24
NA 2 sd.B 7.11 0.684 5.84 8.46
NA 3 sd.resid 3.45 0.155 3.26 3.76# OR expressed as a percentage
(fpsd.p = tidyMCMC(100 * sd.all/rowSums(sd.all), estimate.method = "median",
conf.int = TRUE, conf.method = "HPDinterval"))NA # A tibble: 3 × 5
NA term estimate std.error conf.low conf.high
NA <chr> <dbl> <dbl> <dbl> <dbl>
NA 1 sd.A 22.2 6.21 8.99 32.9
NA 2 sd.B 52.3 4.37 43.9 60.6
NA 3 sd.resid 25.5 2.97 21.0 31.4## we can even plot this as a Bayesian ANOVA table
ggplot(fpsd, aes(y = estimate, x = term)) + geom_hline(yintercept = 0,
linetype = "dashed") + geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
geom_text(aes(label = sprintf("%.2f%%", fpsd.p$estimate), vjust = -1)) +
scale_y_continuous("Finite population standard deviation") + scale_x_discrete() +
coord_flip() + theme_classic()
Approximately \(22.86\)% of the total finite population standard deviation is due to \(x\).
In a frequentist context, the \(R^2\) value is seen as a useful indicator of goodness of fit. Whilst it has long been acknowledged that this measure is not appropriate for comparing models (for such purposes information criterion such as AIC are more appropriate), it is nevertheless useful for estimating the amount (percent) of variance explained by the model. In a frequentist context, \(R^2\) is calculated as the variance in predicted values divided by the variance in the observed (response) values. Unfortunately, this classical formulation does not translate simply into a Bayesian context since the equivalently calculated numerator can be larger than the an equivalently calculated denominator - thereby resulting in an \(R^2\) greater than \(100\)%. Gelman et al. (2019) proposed an alternative formulation in which the denominator comprises the sum of the explained variance and the variance of the residuals.
So in the standard regression model notation of:
\[ y_i \sim \text{Normal}(\boldsymbol X \boldsymbol \beta, \sigma), \]
the \(R^2\) could be formulated as
\[ R^2 = \frac{\sigma^2_f}{\sigma^2_f + \sigma^2_e}, \]
where \(\sigma^2_f=\text{var}(\boldsymbol X \boldsymbol \beta)\), and for normal models \(\sigma^2_e=\text{var}(y-\boldsymbol X \boldsymbol \beta)\).
mcmc <- as.matrix(data.rstan)
Xmat = model.matrix(~A + B, data)
wch = grep("beta", colnames(mcmc))
coefs = mcmc[, wch]
fit = coefs %*% t(Xmat)
resid = sweep(fit, 2, data$Y, "-")
var_f = apply(fit, 1, var)
var_e = apply(resid, 1, var)
R2 = var_f/(var_f + var_e)
tidyMCMC(as.mcmc(R2), conf.int = TRUE, conf.method = "HPDinterval")NA # A tibble: 1 × 5
NA term estimate std.error conf.low conf.high
NA <chr> <dbl> <dbl> <dbl> <dbl>
NA 1 var1 0.905 0.0147 0.875 0.921# for comparison with frequentist
summary(lm(Y ~ A + B, data))NA
NA Call:
NA lm(formula = Y ~ A + B, data = data)
NA
NA Residuals:
NA Min 1Q Median 3Q Max
NA -6.4381 -2.2244 -0.6829 2.1732 8.6607
NA
NA Coefficients:
NA Estimate Std. Error t value Pr(>|t|)
NA (Intercept) 51.00608 1.44814 35.22 < 2e-16 ***
NA AGroup B -16.25472 1.54125 -10.55 6.92e-11 ***
NA AGroup C -20.65596 1.57544 -13.11 5.74e-13 ***
NA B -0.48399 0.04526 -10.69 5.14e-11 ***
NA ---
NA Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
NA
NA Residual standard error: 3.44 on 26 degrees of freedom
NA Multiple R-squared: 0.9149, Adjusted R-squared: 0.9051
NA F-statistic: 93.22 on 3 and 26 DF, p-value: 4.901e-14Generate the data with heterogeneous slope effects.
set.seed(123)
n <- 10
p <- 3
A.eff <- c(40, -15, -20)
beta <- c(-0.45, -0.1, 0.5)
sigma <- 4
B <- rnorm(n * p, 0, 15)
A <- gl(p, n, lab = paste("Group", LETTERS[1:3]))
mm <- model.matrix(~A * B)
data1 <- data.frame(A = A, B = B, Y = as.numeric(c(A.eff, beta) %*% t(mm)) + rnorm(n * p, 0, 4))
data1$B <- data1$B + 20
head(data1)NA A B Y
NA 1 Group A 11.59287 45.48907
NA 2 Group A 16.54734 40.37341
NA 3 Group A 43.38062 33.05922
NA 4 Group A 21.05763 43.03660
NA 5 Group A 21.93932 42.41363
NA 6 Group A 45.72597 31.17787scatterplot(Y ~ B | A, data = data1)
boxplot(Y ~ A, data1)
# OR via ggplot
ggplot(data1, aes(y = Y, x = B, group = A)) + geom_point() + geom_smooth(method = "lm")
ggplot(data1, aes(y = Y, x = A)) + geom_boxplot()
The slopes (\(Y\) vs B trends) do appear to differ between treatment groups - in particular, Group C seems to portray a different trend to Groups A and B.
anova(lm(Y ~ B * A, data = data1))NA Analysis of Variance Table
NA
NA Response: Y
NA Df Sum Sq Mean Sq F value Pr(>F)
NA B 1 442.02 442.02 41.380 1.187e-06 ***
NA A 2 2760.60 1380.30 129.217 1.418e-13 ***
NA B:A 2 285.75 142.87 13.375 0.0001251 ***
NA Residuals 24 256.37 10.68
NA ---
NA Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There is strong evidence to suggest that the assumption of equal slopes is violated.
modelString2 = "
data {
int<lower=1> n;
int<lower=1> nX;
vector [n] y;
matrix [n,nX] X;
}
parameters {
vector[nX] beta;
real<lower=0> sigma;
}
transformed parameters {
vector[n] mu;
mu = X*beta;
}
model {
// Likelihood
y~normal(mu,sigma);
// Priors
beta ~ normal(0,100);
sigma~cauchy(0,5);
}
generated quantities {
vector[n] log_lik;
for (i in 1:n) {
log_lik[i] = normal_lpdf(y[i] | mu[i], sigma);
}
}
"
## write the model to a text file
writeLines(modelString2, con = "ancovaModel2.stan")Arrange the data as a list (as required by Stan). As input, Stan will need to be supplied with: the response variable, the predictor matrix, the number of predictors, the total number of observed items. This all needs to be contained within a list object. We will create two data lists, one for each of the hypotheses.
Xmat <- model.matrix(~A * B, data1)
data1.list <- with(data1, list(y = Y, X = Xmat, nX = ncol(Xmat), n = nrow(data1)))Define the nodes (parameters and derivatives) to monitor and the chain parameters.
params <- c("beta", "sigma", "log_lik")
nChains = 2
burnInSteps = 500
thinSteps = 1
numSavedSteps = 2000 #across all chains
nIter = ceiling(burnInSteps + (numSavedSteps * thinSteps)/nChains)
nIterNA [1] 1500Start the JAGS model (check the model, load data into the model, specify the number of chains and compile the model).
data1.rstan <- stan(data = data1.list, file = "ancovaModel2.stan", chains = nChains, pars = params, iter = nIter, warmup = burnInSteps, thin = thinSteps)NA
NA SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
NA Chain 1:
NA Chain 1: Gradient evaluation took 2.9e-05 seconds
NA Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.29 seconds.
NA Chain 1: Adjust your expectations accordingly!
NA Chain 1:
NA Chain 1:
NA Chain 1: Iteration: 1 / 1500 [ 0%] (Warmup)
NA Chain 1: Iteration: 150 / 1500 [ 10%] (Warmup)
NA Chain 1: Iteration: 300 / 1500 [ 20%] (Warmup)
NA Chain 1: Iteration: 450 / 1500 [ 30%] (Warmup)
NA Chain 1: Iteration: 501 / 1500 [ 33%] (Sampling)
NA Chain 1: Iteration: 650 / 1500 [ 43%] (Sampling)
NA Chain 1: Iteration: 800 / 1500 [ 53%] (Sampling)
NA Chain 1: Iteration: 950 / 1500 [ 63%] (Sampling)
NA Chain 1: Iteration: 1100 / 1500 [ 73%] (Sampling)
NA Chain 1: Iteration: 1250 / 1500 [ 83%] (Sampling)
NA Chain 1: Iteration: 1400 / 1500 [ 93%] (Sampling)
NA Chain 1: Iteration: 1500 / 1500 [100%] (Sampling)
NA Chain 1:
NA Chain 1: Elapsed Time: 0.065 seconds (Warm-up)
NA Chain 1: 0.045 seconds (Sampling)
NA Chain 1: 0.11 seconds (Total)
NA Chain 1:
NA
NA SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
NA Chain 2:
NA Chain 2: Gradient evaluation took 4e-06 seconds
NA Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
NA Chain 2: Adjust your expectations accordingly!
NA Chain 2:
NA Chain 2:
NA Chain 2: Iteration: 1 / 1500 [ 0%] (Warmup)
NA Chain 2: Iteration: 150 / 1500 [ 10%] (Warmup)
NA Chain 2: Iteration: 300 / 1500 [ 20%] (Warmup)
NA Chain 2: Iteration: 450 / 1500 [ 30%] (Warmup)
NA Chain 2: Iteration: 501 / 1500 [ 33%] (Sampling)
NA Chain 2: Iteration: 650 / 1500 [ 43%] (Sampling)
NA Chain 2: Iteration: 800 / 1500 [ 53%] (Sampling)
NA Chain 2: Iteration: 950 / 1500 [ 63%] (Sampling)
NA Chain 2: Iteration: 1100 / 1500 [ 73%] (Sampling)
NA Chain 2: Iteration: 1250 / 1500 [ 83%] (Sampling)
NA Chain 2: Iteration: 1400 / 1500 [ 93%] (Sampling)
NA Chain 2: Iteration: 1500 / 1500 [100%] (Sampling)
NA Chain 2:
NA Chain 2: Elapsed Time: 0.069 seconds (Warm-up)
NA Chain 2: 0.045 seconds (Sampling)
NA Chain 2: 0.114 seconds (Total)
NA Chain 2:print(data1.rstan, par = c("beta", "sigma"))NA Inference for Stan model: anon_model.
NA 2 chains, each with iter=1500; warmup=500; thin=1;
NA post-warmup draws per chain=1000, total post-warmup draws=2000.
NA
NA mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
NA beta[1] 48.04 0.10 1.91 44.19 46.84 48.04 49.29 51.73 375 1
NA beta[2] -10.41 0.12 2.74 -15.56 -12.18 -10.39 -8.72 -4.83 561 1
NA beta[3] -26.32 0.12 2.43 -31.11 -27.90 -26.35 -24.79 -21.39 430 1
NA beta[4] -0.34 0.00 0.08 -0.50 -0.40 -0.35 -0.29 -0.19 416 1
NA beta[5] -0.28 0.00 0.11 -0.49 -0.34 -0.27 -0.21 -0.08 483 1
NA beta[6] 0.26 0.00 0.11 0.04 0.19 0.26 0.33 0.48 511 1
NA sigma 3.36 0.02 0.50 2.55 3.01 3.29 3.65 4.50 1090 1
NA
NA Samples were drawn using NUTS(diag_e) at Mon Jul 22 11:50:52 2024.
NA For each parameter, n_eff is a crude measure of effective sample size,
NA and Rhat is the potential scale reduction factor on split chains (at
NA convergence, Rhat=1).mcmc <- As.mcmc.list(data1.rstan)
denplot(mcmc, parms = c("beta"))
traplot(mcmc, parms = c("beta"))
Trace plots show no evidence that the chains have not reasonably traversed the entire multidimensional parameter space. When there are a lot of parameters, this can result in a very large number of traceplots. To focus on just certain parameters (such as \(\beta\)s).
#Raftery diagnostic
raftery.diag(mcmc)NA [[1]]
NA
NA Quantile (q) = 0.025
NA Accuracy (r) = +/- 0.005
NA Probability (s) = 0.95
NA
NA You need a sample size of at least 3746 with these values of q, r and s
NA
NA [[2]]
NA
NA Quantile (q) = 0.025
NA Accuracy (r) = +/- 0.005
NA Probability (s) = 0.95
NA
NA You need a sample size of at least 3746 with these values of q, r and sThe Raftery diagnostics for each chain estimate that we would require no more than \(5000\) samples to reach the specified level of confidence in convergence. As we have \(10500\) samples, we can be confidence that convergence has occurred.
#Autocorrelation diagnostic
autocorr.diag(mcmc)NA beta[1] beta[2] beta[3] beta[4] beta[5] beta[6]
NA Lag 0 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000
NA Lag 1 0.47965891 0.40490097 0.42921394 0.47213590 0.43786125 0.40348216
NA Lag 5 0.10276601 0.04520477 0.08439589 0.08552624 0.07188554 0.05055728
NA Lag 10 0.03369686 0.04235285 0.03837580 0.04276107 0.06964443 0.05297790
NA Lag 50 -0.05533000 -0.03563501 -0.08740181 -0.04617438 -0.01873684 -0.08254539
NA sigma log_lik[1] log_lik[2] log_lik[3] log_lik[4]
NA Lag 0 1.000000000 1.00000000 1.000000000 1.00000000 1.000000000
NA Lag 1 0.244539989 0.30295557 0.115863972 0.31646347 0.029392467
NA Lag 5 0.034914225 0.03692325 0.025476065 0.04472861 0.009289407
NA Lag 10 0.027104230 0.04375764 0.001071547 0.03530329 0.036276640
NA Lag 50 -0.005837235 -0.04398093 -0.022416916 -0.01700789 -0.017665756
NA log_lik[5] log_lik[6] log_lik[7] log_lik[8] log_lik[9]
NA Lag 0 1.00000000 1.000000000 1.000000000 1.000000000 1.00000000
NA Lag 1 0.05686476 0.279972053 0.333411513 0.393178740 0.26907784
NA Lag 5 0.01335058 0.009737116 0.056417267 -0.001796254 0.02998479
NA Lag 10 0.03736092 0.035609098 0.008349122 0.056906159 0.01515109
NA Lag 50 -0.01570145 -0.012353687 -0.007107282 -0.056604829 -0.03352728
NA log_lik[10] log_lik[11] log_lik[12] log_lik[13] log_lik[14]
NA Lag 0 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000
NA Lag 1 0.192883022 0.068556254 0.204989891 0.003861059 0.105690603
NA Lag 5 0.037648652 0.031608422 0.030527908 0.011161059 0.040450947
NA Lag 10 0.007251782 -0.006369060 0.001890380 0.036208400 0.008267025
NA Lag 50 -0.028984961 -0.005015408 0.005481482 -0.018976796 -0.004975828
NA log_lik[15] log_lik[16] log_lik[17] log_lik[18] log_lik[19]
NA Lag 0 1.000000000 1.000000000 1.0000000000 1.00000000 1.00000000
NA Lag 1 -0.021918535 0.028936497 0.1274525386 0.03355231 -0.01865606
NA Lag 5 -0.006083026 0.025989718 0.0299094329 0.02020100 0.01961257
NA Lag 10 0.008300113 0.014088547 -0.0117960594 0.01141223 0.02524133
NA Lag 50 0.004556124 -0.009052317 -0.0004925802 0.02661374 -0.01229104
NA log_lik[20] log_lik[21] log_lik[22] log_lik[23] log_lik[24]
NA Lag 0 1.00000000 1.000000000 1.00000000 1.000000000 1.000000000
NA Lag 1 0.17819054 0.120052667 0.20632108 -0.024461686 -0.052472858
NA Lag 5 0.03745210 -0.020154436 0.03360932 -0.031301123 -0.009290672
NA Lag 10 0.03345325 0.004498098 0.02884422 -0.002109614 -0.014376687
NA Lag 50 -0.01491661 -0.016212205 -0.02709841 -0.017716486 0.004503930
NA log_lik[25] log_lik[26] log_lik[27] log_lik[28] log_lik[29]
NA Lag 0 1.000000000 1.000000000 1.00000000 1.000000000 1.00000000
NA Lag 1 -0.026932812 -0.010753209 -0.03500057 -0.022972322 0.02432922
NA Lag 5 -0.007038723 -0.018424596 -0.01393167 -0.006548756 -0.03234786
NA Lag 10 0.001166187 0.002621322 0.01265460 0.014571694 -0.00114465
NA Lag 50 -0.007274914 -0.014512013 0.01665142 -0.016130612 -0.01750531
NA log_lik[30] lp__
NA Lag 0 1.000000000 1.000000e+00
NA Lag 1 -0.015077656 4.530274e-01
NA Lag 5 0.006166306 -7.071898e-05
NA Lag 10 0.028726884 8.583829e-03
NA Lag 50 -0.014572075 -8.625792e-02stan_rhat(data1.rstan)
stan_ess(data1.rstan)
Rhat and effective sample size. In this instance, most of the parameters have reasonably high effective samples and thus there is likely to be a good range of values from which to estimate paramter properties.
mcmc = as.data.frame(data1.rstan) %>% dplyr:::select(contains("beta"),
sigma) %>% as.matrix
# generate a model matrix
newdata1 = data1
Xmat = model.matrix(~A * B, newdata1)
## get median parameter estimates
coefs = apply(mcmc[, 1:6], 2, median)
fit = as.vector(coefs %*% t(Xmat))
resid = data1$Y - fit
ggplot() + geom_point(data = NULL, aes(y = resid, x = fit)) + theme_classic()
Residuals against predictors
mcmc = as.data.frame(data1.rstan) %>% dplyr:::select(contains("beta"),
sigma) %>% as.matrix
# generate a model matrix
newdata1 = newdata1
Xmat = model.matrix(~A * B, newdata1)
## get median parameter estimates
coefs = apply(mcmc[, 1:6], 2, median)
fit = as.vector(coefs %*% t(Xmat))
resid = data1$Y - fit
newdata1 = newdata1 %>% cbind(fit, resid)
ggplot(newdata1) + geom_point(aes(y = resid, x = A)) + theme_classic()
ggplot(newdata1) + geom_point(aes(y = resid, x = B)) + theme_classic()
And now for studentised residuals
mcmc = as.data.frame(data1.rstan) %>% dplyr:::select(contains("beta"),
sigma) %>% as.matrix
# generate a model matrix
newdata1 = data1
Xmat = model.matrix(~A * B, newdata1)
## get median parameter estimates
coefs = apply(mcmc[, 1:6], 2, median)
fit = as.vector(coefs %*% t(Xmat))
resid = data1$Y - fit
sresid = resid/sd(resid)
ggplot() + geom_point(data = NULL, aes(y = sresid, x = fit)) + theme_classic()
For this simple model, the studentised residuals yield the same pattern as the raw residuals (or the Pearson residuals for that matter). Lets see how well data simulated from the model reflects the raw data.
mcmc = as.data.frame(data1.rstan) %>% dplyr:::select(contains("beta"),
sigma) %>% as.matrix
# generate a model matrix
Xmat = model.matrix(~A * B, data1)
## get median parameter estimates
coefs = mcmc[, 1:6]
fit = coefs %*% t(Xmat)
## draw samples from this model
yRep = sapply(1:nrow(mcmc), function(i) rnorm(nrow(data1), fit[i,
], mcmc[i, "sigma"]))
newdata1 = data.frame(A = data1$A, B = data1$B, yRep) %>% gather(key = Sample,
value = Value, -A, -B)
ggplot(newdata1) + geom_violin(aes(y = Value, x = A, fill = "Model"),
alpha = 0.5) + geom_violin(data = data1, aes(y = Y, x = A,
fill = "Obs"), alpha = 0.5) + geom_point(data = data1, aes(y = Y,
x = A), position = position_jitter(width = 0.1, height = 0),
color = "black") + theme_classic()
ggplot(newdata1) + geom_violin(aes(y = Value, x = B, fill = "Model",
group = B, color = A), alpha = 0.5) + geom_point(data = data1,
aes(y = Y, x = B, group = B, color = A)) + theme_classic()
The predicted trends do encapsulate the actual data, suggesting that the model is a reasonable representation of the underlying processes. Note, these are prediction intervals rather than confidence intervals as we are seeking intervals within which we can predict individual observations rather than means. We can also explore the posteriors of each parameter.
mcmc_intervals(as.matrix(data1.rstan), regex_pars = "beta|sigma")
mcmc_areas(as.matrix(data1.rstan), regex_pars = "beta|sigma")
First, we look at the results from the additive model.
print(data1.rstan, pars = c("beta", "sigma"))NA Inference for Stan model: anon_model.
NA 2 chains, each with iter=1500; warmup=500; thin=1;
NA post-warmup draws per chain=1000, total post-warmup draws=2000.
NA
NA mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
NA beta[1] 48.04 0.10 1.91 44.19 46.84 48.04 49.29 51.73 375 1
NA beta[2] -10.41 0.12 2.74 -15.56 -12.18 -10.39 -8.72 -4.83 561 1
NA beta[3] -26.32 0.12 2.43 -31.11 -27.90 -26.35 -24.79 -21.39 430 1
NA beta[4] -0.34 0.00 0.08 -0.50 -0.40 -0.35 -0.29 -0.19 416 1
NA beta[5] -0.28 0.00 0.11 -0.49 -0.34 -0.27 -0.21 -0.08 483 1
NA beta[6] 0.26 0.00 0.11 0.04 0.19 0.26 0.33 0.48 511 1
NA sigma 3.36 0.02 0.50 2.55 3.01 3.29 3.65 4.50 1090 1
NA
NA Samples were drawn using NUTS(diag_e) at Mon Jul 22 11:50:52 2024.
NA For each parameter, n_eff is a crude measure of effective sample size,
NA and Rhat is the potential scale reduction factor on split chains (at
NA convergence, Rhat=1).# OR
tidyMCMC(data1.rstan, conf.int = TRUE, conf.method = "HPDinterval", pars = c("beta", "sigma"))NA # A tibble: 7 × 5
NA term estimate std.error conf.low conf.high
NA <chr> <dbl> <dbl> <dbl> <dbl>
NA 1 beta[1] 48.0 1.91 44.2 51.8
NA 2 beta[2] -10.4 2.74 -16.0 -5.34
NA 3 beta[3] -26.3 2.43 -30.9 -21.3
NA 4 beta[4] -0.345 0.0767 -0.492 -0.193
NA 5 beta[5] -0.276 0.105 -0.487 -0.0735
NA 6 beta[6] 0.262 0.110 0.0393 0.479
NA 7 sigma 3.36 0.501 2.52 4.46Conclusions
The intercept of the first group (Group A) is \(48.2\).
The mean of the second group (Group B) is \(-10.6\) units greater than (A).
The mean of the third group (Group C) is \(-26.5\) units greater than (A).
A one unit increase in B in Group A is associated with a \(-0.351\) units increase in \(Y\).
difference in slope between Group B and Group A \(-0.270\).
difference in slope between Group C and Group A \(0.270\).
The \(95\)% confidence interval for the effects of Group B, Group C and the partial slope associated with B do not overlapp with \(0\) implying a significant difference between group A and groups B, C (at the mean level of predictor B) and a significant negative relationship with B (for Group A). The slope associated with Group B was not found to be significantly different from that associated with Group A, however, the slope associated with Group C was found to be significantly less negative than the slope associated with Group A. While workers attempt to become comfortable with a new statistical framework, it is only natural that they like to evaluate and comprehend new structures and output alongside more familiar concepts. One way to facilitate this is via Bayesian p-values that are somewhat analogous to the frequentist p-values for investigating the hypothesis that a parameter is equal to zero.
## since values are less than zero
mcmcpvalue(as.matrix(data1.rstan)[, "beta[2]"]) # effect of (B-A = 0)NA [1] 0.0015mcmcpvalue(as.matrix(data1.rstan)[, "beta[3]"]) # effect of (C-A = 0)NA [1] 0mcmcpvalue(as.matrix(data1.rstan)[, "beta[4]"]) # effect of (slope = 0)NA [1] 0mcmcpvalue(as.matrix(data1.rstan)[, "beta[5]"]) # effect of (slopeB - slopeA = 0)NA [1] 0.0095mcmcpvalue(as.matrix(data1.rstan)[, "beta[6]"]) # effect of (slopeC - slopeA = 0)NA [1] 0.028mcmcpvalue(as.matrix(data1.rstan)[, 2:6]) # effect of (model)NA [1] 0There is evidence that the response differs between the groups.
(full = loo(extract_log_lik(data1.rstan)))NA
NA Computed from 2000 by 30 log-likelihood matrix.
NA
NA Estimate SE
NA elpd_loo -83.0 4.8
NA p_loo 6.9 2.0
NA looic 166.0 9.5
NA ------
NA MCSE of elpd_loo is NA.
NA MCSE and ESS estimates assume independent draws (r_eff=1).
NA
NA Pareto k diagnostic values:
NA Count Pct. Min. ESS
NA (-Inf, 0.7] (good) 29 96.7% 174
NA (0.7, 1] (bad) 1 3.3% <NA>
NA (1, Inf) (very bad) 0 0.0% <NA>
NA See help('pareto-k-diagnostic') for details.# now fit a model without main factor
modelString3 = "
data {
int<lower=1> n;
int<lower=1> nX;
vector [n] y;
matrix [n,nX] X;
}
parameters {
vector[nX] beta;
real<lower=0> sigma;
}
transformed parameters {
vector[n] mu;
mu = X*beta;
}
model {
// Likelihood
y~normal(mu,sigma);
// Priors
beta ~ normal(0,1000);
sigma~cauchy(0,5);
}
generated quantities {
vector[n] log_lik;
for (i in 1:n) {
log_lik[i] = normal_lpdf(y[i] | mu[i], sigma);
}
}
"
## write the model to a stan file
writeLines(modelString3, con = "ancovaModel3.stan")
Xmat <- model.matrix(~A + B, data1)
data1.list <- with(data1, list(y = Y, X = Xmat, n = nrow(data1), nX = ncol(Xmat)))
data1.rstan.red <- stan(data = data1.list, file = "ancovaModel3.stan", chains = nChains,
iter = nIter, warmup = burnInSteps, thin = thinSteps)NA
NA SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
NA Chain 1:
NA Chain 1: Gradient evaluation took 1.8e-05 seconds
NA Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.18 seconds.
NA Chain 1: Adjust your expectations accordingly!
NA Chain 1:
NA Chain 1:
NA Chain 1: Iteration: 1 / 1500 [ 0%] (Warmup)
NA Chain 1: Iteration: 150 / 1500 [ 10%] (Warmup)
NA Chain 1: Iteration: 300 / 1500 [ 20%] (Warmup)
NA Chain 1: Iteration: 450 / 1500 [ 30%] (Warmup)
NA Chain 1: Iteration: 501 / 1500 [ 33%] (Sampling)
NA Chain 1: Iteration: 650 / 1500 [ 43%] (Sampling)
NA Chain 1: Iteration: 800 / 1500 [ 53%] (Sampling)
NA Chain 1: Iteration: 950 / 1500 [ 63%] (Sampling)
NA Chain 1: Iteration: 1100 / 1500 [ 73%] (Sampling)
NA Chain 1: Iteration: 1250 / 1500 [ 83%] (Sampling)
NA Chain 1: Iteration: 1400 / 1500 [ 93%] (Sampling)
NA Chain 1: Iteration: 1500 / 1500 [100%] (Sampling)
NA Chain 1:
NA Chain 1: Elapsed Time: 0.052 seconds (Warm-up)
NA Chain 1: 0.035 seconds (Sampling)
NA Chain 1: 0.087 seconds (Total)
NA Chain 1:
NA
NA SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
NA Chain 2:
NA Chain 2: Gradient evaluation took 5e-06 seconds
NA Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
NA Chain 2: Adjust your expectations accordingly!
NA Chain 2:
NA Chain 2:
NA Chain 2: Iteration: 1 / 1500 [ 0%] (Warmup)
NA Chain 2: Iteration: 150 / 1500 [ 10%] (Warmup)
NA Chain 2: Iteration: 300 / 1500 [ 20%] (Warmup)
NA Chain 2: Iteration: 450 / 1500 [ 30%] (Warmup)
NA Chain 2: Iteration: 501 / 1500 [ 33%] (Sampling)
NA Chain 2: Iteration: 650 / 1500 [ 43%] (Sampling)
NA Chain 2: Iteration: 800 / 1500 [ 53%] (Sampling)
NA Chain 2: Iteration: 950 / 1500 [ 63%] (Sampling)
NA Chain 2: Iteration: 1100 / 1500 [ 73%] (Sampling)
NA Chain 2: Iteration: 1250 / 1500 [ 83%] (Sampling)
NA Chain 2: Iteration: 1400 / 1500 [ 93%] (Sampling)
NA Chain 2: Iteration: 1500 / 1500 [100%] (Sampling)
NA Chain 2:
NA Chain 2: Elapsed Time: 0.043 seconds (Warm-up)
NA Chain 2: 0.031 seconds (Sampling)
NA Chain 2: 0.074 seconds (Total)
NA Chain 2:(reduced = loo(extract_log_lik(data1.rstan.red)))NA
NA Computed from 2000 by 30 log-likelihood matrix.
NA
NA Estimate SE
NA elpd_loo -92.3 4.8
NA p_loo 5.7 2.0
NA looic 184.6 9.7
NA ------
NA MCSE of elpd_loo is NA.
NA MCSE and ESS estimates assume independent draws (r_eff=1).
NA
NA Pareto k diagnostic values:
NA Count Pct. Min. ESS
NA (-Inf, 0.7] (good) 29 96.7% 364
NA (0.7, 1] (bad) 1 3.3% <NA>
NA (1, Inf) (very bad) 0 0.0% <NA>
NA See help('pareto-k-diagnostic') for details.par(mfrow = 1:2, mar = c(5, 3.8, 1, 0) + 0.1, las = 3)
plot(full, label_points = TRUE)
plot(reduced, label_points = TRUE)
The expected out-of-sample predictive accuracy is substantially lower for the model that includes \(x\). This might be used to suggest that the inferential evidence for a general effect of \(x\) on \(y\).
mcmc = as.matrix(data1.rstan)
## Calculate the fitted values
newdata1 = expand.grid(A = levels(data1$A), B = seq(min(data1$B), max(data1$B),
len = 100))
Xmat = model.matrix(~A * B, newdata1)
coefs = mcmc[, c("beta[1]", "beta[2]", "beta[3]", "beta[4]", "beta[5]",
"beta[6]")]
fit = coefs %*% t(Xmat)
newdata1 = newdata1 %>% cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval"))
ggplot(newdata1, aes(y = estimate, x = B, fill = A)) + geom_ribbon(aes(ymin = conf.low,
ymax = conf.high), alpha = 0.2) + geom_line() + scale_y_continuous("Y") +
scale_x_continuous("B") + theme_classic()
As this is simple single factor ANOVA, we can simple add the raw data to this figure. For more complex designs with additional predictors, it is necessary to plot partial residuals.
## Calculate partial residuals fitted values
fdata1 = rdata1 = data1
fMat = rMat = model.matrix(~A * B, fdata1)
fit = as.vector(apply(coefs, 2, median) %*% t(fMat))
resid = as.vector(data1$Y - apply(coefs, 2, median) %*% t(rMat))
rdata1 = rdata1 %>% mutate(partial.resid = resid + fit)
ggplot(newdata1, aes(y = estimate, x = B, fill = A)) + geom_point(data = rdata1,
aes(y = partial.resid, x = B, color = A)) + geom_ribbon(aes(ymin = conf.low,
ymax = conf.high), alpha = 0.2) + geom_line() + scale_y_continuous("Y") +
scale_x_continuous("B") + theme_classic()