Multilevel Model with Heterogeneous Variance (Location-Scale Models)

Introduction to Concepts and Model (using nlme)

Overview

In the previous tutorials we covered how the multilevel model is used to examine intraindividual covariability.

In this tutorial, we outline how an extension, the multilevel model with heterogeneous variance (sometimes called location-scale models) can be used to examine differences in intraindividual variability - which we had previously done in a 2-step way using the iSD.

Prelim - Loading libraries used in this script

library(psych)    #data descriptives
library(ggplot2)  #plotting
library(nlme)     #multilevel models
library(brms)     #bayesian multilevel models
library(dplyr)    #data manipulation

Introduction

We have used two “separate” sets of methods to examine…
1. Intraindividual Variation (calculation of within-person summaries; iSD, iEntropy, iMSSD, etc. following Ram & Gerstorf, 2009)
2. Intraindividual Covariation (multilevel models - following Bolger & Laurenceau, 2013)

In 1., our interest was in the interindividual differences in within-person variability in an outcome variable, $\sigma_{yi}^2$.

In 2., our interest was in the interindividual differences in the within-person association between a predictor and outcome $\beta_{1i}$.

Here we outline a framework where these two interests can be combined together into a single model: Multilevel Model with Heterogeneous Variance.

Like with multilevel modeling itself, the framework is referred to by different names in different fields. The larger class of models includes models referred to as Variance Heterogeneity models, Location-Scale models, Generalized Additive Models for Location, Scale and Shape (GAMLSS, http://www.gamlss.org). It is not totally clear if there are some nuanced differences among the various different models. However, as far as I can tell, all of these model variants are applied to longitudinal panel-type data (multiple persons, multiple occasions).

A parallel form of the model exists for time-series data, ARCH models (auto-regressive models with conditional heteroscedasticity; https://en.wikipedia.org/wiki/Autoregressive_conditional_heteroskedasticity). There are lots of variants there including GARCH, NGARCH, IGARCH, EGARCH, GARCH-M, QGARCH, TGARCH. All of these model variants are applied to time-series data (single entity, many many occasions). As far as I can tell, there have not been many attempts to link the two different traditions, but there are a few.

Much of what we explore today is covered in Hedeker, Mermelstein, Berbaum, & Campbell (2009), Hedeker & Mermelstein (2007), and Hoffman (2007). Some more recent extensions are in Hedeker, Mermelstein & Demirtas (2012), and Rast, Hofer, & Sparks (2012).

Some Basic Definitions

A few definitions from Wikipedia (which may or may not be useful today, but are part of our larger effort to understand how we can make use of all these models).

Location parameter: https://en.wikipedia.org/wiki/Location_parameter
In statistics, a location family is a class of probability distributions that is parameterized by a scalar- or vector-valued parameter $x_0$, which determines the “location” or shift of the distribution. Examples of location parameters include the mean, the median, and the mode.

Scale parameters: https://en.wikipedia.org/wiki/Scale_parameter
In probability theory and statistics, a scale parameter is a special kind of numerical parameter of a parametric family of probability distributions. The larger the scale parameter, the more spread out the distribution.

The normal distribution has two parameters: a location parameter $\mu$ and a scale parameter $\sigma$. In practice the normal distribution is often parameterized in terms of the squared scale $\sigma^2$, which corresponds to the variance of the distribution.

Here you can see location and scale of the distribution changing over time (image from http://www.gamlss.com)

Shape parameters: https://en.wikipedia.org/wiki/Shape_parameter
A shape parameter is any parameter of a probability distribution that is neither a location parameter nor a scale parameter (nor a function of either or both of these only, such as a rate parameter). Such a parameter must affect the shape of a distribution rather than simply shifting it (as a location parameter does) or stretching/shrinking it (as a scale parameter does). Distributions with a shape parameter include skew normal, gamma, Weibull, ExGaussian distributions.

The Usual (Homogeneous Variance) Multilevel Model

Typically, the multilevel models we use (and that are covered in B&L) make a homogeneity of variance assumption. For example, lets look at the basic “unconditional means” model.

\[y_{it} = \beta_{0i} + e_{it}\] \[\beta_{0i} = \gamma_{00} + u_{0i}\]

where \[\mathbf{e_{it}} \sim \mathcal{N}(0,\mathbf{R})\], and \[\mathbf{u_{i}} \sim \mathcal{N}(0,\mathbf{G})\].

For clarity, lets write out the full variance covariance matrix of the between-person residuals, $\mathbf{G}$. In this case it is just a 1 x 1 matrix. \[\mathbf{G} = \left[\begin{array} {r} \sigma^{2}_{u0} \end{array}\right]\] Embedded here is the assumption that the between-person variance characterizes the entire sample (population). There are no subgroups. Later, we will relax this assumption, and allow that $\sigma^{2}_{u0g}$ can differ across “groups” of individuals (e.g., Reds, Greens).

For clarity, lets also write out the full variance covariance matrix of the within-person residuals, $\mathbf{R}$ (spanning across the t = 8 days of observation in the example data), \[\mathbf{R} = \mathbf{I} \left[\begin{array} {r} \sigma^{2}_{e} \end{array}\right] = \left[\begin{array} {rrrrrrrr} \sigma^{2}_{e} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \sigma^{2}_{e} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \sigma^{2}_{e} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \sigma^{2}_{e} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \sigma^{2}_{e} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \sigma^{2}_{e} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & \sigma^{2}_{e} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \sigma^{2}_{e} \end{array}\right]\], The homoscedasticity of errors applies across persons and occasions. There are no subgroups. Later, we will relax this assumption, and allow that $\sigma^{2}_{eg}$ can differ across g “groups” of individuals (e.g., Reds, Greens).

Here is a little plot of some simulated data that follow this model.

##         vars    n  mean    sd median trimmed   mad   min    max range skew
## x_i        1 2000 50.50 28.87  50.50   50.50 37.06  1.00 100.00 99.00 0.00
## group_i    2 2000  0.50  0.50   0.50    0.50  0.74  0.00   1.00  1.00 0.00
## b0_i       3 2000  0.12  1.26   0.00    0.04  0.01 -3.13   4.41  7.54 0.80
## sig2_ei    4 2000  0.15  0.05   0.15    0.15  0.07  0.10   0.20  0.10 0.00
## y_ti       5 2000  0.13  1.32   0.01    0.05  0.60 -3.92   5.02  8.94 0.69
## time_ti    6 2000 10.50  5.77  10.50   10.50  7.41  1.00  20.00 19.00 0.00
## id         7 2000 50.50 28.87  50.50   50.50 37.06  1.00 100.00 99.00 0.00
##         kurtosis   se
## x_i        -1.20 0.65
## group_i    -2.00 0.01
## b0_i        2.02 0.03
## sig2_ei    -2.00 0.00
## y_ti        1.83 0.03
## time_ti    -1.21 0.13
## id         -1.20 0.65

From this plot we can see …
* between-person differences in $u_{0i}$,
* homogeneity of within-person variance $\sigma^{2}_{e}$

The Heterogeneous Variance (Location-Scale) Multilevel Model

The objective of this model is to relax the homogeneity of variance assumption, and we can do this at both the between-person level, allowing for group or individual differences in $\sigma^{2}_{u0g}$, and/or at the within-person level, allowing for group or individual differences in $\sigma^{2}_{eg}$. We expand the basic “unconditional means” model to illustrate.

\[y_{it} = \beta_{0i} + e_{it}\] \[\beta_{0i} = \gamma_{00} + u_{0i}\] where \[\mathbf{e_{it}} \sim \mathcal{N}(0,\mathbf{R})\], and \[\mathbf{u_{i}} \sim \mathcal{N}(0,\mathbf{G})\]. That is all the same as usual.

But now lets add some complexity. Rather than the previous assumption that $\sigma^{2}_{u0g} = \sigma^{2}_{u0}$, lets allow that $\sigma^{2}_{u0g} = f(g_{i})$, where $g_{i}$ is a grouping variable (Reds, Greens).
This allows that $\sigma^{2}_{u0g=1} \neq \sigma^{2}_{u0g=2}$.

Formally, the full variance covariance matrix of the between-person residuals is now,
\[\mathbf{G} = \left[\begin{array} {rr} \sigma^{2}_{u0g=1} & 0 \\ 0 & \sigma^{2}_{u0g=2} \end{array}\right]\] There are now, explicitly, two variances $\sigma^{2}_{u0g=1}$ and $\sigma^{2}_{u0g=2}$, with the off-diagonal covariance being zero because individuals are only in one group. Lets look at some example data.

Note how the extent of between-person differences in $u_{0i}$ is different in the g = 1 = Red and g = 2 = Green groups. Visually, the difference in the group variances map to
\[\mathbf{G} = \left[\begin{array} {rr} \sigma^{2}_{u0g=1}=0 & 0 \\ 0 & \sigma^{2}_{u0g=2}=1 \end{array}\right]\]

A Series of Analysis Models

Let’s set up an analysis model. We use the same lme() function from the nlme package that we have been using, and simply make some adjustments (some are similar to what we did for the dyadic models).

Specific for the heterogeneity in between-person differences we make adjustments to the random = formula part of the script.

Homogeneity of Variances - The usual model

This is a model of “common” between-person variance and “common” within person variance, where “common” means that these variances are equivalent across groups, $\sigma^{2}_{u0g} = \sigma^{2}_{u0}$.

# Common models for between and within person variance
model_01_lme = lme(fixed = y_ti ~ 1,  
                   random = ~ 1 | id,
                           data = sim_nogrow,
                       method = 'REML')
summary(model_01_lme)

## Linear mixed-effects model fit by REML
##   Data: sim_nogrow 
##        AIC     BIC    logLik
##   2351.889 2368.69 -1172.944
## 
## Random effects:
##  Formula: ~1 | id
##         (Intercept)  Residual
## StdDev:    1.271246 0.3797654
## 
## Fixed effects:  y_ti ~ 1 
##                 Value Std.Error   DF   t-value p-value
## (Intercept) 0.1258579 0.1274079 1900 0.9878346  0.3234
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.489280944 -0.612660509  0.004013223  0.621295100  4.750794513 
## 
## Number of Observations: 2000
## Number of Groups: 100

Now, without changing the model, we are going to write the random = formula as a list(). The exact same model is … where pdSymm() indicates that the random effects are structured as a positive-definite Symmetric matrix.

# Common models for between and within person variance
model_01_lme = lme(fixed = y_ti ~ 1,  
                   random = list(id = pdSymm(form = ~ 1)),
                           data = sim_nogrow,
                       method = 'REML')
summary(model_01_lme)

## Linear mixed-effects model fit by REML
##   Data: sim_nogrow 
##        AIC     BIC    logLik
##   2351.889 2368.69 -1172.944
## 
## Random effects:
##  Formula: ~1 | id
##         (Intercept)  Residual
## StdDev:    1.271246 0.3797654
## 
## Fixed effects:  y_ti ~ 1 
##                 Value Std.Error   DF   t-value p-value
## (Intercept) 0.1258579 0.1274079 1900 0.9878346  0.3234
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.489280944 -0.612660509  0.004013223  0.621295100  4.750794513 
## 
## Number of Observations: 2000
## Number of Groups: 100

Looking specifically at the variance structures …

VarCorr(model_01_lme)

## id = pdSymm(1) 
##             Variance  StdDev   
## (Intercept) 1.6160662 1.2712459
## Residual    0.1442217 0.3797654

We see that $\sigma^{2}_{u0g} = \sigma^{2}_{u0} = 1.6160662$,
and that $\sigma^{2}_{eg} = \sigma^{2}_{e} = 0.1442217$.

The true values match the descriptives above. Yay!

Heterogeneity in Between-person residuals (Level-2 random effects).

This is a model of heterogeneous between-person variance and common within person variance. Note that in addition to changing the form = we also explicitly specify that the $\mathbf{G}$ matrix should be diagonal with the pdDiag() specification (positive-definite Diagonal). This gets us the set-up we want (see above), \[\mathbf{G} = \left[\begin{array} {rr} \sigma^{2}_{u0g=0} & 0 \\ 0 & \sigma^{2}_{u0g=1} \end{array}\right]\]

# Heterogeneity for between variance and common within person variance
model_01b_lme = lme(fixed = y_ti ~ 1,  
                    random = list(id = pdDiag(form = ~ factor(group_i))),
                            data = sim_nogrow,
                        method = 'REML')
summary(model_01b_lme)

## Linear mixed-effects model fit by REML
##   Data: sim_nogrow 
##        AIC      BIC    logLik
##   2104.694 2127.096 -1048.347
## 
## Random effects:
##  Formula: ~factor(group_i) | id
##  Structure: Diagonal
##         (Intercept) factor(group_i)1  Residual
## StdDev: 2.44783e-05         1.798439 0.3780477
## 
## Fixed effects:  y_ti ~ 1 
##                     Value  Std.Error   DF     t-value p-value
## (Intercept) -0.0002883882 0.01194176 1900 -0.02414955  0.9807
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.525619386 -0.631708514  0.003199987  0.624580177  4.787642278 
## 
## Number of Observations: 2000
## Number of Groups: 100

Looking specifically at the variance and “correlation” structures …

VarCorr(model_01b_lme)

## id = pdDiag(factor(group_i)) 
##                  Variance     StdDev      
## (Intercept)      5.991871e-10 0.0000244783
## factor(group_i)1 3.234381e+00 1.7984385740
## Residual         1.429200e-01 0.3780476785

We see that …
$\sigma^{2}_{u0g=0} = 0.00$, $\sigma^{2}_{u0g=1} = 3.23$,
and that $\sigma^{2}_{eg} = \sigma^{2}_{e} = 0.149$.

Ok - so hopefully the idea of variance heterogeneity in the between-person random effects is clear.

Lets turn to the within-person random effects.

Heterogeneity in within-person residual variance (Level-1 random effects).

This is a model of heterogeneous between-person variance and heterogeneous within person variance. We do this using the weights = formula in the same way that we did for the dyadic models, in order that we can get the desired structure … \[\mathbf{R} = \left[\begin{array} {rr} \sigma^{2}_{eg=0} & 0 \\ 0 & \sigma^{2}_{eg=1} \end{array}\right]\]

Again note the diagonal structure, meaning that we need the weights =part of the code, but we do not need the correlation = part of the code we use in the dyadic models.

Our expanded model is specified as …

# Heterogeneity for both between and within person variance
model_01c_lme = lme(fixed = y_ti ~ 1,  
                    random = list(id = pdDiag(form = ~ factor(group_i))),
                    weights = varIdent(form = ~ 1 | factor(group_i)),
                            data = sim_nogrow,
                        method = 'REML')
summary(model_01c_lme)

## Linear mixed-effects model fit by REML
##   Data: sim_nogrow 
##        AIC      BIC   logLik
##   2005.756 2033.758 -997.878
## 
## Random effects:
##  Formula: ~factor(group_i) | id
##  Structure: Diagonal
##          (Intercept) factor(group_i)1  Residual
## StdDev: 4.500022e-05         1.797808 0.3145904
## 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | factor(group_i) 
##  Parameter estimates:
##       0       1 
## 1.00000 1.38244 
## Fixed effects:  y_ti ~ 1 
##                     Value   Std.Error   DF     t-value p-value
## (Intercept) -0.0004595558 0.009940638 1900 -0.04623001  0.9631
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.058631812 -0.651976600  0.003777448  0.649732254  4.157940316 
## 
## Number of Observations: 2000
## Number of Groups: 100

The between-person variances are given by the upper portion of …

VarCorr(model_01c_lme)

## id = pdDiag(factor(group_i)) 
##                  Variance     StdDev      
## (Intercept)      2.025020e-09 4.500022e-05
## factor(group_i)1 3.232113e+00 1.797808e+00
## Residual         9.896710e-02 3.145904e-01

For the within-person variances we need to do some multiplication (and squaring to get variances). The Residual StdDev is here ..

summary(model_01c_lme)$sigma

## [1] 0.3145904

and the weights are here … (yeah, I’m not sure why they are so hard to extract)

coef(model_01c_lme$modelStruct$varStruct, unconstrained=FALSE)

##       1 
## 1.38244

So the estimated Residual variance for group = 0, $\sigma^{2}_{eg=0}$ is … (true value = .10)

(summary(model_01c_lme)$sigma*1.0000)^2

## [1] 0.0989671

and the estimated variance for group = 1, $\sigma^{2}_{eg=1}$ is … (true value = .20)

(summary(model_01c_lme)$sigma*coef(model_01c_lme$modelStruct$varStruct, uncons=FALSE))^2

##       1 
## 0.18914

Yay! - now we expanded the original model to include both heterogeneous between-person variances and heterogeneous within-person variances.

To formally test if they heterogeneous variances are needed we would do a log-likelihood test, e.g., using the anova() function.

anova(model_01_lme, model_01b_lme)

##               Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## model_01_lme      1  3 2351.889 2368.690 -1172.944                        
## model_01b_lme     2  4 2104.694 2127.096 -1048.347 1 vs 2 249.1944  <.0001

anova(model_01b_lme, model_01c_lme)

##               Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## model_01b_lme     1  4 2104.694 2127.096 -1048.347                        
## model_01c_lme     2  5 2005.756 2033.758  -997.878 1 vs 2 100.9382  <.0001

Those tests indicate that the more complex models are significantly better for these data.

Lets calculate the within-person variances for the two groups

#for group = 0
exp(-1.1581 + 0*0.3268)^2

## [1] 0.09864774

#for group = 1
exp(-1.1581 + 1*0.3268)^2

## [1] 0.1896453

True values are 0.10 and 0.20, so the model did pretty good!

Continuous Predictor of Heterogeneous Within-person Variance

In the above we accommodated group (i.e., categorical) differences in between-person and within-person variances. We can also model the variances as a function of a continuous variable. For example, \[\sigma^{2}_{ei} = f(x_{i})\], where $x_{i}$ is an individual differences variable (e.g., sex, neuroticism).

Lets look at what some data might look like.

Formally, the model is \[y_{it} = \beta_{0i} + e_{it}\] \[\beta_{0i} = \gamma_{00} + \gamma_{01}x_i + u_{0i}\]

where \[\mathbf{e_{it}} \sim \mathcal{N}(0,\mathbf{R})\], and \[\mathbf{u_{i}} \sim \mathcal{N}(0,\mathbf{G})\].

We assume the between-person residuals in intercepts are homogenous, $\mathbf{G}$. In this case it is just a 1 x 1 matrix. \[\mathbf{G} = \left[\begin{array} {r} \sigma^{2}_{u0} \end{array}\right]\] and the within-person residuals are heterogeneous, so that
\[\mathbf{R} = \mathbf{I} \left[\begin{array} {r} \sigma^{2}_{ei} \end{array}\right]\] where, here in lme() … see ?varExp \[\sigma^{2}_{ei} = \alpha_{0}^{2}exp(2\alpha_{1}x_i)\] (Yes, the details of the exact implementations are confusing, and not very clear in any book or paper. One has to dig a bit to discover that in lme() there is a multiplication by 2. A good, practical overview piece is needed. I’ve found one paper that is … ok.)

Our expanded model changes the specification in the weights = option to …

# Continuous predictor for within person variance, common between-person differences
model_02a_lme = lme(fixed = y_ti ~ 1 + x_i,  
                    random = list(id = pdSymm(form = ~ 1)),
                    weights = varExp(form = ~ x_i),
                            data = sim_nogrow[which(sim_nogrow$id <=100),],
                        method = 'REML')
summary(model_02a_lme)

## Linear mixed-effects model fit by REML
##   Data: sim_nogrow[which(sim_nogrow$id <= 100), ] 
##        AIC      BIC logLik
##   41399.99 41427.99 -20695
## 
## Random effects:
##  Formula: ~1 | id
##         (Intercept)  Residual
## StdDev:  0.05305307 0.2976798
## 
## Variance function:
##  Structure: Exponential of variance covariate
##  Formula: ~x_i 
##  Parameter estimates:
##     expon 
## 0.2007523 
## Fixed effects:  y_ti ~ 1 + x_i 
##                Value  Std.Error   DF  t-value p-value
## (Intercept) 0.042113 0.08369139 1900   0.5032  0.6149
## x_i         4.007326 0.01993610   98 201.0085  0.0000
##  Correlation: 
##     (Intr)
## x_i -0.783
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.263324084 -0.642496633 -0.004136157  0.665231737  4.327043763 
## 
## Number of Observations: 2000
## Number of Groups: 100

The residual between-person variances are given by the upper portion of …

VarCorr(model_02a_lme)

## id = pdSymm(1) 
##             Variance    StdDev    
## (Intercept) 0.002814628 0.05305307
## Residual    0.088613246 0.29767977

Specifically, residual variance in intercepts unrelated to $x_i* are … (true value = 0)

VarCorr(model_02a_lme)[1]

## [1] "0.002814628"

For the within-person variances we need to do some multiplication (and squaring to get variances). The Residual Std Deviation, $\alpha_{0}$ is here … (true value for $\alpha_{0}^2$ = 0.1)

summary(model_02a_lme)$sigma

## [1] 0.2976798

And the parameter in the exponential variance function is here … (true value = 0.2)

coef(model_02a_lme$modelStruct$varStruct, uncons=FALSE)

##     expon 
## 0.2007523

So for a given value of x_i we expect the residual within-person variance to be … \[\sigma^{2}_{ei} = \alpha_{0}^{2}exp(2\alpha_{1}x_i)\]

x_i <- as.matrix(1:50)
varpred <- summary(model_02a_lme)$sigma^2 * exp(coef(2*model_02a_lme$modelStruct$varStruct, uncons=FALSE)*x_i)
head(varpred,10)

##            [,1]
##  [1,] 0.1323945
##  [2,] 0.1978068
##  [3,] 0.2955374
##  [4,] 0.4415538
##  [5,] 0.6597128
##  [6,] 0.9856578
##  [7,] 1.4726428
##  [8,] 2.2002330
##  [9,] 3.2873044
## [10,] 4.9114664

And how does that compare to the true values in our simulation?

head(personframe, 10)

##    x_i group_i b0_i   sig2_ei
## 1    1       0    4 0.1290658
## 2    2       0    8 0.2300880
## 3    3       0   12 0.3781546
## 4    4       0   16 0.3737871
## 5    5       0   20 0.7779523
## 6    6       0   24 1.1713321
## 7    7       0   28 1.5348707
## 8    8       0   32 2.2974943
## 9    9       0   36 3.4201372
## 10  10       0   40 4.9067429

A bit underestimated, but pretty close.

We can make a plot of how the residual variances are related to the predictor

#make a data frame with predictions and log predictions
plotvars <- as.data.frame(cbind(x_i, varpred))
plotvars <- plotvars %>%
  rename(x_i = V1,
         varpred = V2) %>%
  mutate(logvarpred = log(varpred))

#plotting
plotvars %>%
  filter(x_i <= 10) %>%
  ggplot(aes(x=x_i,y = varpred), legend=FALSE) +
  geom_point() +
  geom_line() + 
  xlab("x") + 
  ylab("Predicted IIV (residual variance)") + #ylim(-5,5) +
  scale_x_continuous(breaks=seq(0, 10, by=1)) +
  theme_minimal()

Side note: Exponential function

Why did we do this with an exponential function? Same as Poisson, the predicted scores should not go below zero (we cannot have negative variances) So we transform to obtain a linear equation. \[log(\sigma^{2}_{ei}) = \alpha_{0}^{2} + 2\alpha_{1}x_i\] In the log space we get …

plotvars %>%
  filter(x_i <= 10) %>%
  ggplot(aes(x=x_i,y = logvarpred), legend=FALSE) +
  #geom_rect(mapping=aes(xmin=day-.5, xmax=day+.5, ymin=1, ymax=7, fill=stress), alpha=0.8) +
  geom_point() +
  geom_line() + 
  xlab("x") + 
  ylab("Predicted IIV (residual variance)") + #ylim(-5,5) +
  scale_x_continuous(breaks=seq(0, 10, by=1)) +
  theme_minimal()

Interim Thoughts

Now have the possibility to relax the homogeneity of variance assumption. And this leads us into a world where we may be able to do some more precise testing of how between-person difference characteristics are related to IIV.

We have to be careful with how we construct the models, and although lme() does not make it super easy to test for significance we can probably figure all the details out. Notably, lme() actually allows quite a few different variance functions (varFunc), so there is a rich set of potential models to be explored. The most detailed information I have found so far on how to do this in R is Chapter 7 of Galecki, A. & Burzykowski, T. (2013). Linear mixed-effects models using R. A step-by-step approach. Springer, New York. A few papers have come out in the last few years where these models are applied to experience sampling data - but there are certainly more papers that need to be written to bring the model to “mainstream” use.

Of note: There are reasons to do this in a Bayesian Framework - and we have explored the models using the brms package. It all works very nicely, is quite easy to code, and provides for additional precision through inclusion of covariances and residuals. It’s very cool! - and covered in another tutorial.

Conclusion

Ok - we have attempted to push into the variance heterogeneity models. These models are not yet widely used in the analysis of experience sampling data. There are still a variety of details that must be worked out and understood before publishing analyses with these models. For example, we still have to figure out how to display the raw data and the results in ways that facilitate intuitive understanding.
But, it is all within reach! And getting closer as computation and implementation get easier and easier!

Be Careful! It is easy to get into trouble as the models get more complicated! Make sure to do some simulations so that you are confident in your interpretation of the parameters. It is not always 2*straightforward. =:-]

Citations

Autoregressive conditional heteroskedasticity. (2025). In Wikipedia. https://en.wikipedia.org/w/index.php?title=Autoregressive_conditional_heteroskedasticity&oldid=1298090386

Bürkner, P.-C. (2017). brms: An R Package for Bayesian Multilevel Models Using Stan. Journal of Statistical Software, 80, 1–28. https://doi.org/10.18637/jss.v080.i01

Bolger, N., & Laurenceau, J.-P. (2013). Intensive longitudinal methods: An introduction to diary and experience sampling research (pp. xv, 256). The Guilford Press.

Gałecki, A., & Burzykowski, T. (2013). Linear Models with Heterogeneous Variance. In A. Gałecki & T. Burzykowski (Eds.), Linear Mixed-Effects Models Using R: A Step-by-Step Approach (pp. 123–147). Springer. https://doi.org/10.1007/978-1-4614-3900-4_7

Gamlss—For statistical modelling. (2023, June 20). https://www.gamlss.com/

Hedeker, D., & Mermelstein, R. J. (2007). Mixed-effects regression models with heterogeneous variance: Analyzing ecological momentary assessment (EMA) data of smoking. In Modeling contextual effects in longitudinal studies (pp. 183–206). Lawrence Erlbaum Associates Publishers.

Hedeker, D., Mermelstein, R. J., Berbaum, M. L., & Campbell, R. T. (2009). Modeling mood variation associated with smoking: An application of a heterogeneous mixed-effects model for analysis of ecological momentary assessment (EMA) data. Addiction (Abingdon, England), 104(2), 297–307. https://doi.org/10.1111/j.1360-0443.2008.02435.x

Hedeker, D., Mermelstein, R. J., & Demirtas, H. (2012). Modeling between-subject and within-subject variances in ecological momentary assessment data using mixed-effects location scale models. Statistics in Medicine, 31(27), 3328–3336. https://doi.org/10.1002/sim.5338

Hoffman, L. (2007). Multilevel Models for Examining Individual Differences in Within-Person Variation and Covariation Over Time. Multivariate Behavioral Research, 42(4), 609–629. https://doi.org/10.1080/00273170701710072

Location parameter. (2025). In Wikipedia. https://en.wikipedia.org/w/index.php?title=Location_parameter&oldid=1295021758

Pinheiro, J. C., Bates, D., & R Core Team. (2025). nlme: Linear and Nonlinear Mixed Effects Models (Version 3.1-167) [Computer software]. https://doi.org/10.1007/b98882

Ram, N., & Gerstorf, D. (2009). Time-structured and net intraindividual variability: Tools for examining the development of dynamic characteristics and processes. Psychology and Aging, 24(4), 778–791. https://doi.org/10.1037/a0017915

Revelle, W. (2024). psych: Procedures for Psychological, Psychometric, and Personality Research. Northwestern University. https://CRAN.R-project.org/package=psych

Scale parameter. (2025). In Wikipedia. https://en.wikipedia.org/w/index.php?title=Scale_parameter&oldid=1281036107

Shape parameter. (2023). In Wikipedia. https://en.wikipedia.org/w/index.php?title=Shape_parameter&oldid=1172404465

Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag.

Wickham, H., François, R., Henry, L., Müller, K., & Vaughan, D. (2023). dplyr: A Grammar of Data Manipulation (Version 1.1.4). https://CRAN.R-project.org/package=dplyr