# Using mixed effects modelling to estimate and compare variability

Say I observe two groups of 10 people, measuring some quantity 100 times in each person. There will presumably be some variability across these 100 measures in each person. Can I use mixed effects analysis to assess whether this within-person variability is, on average, different between the two groups? For example, using traditional statistics, I could compute the standard deviation (SD) within each person then submit these SDs to an anova comparing the groups, but I wonder if this two-stage process can be replaced by a single mixed effects model, consequently obtaining the various advantages of mixed effects modelling (shrinkage, accounting for different numbers of observations per person, etc) as well.

To be clear, here is R code depicting the scenario described above and the SD/anova-based approach:

set.seed(1)

group_A_base_sd = 1
group_B_base_sd = 2
within_group_sd_of_sds = .1

n_per_group = 10
obs_per_id = 100

temp = data.frame(
id = 1:(n_per_group*2)
, group = rep(c('A','B'))
)

#generate example data
library(plyr) #to avoid loops (for coding convenience only)
obs_data = ddply(
.data = temp
, .variables = .(id,group)
, .fun = function(x){
#generate a unique sd for this individual
# based on their group's sd plus some
# within-group variability
id_sd = ifelse(
x$group=='A' , rnorm( 1 , group_A_base_sd , within_group_sd_of_sds ) , rnorm( 1 , group_B_base_sd , within_group_sd_of_sds ) ) #generate data points with the above generated # variability to_return = data.frame( obs_num = 1:obs_per_id , measurement = rnorm(obs_per_id,0,id_sd) ) return(to_return) } ) #first step of an anova-based approach: # compute SDs within each Ss obs_sds = ddply( .data = obs_data , .variables = .(id,group) , .fun = function(x){ to_return = data.frame( obs_sd = sd(x$measurement)
)
}
)

#second step of an anova-based approach:
# compute the anova on the SDs
summary(
aov(
formula = obs_sd~group
, data = obs_sds
)
)


You can structure the model along the following lines. Let,

$j = 1, 2$ be the two groups and

$i$ index the individuals in the two groups.

$y_{ij} \sim N(\mu_j,\sigma_j^2)$ $\forall \ i, j$

$\sigma_j^2 \sim IG(v,1)$ $\forall \ j$

(Note: $IG(.)$ is the inverse gamma distribution.)

Priors

$\mu_j \sim N(\bar{\mu},\sigma_{\mu}^2)$

$v \sim IG(\bar{v},1)$ is the prior for $v$.

The above structure will let you shrink the error variances ($\sigma_j^2$) appropriately. You can then evaluate whether the within group variability is different by looking at the credible intervals associated with the group variabilities.

You can calculate SD for each person and its standard error (for example using bootstrap). Then you can use rmeta package to do analysis. I think you should use some transformation of SD for example log (or maybe better log of variance).