Test for significance of unexpected non-linear trend

Question

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.

Answer 1

(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.