9
$\begingroup$

I used the following R code:

ggscatter(DF_w, x = "Age", y = "Value", 
          add = "loess", conf.int = TRUE 
          ) 

to get this plot: loess

However, I don't understand how I can report this in a scientific paper or presentation. I understand the graph suggests a non-linear relationship between X Age and Y Value, where there is some variability until late life when Value declines with increased Age. Is that all there is to it? Is there a way to show this is a significant relationship with loess or do I need to carry out a regular regression / correlation to determine significance?

Any information on interpreting such plots is appreciated.

$\endgroup$
10
  • 4
    $\begingroup$ Loess is for exploration and smoothing. It's a bit of a stretch to use it to demonstrate "significance," because what would the null hypothesis be? $\endgroup$
    – whuber
    Dec 16, 2021 at 20:02
  • 1
    $\begingroup$ @whuber: one could do something ANOVA-like and compare the loess fit to an intercept-only model, or to a simple linear regression. Perhaps some kind of permutation test could be used? $\endgroup$ Dec 16, 2021 at 20:15
  • $\begingroup$ @Stephan It's unclear how you would compute the DF for an ANOVA-like solution. Permutation tests and bootstrapping are always available, of course--but they still require formulating a test statistic. I wonder what statistic someone would select before examining the data? I also wonder whether and to what extent this approach might be superior to standard nonparametric trend tests, if at all. $\endgroup$
    – whuber
    Dec 16, 2021 at 21:31
  • $\begingroup$ As whuber has pointed out, you can use it for exploring apparent trends in the data. In this case, your response variable appears more variable for smaller values of Age. Additionally, a cubic or quintic trend is shown $\endgroup$
    – user277126
    Dec 16, 2021 at 22:17
  • $\begingroup$ @whuber do you have any suggestions on how to further explore this data? I am trying to show how value is changing with Age but a normal regression doesn't seem appropriate $\endgroup$
    – CanyonView
    Dec 17, 2021 at 0:50

2 Answers 2

12
$\begingroup$

Bootstrapping or a permutation test (as suggested by Stephan Kolassa) will help you assess the significance of the apparent (but complex) association in the plot.

You need to adopt a reasonable measure of what the Loess fit has accomplished. One is the mean squared residual. (Many others are possible, but a discussion of that is not pertinent here.) Let's call this the "loss."

You also need to formulate a specific null hypothesis. The simplest is that the data exhibit no trend: that is, they vary randomly and independently around a common value. The best way to estimate this common value with mean squared error loss is to compute its arithmetic mean. With a sufficient amount of data--often ten or more observations suffice--the differences between the response ("Value" in the plot) and this mean will be only slightly (negatively) correlated and can stand in as surrogates for the true random errors.

The permutation distribution in this setting is the distribution of the losses associated with all possible permutations (reorderings) of the residuals after conducting a Loess fit. Under the null hypothesis, all these permutations are equally likely.

The permutation test compares the actual loss to the permutation distribution of losses. As a practical matter, the latter is estimated by means of a few random permutations. (There are far too many permutations to permit generating them all.)


Here, to illustrate, are data generated with no inherent trend along with a permutation distribution estimated from 500 draws. The vertical red line shows the loss for these data: it is close to the middle. Its p-value is computed as usual for a two-sided test: the red line splits the histogram left and right into two areas and the p-value is twice the smaller area. Very small p-values are called "significant" and taken as evidence of some kind of trend. The shape of the Loess plot (shown at left) helps you interpret just what that trend might be.

Figure 1

The large, anodyne p-value is consistent with the trend-free method of generating these data.

For data closer to those in the question, the result is different:

Figure 2

The actual mean squared error of $0.000735$ is inconsistent with the squared errors typical of the permutation distribution: this is a significant trend. The data plot at the left suggests the trend is principally a decrease in mean values from $0.31$ at age 60 to $0.25$ at ages 80 and over.

BTW, it is no surprise that the actual statistics in both cases have approximately the same values: they both estimate the error variance, which was equal to $0.025^2 = 0.000625$ in both cases. The curvilinear trends in the second instance, though, cause the simple fit (under the null hypothesis) to be poorer, thereby shifting the permutation distribution to higher values, as you can see by comparing the two figures.


The R code needed is simple, clear, and efficient. fit performs the Loess fit while stat uses that to compute the mean squared error.

fit <- function(y, x, ...) lowess(x, y, ...)
stat <- function(y, x) mean((y -  fit(y, x)$y)^2) # Mean squared error loss

Given a data frame object X with Value and Age columns to store the response and explanatory variables, respectively, the permutation distribution is estimated by computing predicted values and residuals under the null hypothesis and then iteratively permuting the residuals (with the sample function) and recomputing the loss.

  predicted <- mean(X$Value)
  residuals <- X$Value - predicted
  dsample <- replicate(5e3, stat(predicted + sample(residuals), X$Age))

In this case, after about one second of computation, dsample winds up with 5e3 ($5000$) values randomly drawn from the permutation distribution. The figures are then created by applying hist to dsample to show these values.

  # Compute the p-value
  actual <- with(X, stat(Value, Age))
  stats <- c(actual, dsample)
  p <- mean(stats <= actual)
  p <- 2 * min(1/2, p, 1-p)
  # Display the results
  hist(dsample, freq=FALSE, xlim=range(stats),
       col=gray(.95),
       sub=paste("p-value is approximately", signif(p, 2)),
       main="Simulated Null Permutation Distribution",
       xlab = "Mean Squared Difference")
  abline(v = actual, lwd=2, col="Red") # The statistic for the data

One caveat: I had to tune this Loess fit by specifying a relatively short width for its search radius. This frequently is the case. An honest permutation test must implement this fine-tuning in some automatic way and apply it to each permutation. Otherwise, the original (hand-tuned) fit will be too good and the resulting p-value will be too small--perhaps far too small. People often use some sort of cross-validation technique for such automatic tuning.

$\endgroup$
8
  • $\begingroup$ So under the null of no trend, the sample mean is taken. Then the residuals are found, and these are “permuted” by sampling from them all without replacement, and adding them back to the mean. These are then run through loess and the MSE found, and then this is compared with the actual loess MSE. Is this a correct summary? There is a very slight readability reliance on knowing that the sample function with all values and without replacement is a permutation. $\endgroup$ Dec 17, 2021 at 18:56
  • 1
    $\begingroup$ @Single Aware that many readers will not know R well, I explained the code, including saying explicitly what sample does. Your summary is accurate, although a bit oversimplified (which is fair for any very short summary). $\endgroup$
    – whuber
    Dec 17, 2021 at 18:59
  • 1
    $\begingroup$ I read the text then went through the code neglecting to then remember the commentary that function sample was used to do the permutation. A permutation test is probably an underrated tool—after viewing this question I had not considered it as a solution to it. Possibly and speculatively, its under use is because it requires some tailoring to the context, and maybe off the shelf functions to perform it are not so common. $\endgroup$ Dec 17, 2021 at 19:07
  • 1
    $\begingroup$ @Single Those are interesting thoughts (I mean that truly, not as dismissive hyperbole). But since bootstrapping can sort of be automated, as in the boot function in R, permutation tests could be automated just as readily, one would think. As always, though, there is no substitute for thinking about the problem... . $\endgroup$
    – whuber
    Dec 17, 2021 at 19:28
  • 2
    $\begingroup$ @SingleMalt: I sort of agree. Permutation tests for simple situations are easy enough (e.g. for the equivalent of a t test for the means in two groups), and indeed some such are implemented in R and other software. But more complex situations may require a lot of thought, indeed, as to the test statistic and which permutation scheme represents the null hypothesis best. One could argue that having to think about this is a Good Thing, compared to blindly applying a test without understanding these factors. Incidentally, I liked Good's Permutation, Parametric, and Bootstrap Tests of Hypotheses. $\endgroup$ Dec 18, 2021 at 16:46
3
$\begingroup$

First off, it's extremely good practice to also show the original data, which puts the loess plot into context. Here, the context is that there is still a lot of variation in the data. For instance, the initial dip looks rather strange and could easily be due to noise - and we only see this because we see the full point cloud, not only the loess line with the confidence band. So this is good.

You ask two questions: one about interpretation, the other about significance. In terms of interpretation, I would discuss the downward slope at the right end, but per above, not really trust the dip in the middle.

In terms of significance, this is harder. You could try to assess whether your loess model explains significantly more variation in the data than a comparison model, like an intercept-only model (a horizontal flat line), or a simple linear regression (a slanted straight line). This is what ANOVA does. The problem is that the standard F test in ANOVA requires knowing how many parameters (degrees of freedom) your model used - and that is notoriously difficult to know in the case of a loess model. Greg Snow's answer to How do I find a p-value of smooth spline / loess regression? gives a few tentative ideas in his penultimate paragraph (though not in terms of variance explained, but his ideas could be adapted to this test statistic).

However, as Greg also notes in his answer, assessing the significance of a spline fit can be done in an ANOVA framework. Given the suspicious behavior of your loess fit, I would suggest that you try a fit using natural or restricted cubic splines with few knots, then test this spline fit against the more parsimonious linear fit.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.