Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

More text about autocorrelations. #13

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

More text about autocorrelations. #13

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
60 changes: 34 additions & 26 deletions vignettes/references.bib
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,37 +121,45 @@ @ARTICLE{Yeatman2018AFQB
}

@ARTICLE{VanRij2019,
title = {Analyzing the {Time} {Course} of {Pupillometric} {Data}},
title = {Analyzing the Time Course of Pupillometric data},
author = {Van Rij, Jacolien and Hendriks, Petra and Van Rijn,
Hedderik and Baayen, R. Harald and Wood, Simon N.},
abstract = {This article provides a tutorial for analyzing pupillometric
data. Pupildilation has become increasingly popular in
psychological and psycholinguistic research as a measure to
trace language processing. However, there is no general consensus
about procedures to analyze the data, with most studies analyzing
extracted features from the pupil dilation data instead of analyzing
the pupil dilation trajectories directly. Recent studies have started
to apply nonlinear regression and other methods to analyze the pupil
dilation trajectories directly, utilizing all available information
in the continuously measured signal. This article applies a nonlinear
regression analysis, generalized additive mixed modeling, and
illustrates how to analyze the full-time course of the pupil dilation
signal. The regression analysis is particularly suited for analyzing
pupil dilation in the fields of psychological and psycholinguistic
research because generalized additive mixed models can include complex
nonlinear interactions for investigating the effects of properties of
stimuli (e.g., formant frequency) or participants (e.g., working memory
score) on the pupil dilation signal. To account for the variation due to
participants and items, nonlinear random effects can be included. However,
one of the challenges for analyzing time series data is dealing with the
autocorrelation in the residuals, which is rather extreme for the pupillary
signal. On the basis of simulations, we explain potential causes of this
extreme autocorrelation, and on the basis of the experimental data, we
show how to reduce their adverse effects, allowing a much more coherent
interpretation of pupillary data than possible with feature-based techniques.},
journal = "Trends in Hearing",
volume = 23,
pages = 1-22,
month = jan,
year = 2019,
}


@ARTICLE{Olszowy2019-ge,
title = "Accurate autocorrelation modeling substantially improves {fMRI}
reliability",
author = "Olszowy, Wiktor and Aston, John and Rua, Catarina and Williams,
Guy B",
journal = "Nat. Commun.",
publisher = "Springer Science and Business Media LLC",
volume = 10,
number = 1,
pages = 1220,
abstract = "Given the recent controversies in some neuroimaging statistical
methods, we compare the most frequently used functional Magnetic
Resonance Imaging (fMRI) analysis packages: AFNI, FSL and SPM,
with regard to temporal autocorrelation modeling. This process,
sometimes known as pre-whitening, is conducted in virtually all
task fMRI studies. Here, we employ eleven datasets containing 980
scans corresponding to different fMRI protocols and subject
populations. We found that autocorrelation modeling in AFNI,
although imperfect, performed much better than the
autocorrelation modeling of FSL and SPM. The presence of residual
autocorrelated noise in FSL and SPM leads to heavily confounded
first level results, particularly for low-frequency experimental
designs. SPM's alternative pre-whitening method, FAST, performed
better than SPM's default. The reliability of task fMRI studies
could be improved with more accurate autocorrelation modeling. We
recommend that fMRI analysis packages provide diagnostic plots to
make users aware of any pre-whitening problems.",
month = dec,
year = 2019,
language = "en"
}
96 changes: 75 additions & 21 deletions vignettes/tractable-autocorrelations.Rmd
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,15 +15,15 @@ knitr::opts_chunk$set(
)
```

This vignette demonstrates the use of an AR1 model within a GAM to account for
autocorrelation of the error terms of the model. We will use data from the
Sarica dataset [@Sarica2017] to demonstrate the how to include an autocorrelation
term and its impact on the model statistics.
This vignette demonstrates the use of an AR1 model within the GAM to account for
the spatial autocorrelation of the errors of the model. The ideas that we use closely follow the work of Van Rij and colleagues [@VanRij2019]. While their work used data that is different in many respects from our data (time-series of eye pupil size in psycholingustic experiments), there are also some similarities to the tract profile data that we are analyzing in `tractable`. For example, the signals tend to change slowly over time in their analysis, and tend to change rather slowly in space in our analysis. In both cases, this means that the models fail to capture some characteristics of the data, and that some inferences from the GAM models tend to be anti-conservative, unless a mitigation strategy and careful model checking are implemented. In particular, in both cases, the residuals of the model may exhibit substantial auto-correlation. They may also end up being centered not on zero. It is always worth checking the underlying assumption of normally-distributed residuals by plotting a so-called QQ plot. Finally, we can use formal model comparison methods to adjudicate between alternative models.

We start by loading the `tractable` library:
As an example of this approach, we will use data from the Sarica dataset [@Sarica2017] to demonstrate how to include an autocorrelation term in the GAM, and its impact on model statistics. We start by loading the `tractable` library, as well as the [itsadug](https://www.rdocumentation.org/packages/itsadug/versions/2.4.1/topics/itsadug) and [gratia](https://gavinsimpson.github.io/gratia/) libraries, which both provide functionality to assess GAMs fit by `mgcv` (our workhorse for GAM fitting).

```{r setup}
library(tractable)
library(itsadug)
library(gratia)
```

Next, we will use a function that is included in `tractable` to read this dataset
Expand All @@ -37,11 +37,15 @@ sarica$group <- factor(sarica$class)
sarica$subjectID <- unclass(factor(sarica$subjectID))
```

We will first fit a GAM model without any autocorrelation structure using the
`tractable_single_bundle` function. This model will use "group" and "age" to
predict FA, while smoothing over the tract nodes. We will also use the automated
procedure implemented in `tractable_single_bundle` to determine the ideal value
for `k`, a parameter used to determine the number of spline functions. The default behavior for `tractable_single_bundle` is to include the AR1 model. To avoid this, we set the parameter `autocor` to `FALSE`.
We will first fit a GAM model that does not account for autocorrelation structure
in the residuals using the `tractable_single_bundle` function. This model will use
"group" and "age" to account for the measured FA, while smoothing over the tract
nodes. We will also use the automated procedure implemented in
`tractable_single_bundle` to determine the ideal value for `k`, a parameter used
to determine the number of spline functions. The default behavior for
`tractable_single_bundle` is to include an AR1 model to account for
autocorrelations, as we will see below. But for now, to avoid this, we first set
the parameter `autocor` to `FALSE`.


```{r fit no ac model}
Expand All @@ -67,22 +71,46 @@ In addition, we see that there is a statistically significant effect of group
summary(gam_fit_cst_no_ac)
```

In this model, no steps were taken to account for autocorrelated residuals. We will use a few model diagnostics to evaluate the model. First, we will use the `gratia::appraise` function, which presents a few different visuals about the model:

```{r appraise gam_fit_cst_no_ac, echo=TRUE, fig.height=8, fig.width=8}
gratia::appraise(gam_fit_cst_no_ac)
```


The top left plot is a QQ plot, it shows the residuals as a function of their quantile. In a perfect world (or at least a world in which model assumptions are valid), these points should fall along the equality line. This plot is not terrible, but there are some deviations. Another plot to look at is the plot of the residuals as a function of the prediction. Here, we look both at the overall location and shape of the point cloud, as well as for any signs of clear structure. In this case, we see some signs of trouble: the whole cloud is shifted away from zero, and there are what appear as trails of points, suggesting that some points form patterns. The first of these issues is also apparent in the bottom left, where residuals are not zero centerd in the histogram of model residuals.

Another diagnostic plot that is going to be crucial as a diagnostic is the plot of the autocorrelation function of the residuals. A plot of this sort is provided by the `itsadug` library:

```{r plot acf residuals without ac, echo=TRUE, fig.height=6, fig.width=6}
itsadug::acf_resid(gam_fit_cst_no_ac)
```

The dashed blue lines indicate the 95% confidence interval for the auto-correlation
function of white noise. Here, we see that the auto-correlation at many of the lags
(and particularly in the immediate neighbor, lag-1) is substantially larger than
would be expected for a function with no autocorrelations. These autocorrelations
pose a danger to inference not only because of mis-specification of the model, but
also because we are going to under-estimate the standard error of the model in
this setting and this will result in false positives (this is what Van Rij et al.
elegantly refer to as "anti-conservative" models)

To account for potential spatial autocorrelation of FA values along the length
of the tract profile, we can incorporate a AR1 model into our GAM. Briefly, this
AR1 model is a linear model that estimates the amount of influence of the FA value of the preceding node on the FA value of the current node (see [@VanRij2019](https://journals.sagepub.com/doi/pdf/10.1177/2331216519832483) for an overview of accounting for autocorrelation using the `mgcv` package).
of the tract profile, we can incorporate an AR1 model into our GAM. Briefly, the
AR1 model is a linear model that estimates and accounts for the amount of influence
of the model residual FA value of each node on the residual FA value of
their neighbor node. This is somewhat akin to "pre-whitening" that fMRI researchers undertake, to account for temporal auto-correlations in the time-series measured with fMRI (see e.g. [@Olszowy2019-ge]).

The AR1 model takes a parameter $\rho$ to estimate autocorrelation effects. We
can pass our initial model into the function `itsadug::start_value_rho` to
automatically determine the value of $\rho$. We can also plot the autocorrelation
by setting `plot=TRUE` within that function (pictured below).
automatically determine the value of $\rho$.

```{r calculate rho}
rho = itsadug::start_value_rho(gam_fit_cst_no_ac)
itsadug::start_value_rho(gam_fit_cst_no_ac, plot=T)
```{r estimate rho, echo=TRUE, fig.height=6, fig.width=6}
rho1 <- itsadug::start_value_rho(gam_fit_cst_no_ac)
print(rho1)
```

By default, the `tractable_single_bundle` function will determine the value of
$\rho$ and incorporate the AR1 model into the GAM estimation.
By default, the `tractable_single_bundle` function empirically determines the value of $\rho$ based on the data and uses it to incorporate the AR1 model of the residuals into the GAM estimation.

```{r fit model including AR structure}
gam_fit_cst <- tractable::tractable_single_bundle(
Expand All @@ -92,17 +120,43 @@ gam_fit_cst <- tractable::tractable_single_bundle(
group_by = "group",
covariates = c("age", "group"),
dwi_metric = "fa",
k = "auto"
k = 16
)
```

Examining the summary of the resulting GAM fit object shows us that the inclusion
of the AR1 model changes the resulting statistics of our model. Although there
is still a statistically significant effect of group (p=0.044), the value of the
t-statistic on this term has changed from 6.243 to 2.015.
t-statistic on this term has changed from 6.243 to 2.015, suggesting that the model
has become substantially more conservative.

```{r summarise AR model}
summary(gam_fit_cst)
```

Here as well, we can appraise the model with gratia:

```{r appraise gam_fit_cst, echo=TRUE, fig.height=8, fig.width=8}
gratia::appraise(gam_fit_cst)
```

Notice some improvements to model characteristics: the residuals are more centered around zero and the QQ plot is somewhat improved. There is some residual structure in the scatter plot of residuals as a function of prediction. We can ask how bad this structure is in terms of the residual autocorrelation:

```{r plot ACF with AR1, echo=TRUE, fig.height=6, fig.width=6}
rho2 <- itsadug::acf_resid(gam_fit_cst)
print(rho2)
```

This shows that the lag-1 autocorrelation has been reduced from approximately 0.94 to approximately 0.29.

Finally, formal model comparison can tell us which of these models better fit the
data. Using the `itsadug` library this can be done using the Akaike Information
Criterion as a comparator. In this case, this also indicates that the model that
accounts for autocorrelations also has smaller residuals considering the number of
parameters, suggesting that it is overall a better model of the data.

```{r compare models}
itsadug::compareML(gam_fit_cst_no_ac, gam_fit_cst)
```

## References
Loading