Report non-linear term in a full regression model

Question

Simple question about statistical reporting best practices. I was asked to report the full model as a supplementary material of my paper. With the help of a toy example, let's say I have two variables in the model: var1 (linear) and var2 (modelled as a restricted cubic spline with 4 knots).

Here is the output:

  Coef    S.E.    Wald Z Pr(>|Z|)

 var1   -0.8754  0.0158  8.79  <0.0001 
 var2    4.5077  1.5089  2.99  0.0028  
 var2'  -0.1702  4.8386 -0.04  0.9719  
 var2'' -5.4636 12.5472 -0.44  0.6632 

So the question is: how am I suppose to report this model in full, and particularly how to deal with non-linear terms?

My guess was to produce a table like this one:

Variable	Coefficient	SE	p
Var1	-0.8754	0.0158	<0.0001
Var2	4.5077	1.5089	0.0028
Var2 (non-linear term 1)	-0.1702	4.8386	0.9719
Var2 (non-linear term 2)	-5.4636	12.5472	0.6632

But I am not sure this would be appropriate (especially referring to "non-linear term 1" and "non-linear term 2").

Can anyone help?

Answer 1

I would keep the two sets of predictors separate as they convey different messages.

In Var1 and Var2 perhaps you may want to test if the relation between response and predictors is linear and significant.

The nonlinear spline terms, which are typically included for model flexibility, have more parameters and thus require specific explanation; for instance, here it's not obvious what Coefficient, se and $p$-values are referring to.

Answer 2

The individual coefficients for the spline can be difficult to interpret, and depend on the locations specified by the knots. The table that you propose might be included in supplemental material, but it's not going to be very helpful for presenting your results to your audience.

You need to document the magnitude and significance of the predictor as a whole and its combined non-linear terms.

For describing the magnitude of the association of var2 with outcome, a plot of the predicted outcome values over a range of var2 values, along with confidence intervals, is probably most informative. In your model, without interactions, you could choose any single value of var1 to document that. That provides a visual display of any deviation from linearity.

For describing the significance of the association of var2 with outcome overall and with respect to its non-linear terms, a set of "chunk" Wald tests on combinations of coefficients is a good choice. That takes into account both the individual variances of the coefficient estimates and their covariances.

Overall, you test whether any of the 3 coefficients associated with var2 is significantly different from zero. You similarly test whether either of the 2 non-linear coefficients (var2', var2'') is different from zero.

As these seem to be the types of coefficients reported by the rcs() function in Harrell's rms package, you can can do both of the above simply if you work with models produced by that package. You can use its Predict() function (capital "P") to get predictions for outcome as a function of var2, and plot with either base R graphics or with ggplot(). The anova() function applied to models from that package directly reports overall significance and significance of non-linear or interaction terms involving a predictor.