I have following model using mtcars dataset:

> library(rms)
> modr = ols(mpg~rcs(wt), mtcars)
> modr

Linear Regression Model

ols(formula = mpg ~ rcs(wt), data = mtcars)

                Model Likelihood     Discrimination    
                   Ratio Test           Indexes        
Obs       32    LR chi2     57.75    R2       0.835    
sigma 2.6196    d.f.            4    R2 adj   0.811    
d.f.      27    Pr(> chi2) 0.0000    g        6.210    


    Min      1Q  Median      3Q     Max 
-4.6560 -1.7268 -0.3086  1.3239  6.6445 

          Coef     S.E.     t     Pr(>|t|)
Intercept  48.6788   5.0076  9.72 <0.0001 
wt        -10.4203   2.3532 -4.43 0.0001  
wt'        19.3800  12.2277  1.58 0.1246  
wt''      -68.7537  44.6474 -1.54 0.1352  
wt'''     231.1378 145.3974  1.59 0.1235  

The model is created using rms package in R and using rcs function in regression equation. How can I plot or what kind of plot can be made to show above model with all 4 components of wt shown in the output above? The first component is clearly significant while other 3 components also have a trend (P about 0.1). Any links to such graphs will also be helpful. Thanks in advance.

  • 1
    $\begingroup$ What do you mean by "show" a model? Plot predicted values against the predictor values? $\endgroup$ Jun 5 '15 at 11:59
  • $\begingroup$ Can we display wt, wt', wt'', wt''' clearly on one plot. $\endgroup$
    – rnso
    Jun 5 '15 at 12:56
  1. Choose a set of predictor values over which to calculate the expected response.

  2. Evaluate the natural spline basis using the same knot locations you used for the predictor in your model.

  3. Cross-multiply the spline basis by the coefficient estimates obtained from your model to get predicted values of the response for each predictor value.

  4. Plot the predictions against the predictor values.

(In the rms package the Predict function conveniently does all this for you.)

So this is just the same as what you'd do when using a polynomial basis for a predictor $x$, except that the basis functions are more complicated than $x, x^2, x^3, \ldots$. There's nothing to stop you looking at the effect of each basis function individually, but no reason to do so—one can't be held constant while the others vary, even in principle.

  • $\begingroup$ How would you explain in simple, layman terms what rcs() function does? I still have doubts regarding this. $\endgroup$
    – rnso
    Jun 5 '15 at 16:35

In simpler terms, what the rcs output shows are the betas for parts of the fit (each segment between the knots). If we plot the fit you have:

mtcars$fit = predict(modr)
ggplot(mtcars, ) +
  geom_point(aes(wt, mpg)) +
  geom_line(aes(wt, fit), color = "orange")

enter image description here

Since you have 4 betas, there are 3 knots. However, note that you need to add them up to get the one for a particular segment.

So, in your case:

  • Segment 1 has beta of -10.4203
  • Segment 2 has beta of -10.4203 + 19.3800 = 8.96
  • Segment 3 has beta of -10.4203 + 19.3800 + -68.7537 = -59.8
  • Segment 4 has beta of -10.4203 + 19.3800 + -68.7537 + 231.1378 = 171

I don't really know how to interpret these values, but essentially they give the idea that the first part (about 1.5 to 2.5) has negative slope, the next part (2.5 to 3.25) is less negative, then more negative (3.25 to 3.75), and then less negative (3.75 to 5.5). This fits the plot of the function fit, but does not really fit in terms of upwards vs. downwards slope.

I eyeballed the knot positions given above, but you can get them by:

> modr$Design$parms$wt
[1] 2.14 2.70 3.33 3.57 3.85

I don't really understand how these work in detail, I treat them as semi-understandable magic that approximates nonlinear functions without spending too many degrees of freedom.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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