# Test for significance of unexpected non-linear trend

How can I test whether there is a significant non-linear effect of a continuous predictor on a continuous DV?

Here's my situation. I ran a simple linear regression and found a significant effect of a predictor (predictor2 below). However, when I plot my dv against that predictor, the trend does not look linear - in fact it looks like the dv has two local maxima. When I submitted these findings to a journal, reviewers commented that I should explore the possible statistical significance of this apparent non-linear effect. The effect was not expected and so I do not have a particular non-linear model to test here. I am just looking for an appropriate off-the-shelf technique of first call to use in a situation like this. Ideas?

For reference, my data look like this:

> D[1:10,]
sid      posttest      pretest       predictor2      predictor3       predictor4        predictor5
1    2          0.75         0.75        0.8823529      1.00000000        0.8888889                18
2    3          1.00         1.00        0.8666667      0.06666667        1.0000000                16
3    5          1.00         1.00        0.9000000      1.00000000        0.9047619                21
4    6          0.75         0.75        0.3333333      0.93333333        0.9375000                16
5    7          0.50         0.75        0.8750000      0.87500000        0.8333333                18
6    8          0.75         1.00        0.5000000      0.50000000        0.8823529                17
7    9          1.00         0.50        0.8823529      0.58823529        0.4444444                18
8   10          1.00         0.50        0.8750000      0.87500000        0.5882353                17
9   13          0.75         0.75        0.7777778      1.00000000        0.8947368                19
10  15          0.50         0.25        0.8666667      0.46666667        0.8125000                16
...


The linear regression model and its output look like this:

> lmA <- lm( posttest ~ pretest + predictor2 + predictor3 + predictor4 + predictor5, data=D )
> summary( lmA )

Call:
lm(formula = posttest ~ pretest + predictor2 + predictor3 + predictor4 + predictor5, data = D)

Residuals:
Min       1Q   Median       3Q      Max
-0.68810 -0.11853  0.05159  0.13589  0.36904

Coefficients:
Estimate  Std. Error   t value  Pr(>|t|)
(Intercept)     0.773645    0.106033     7.296  8.45e-13 ***
pretest         0.292674    0.030905     9.470   < 2e-16 ***
predictor2      0.063541    0.024557     2.588   0.00988 **
predictor3      0.005529    0.026369     0.210   0.83398
predictor4      -0.026421   0.065049    -0.406   0.68474
predictor5      -0.009412   0.005010    -1.879   0.06070 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1801 on 665 degrees of freedom
Multiple R-squared:  0.1298,    Adjusted R-squared:  0.1232
F-statistic: 19.83 on 5 and 665 DF,  p-value: < 2.2e-16


Finally, here is the plot that gave rise to the reviewer comment. Note that binned data with bubbles were used, with bubble diameter indicating number of data rows in each bin, instead of a scatterplot because there are only 5 levels of the dv (0.0, .25, .5, .75, 1.0) so you can't see anything at all with a scatterplot. Reviewers commented that bins could be arbitrarily chosen to make the linear trend look more persuasive, which is true, but binned data still has to be used for the above reason I think.

• Why are you fitting linear regression with a non-continuous dependent variable? In any case the deviation from linear looks very weak to me. You might do a bootstrap test or permutation test for a local linear regression over a straight line, perhaps, or you might perhaps test for dependence in the residuals (say by a runs test). Commented Nov 26, 2013 at 5:50
• @Glen_b: With over 600 data points, the discreteness of the response should not cause much trouble as long as variance homogeneity (among other assumptions) looks acceptable. Commented Nov 26, 2013 at 6:41
• @MichaelMayer There's not just discreteness; the variable is bounded. Clearly the boundedness has the potential to impact linearity (it can't remain linear as you approach the bounds), equality of variance (variance must shrink if the mean approaches one of the bounds) and anything approaching symmetry (it will tend to become skew as you approach bounds). This seems made for GLMs. Commented Nov 26, 2013 at 7:31
• @Glen_b: That's indeed true. Starting from scratch, I would also try to go for an appropriate GLM with S-shaped link function. But given the review status of the manuscript, a quick and dirty solution also has its advantages. Commented Nov 26, 2013 at 8:04
• I was actually rejected and am going to resubmit to another journal, so I have leeway to make bigger changes if needed. :) "An appropriate GLM" would be something like glmer in R's lme4 package with family='binomial', using the correct/incorrect data for each trial (instead of aggregate data by subject) as my dv, right? But even so, that model wouldn't be able to accommodate a non-monotonic relationship, would it? How then could I do that in the context of a GLM? Commented Nov 26, 2013 at 15:58

## 1 Answer

(1) Do you want to use linear regression model with an outcome variable restricted to 5 levels?

Assuming you want to stick with the linear regression model you have an abundance of options, including but not limited to:
(2) Create new variable: 5 level categorical variable of predictor2. Add this variable to the model. If predictor2_categorical is significant this represents a significant deviation from linear trend. (just need to be careful how your stats program deals with the collinearity, it should remove on level from categorical variable for collinearity and keep continuous variable, and this only tests deviation from linear trend. coefficients aren't interpretable)
(3) Restricted cubic spline with 5 knots of predictor2. The first of the four spline variables will be the untransformed predictor. So can test departure from linearity by testing the statistical significance of other three spline variables together.
(4) The MARS/Earth package can be used to look for deviations from linearity (using GCV).
(5) GAM will give you sort of a statistical significance of the smoothing term.
(6) --as mentioned in the comments -- add square and cubic term and test for joint significance. This is likely the easiest to implement.

Your choice will likely depend on what is easiest to implement in your statistical package and how concerned you are about potential criticism from reviewers for using data driven approach.

(7) Number 6 was focused on checking for deviation from linearity. I mostly use Stata and this is incredible easy to do in that program:

gen pred2_2 = predictor2^3
gen pred2_3 = predictor2^3
quietly regress outcome predictor2 pred2_2 pred2_3
test pred2_2  presd2_3


I’m fairly weak with R. My best guess is it would be something like this:

data$pred2_2 = predictor2^2 data$pred2_3 = predictor2^3
attach(data)
full.model=lm(outcome ~ predictor2+pred2_2+pred2_3)
nested.model=lm(outcome ~ predictor2)
anova(full.model, nested.model)


### or using -car- package: results should be identical

library(car)
crPlots(full.model)
crPlots(nested.model)
linearHypothesis full.model, "pred2_2+pred2_3")


The -car- package has the appeal that it can be used with most models. Not convinced this is the right syntax though - should be something like this.

(8) I think there is some confusion about glms here. There is a concern that your outcome looks like a proportion and thus models for proportions might fit the data better (beta models etc.). This has been well covered before: http://www.stata.com/meeting/germany10/germany10_buis.pdf Implementation of log models for proportions here:http://www.ats.ucla.edu/stat/stata/faq/proportion.htm
Recent stackexchange thread: Generalized linear models with continuous proportions

(8) I’ve no idea about glmer. I know lme4 is a random effects package. But this is a whole different question from the question of the distribution of the outcome variable. This is: how should change scores be analyzed?

There are three ways
(a) as change score – subtract pre from post to get outcome
(b) repeat measures models such as GEE or mixed-effects models, and
(c) using pre value as predictor.

You’ve chosen (c) above. This is admittedly the least attractive of the options and is generally avoided in observational studies. In randomized studies, where there should be no dependence between the effect of the predictor and the baseline value of the outcome, this may be more reasonable.

This is all based on my understanding that you only have two readings on your subjects, if you have more then that changes things.

(9) Most models here (except for GAM I think) are linear in the predictors (even if they have different link functions), so the methods of allowing for non-linearity are fairly broadly applicable.

• Why not simply add a squared and cubic version of the IV and then check the p value of the corresponding partial F test? The referee's comment is data driven anyway. Commented Nov 26, 2013 at 6:37
• that is most probably the best/easiest approach. was getting carried away with all the options available. Commented Nov 26, 2013 at 14:59
• For (6) - could you give me some clues or refs about how to implement this in R? I have no familiarity with this. For (1) well ... yes, I want to, but I could be persuaded otherwise - alternatives were proposed in the comments to the OP. Would a GLM like that described in those comments work together with (6)? Commented Nov 26, 2013 at 16:01
• To (6): Compare the current model by 'anova' with the model adding the squared and the cubed regressor in question. But better go with a GLM. Commented Nov 26, 2013 at 19:35
• OK ... in that case, any idea how to do it with a GLM? As far as I am aware, anova cannot be used to compare the output of glmer from lme4 with family='binomial', the model I floated in comments to the OP. -- also, I was mainly asking how to implement the model itself in (6), not the model comparison, sorry for not being clear. Commented Nov 26, 2013 at 20:27