4
$\begingroup$

I am interested in comparing three models

  1. a linear regression model: mod_linear <- lm (dv ~ iv)
  2. a polynomial regression model: mod_polynomial <- lm (dv ~ iv + I(iv^2))
  3. a one-breakpoint piecewise/segmented linear model: mod_segmented <- segmented(lm(dv~iv), seg = ~iv, psi=NA)

iv=independent variable, dv= dependent variable

Note for the segmented model: It estimates a breakpoint/changepoint/shift/knot in the slope, i.e., two lines and one inflection point in which the slope bends. This inflection point is itself a parameter (psi); it is not pre-specified, but is estimated via a robust fitting algorithm (built into the segmented package in R).

My question is: How do I compare the models to each other?

A. Is there a formal test that I can use? I know that nested models can be compared using the anova() function for a chi-square test in R, but I am not sure if it is appropriate to use anova() to compare the polynomial and piecewise models to each other.

B. Is it appropriate to use the R-squared values, residual standard error, AIC, and BIC, to compare the models?

Note: I am basing my reasoning for (A) on a paper where they compare piecewise and polynomial fits, and they state that "there is no formal test available that allows the direct comparison of a piecewise model to a power[polynomial] model" (pg 10, Ryan & Porth, 2007)

Ryan, S. E. & Porth, L. S. (2007). A tutorial on the piecewise regression approach applied to bedload transport data. U.S. Department of Agriculture & Forest Service Rocky Mountain Research Station, Fort Collins, CO, United States. https://www.arlis.org/docs/vol1/173845244.pdf

Below are fitted lines from each model polynomial=blue, segmented/piecewise=red, data points=gray

Fitted lines from each model, polynomial=blue, segmented/piecewise=red

Reproducible code and sample data

## load packages required
library(segmented)
library(lme4)

## Sample data
iv <-c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,232,233,234,235,236,237,238,239,240,241,242,243,244,245,246,247,248,249,250,251,252,253,254,255,256,257,258,259,260,261,262,263,264,265,266,267,268,269,270,271,272,273,274,275,276,277,278,279,280,281,282,283,284,285,286,287,288,289,290,291,292,293,294,295,296,297,298,299,300,301,302,303,304,305,306,307,308,309,310,311,312)
dv <-c(0.8532,0.8532,0.8532,0.8532,0.8532,0.8532,0.8532,0.8532,0.8916,0.9,0.9,0.9,0.9,0.9,0.9,0.9,0.9372,0.9,0.9,0.9492,0.9492,0.9492,0.9492,0.9492,0.9984,0.9984,0.984,0.9984,0.984,0.9984,0.9984,0.9984,1.05,1.05,1.05,1.05,1.0812,1.0812,1.0812,1.0812,1.05,1.05,1.1016,1.1016,1.1016,1.1016,1.1016,1.1316,1.1316,1.1832,1.1832,1.1832,1.1832,1.1832,1.1544,1.1544,1.1544,1.2072,1.2072,1.2072,1.2072,1.236,1.2888,1.2888,1.2888,1.2888,1.236,1.2612,1.2612,1.2612,1.2612,1.2612,1.3164,1.3164,1.3416,1.3416,1.3416,1.3416,1.3956,1.3416,1.3704,1.3704,1.3704,1.3704,1.3704,1.4256,1.4256,1.4256,1.4256,1.4496,1.4496,1.4496,1.4496,1.482,1.482,1.482,1.482,1.482,1.482,1.482,1.5384,1.5384,1.5384,1.5384,1.5384,1.5384,1.5384,1.5384,1.5948,1.5948,1.5948,1.5948,1.5948,1.5948,1.5948,1.5948,1.6512,1.6512,1.6512,1.6512,1.6512,1.6512,1.6512,1.6512,1.6512,1.7076,1.7076,1.7076,1.7076,1.7076,1.7076,1.7076,1.7076,1.7076,1.7652,1.728,1.728,1.8408,1.7844,1.7844,1.7844,1.7652,1.8216,1.8216,1.8216,1.8792,1.8792,1.8792,1.8792,1.8792,1.9212,1.9368,1.9368,1.9788,1.9944,1.9788,1.9944,1.9944,2.052,2.052,2.1108,2.1108,2.1108,2.1684,2.1684,2.1684,2.1684,2.226,2.226,2.226,2.2848,2.2848,2.2848,2.2848,2.3436,2.3436,2.3436,2.4012,2.4012,2.4012,2.46,2.46,2.46,2.46,2.5188,2.5188,2.5188,2.5776,2.5776,2.5776,2.6244,2.6364,2.6832,2.6832,2.6832,2.6952,2.742,2.742,2.802,2.802,2.802,2.8608,2.8608,2.9196,2.9196,2.9196,2.9784,2.9784,2.9784,3.0384,3.0384,3.0972,3.0972,3.0972,3.0972,3.1572,3.1572,3.1572,3.216,3.216,3.2748,3.2676,3.2748,3.2748,3.3348,3.3264,3.3864,3.3864,3.4452,3.4452,3.5052,3.5052,3.5052,3.5652,3.624,3.624,3.624,3.684,3.6912,3.6912,3.744,3.8028,3.8028,3.87,3.87,3.8628,3.9228,3.9228,3.9828,3.9888,4.0416,4.056,4.1016,4.1076,4.1616,4.1676,4.2204,4.2276,4.2864,4.2936,4.3464,4.3464,4.4388,4.4064,4.4664,4.4664,4.4724,4.5324,4.5672,4.5996,4.626,4.7352,4.7352,4.7184,4.7784,4.794,4.8456,4.8372,4.8972,4.9728,5.016,5.0232,5.0832,5.142,5.0916,5.1348,5.1948,5.202,5.2692,5.2548,5.3292,5.3136,5.3736,5.3736,5.4336,5.4336,5.4996,5.4336,5.5068,5.4924,5.5752,5.6856,5.694,5.6184,5.7384,5.7972,5.7972,5.7972,5.8572,5.8716,5.9316,5.9172)

## Fit models
mod_linear <- lm (dv ~ iv)

mod_polynomial <- lm (dv ~ iv + I(iv^2))

mod_segmented <- segmented(mod_linear, seg.Z = ~iv, psi=NA,
            it=20, control=seg.control(n.boot=50, size.boot=1000, K=1)) 

## Model comparison
# Chi-square test
anova(mod_linear, mod_polynomial, mod_segmented)

# Parsimony criteria: AIC and BIC 
AIC(mod_linear, mod_polynomial, mod_segmented)
BIC (mod_linear, mod_polynomial, mod_segmented)

## Plot fitted models
plot(dv ~ iv , col="grey")
lines(predict(mod_polynomial), col="blue", lwd=3)
plot(mod_segmented, add=TRUE, col="red", lwd=3)

Results from model comparison

> # Chi-square test
> anova(mod_linear, mod_polynomial, mod_segmented)

Analysis of Variance Table

Model 1: dv ~ iv
Model 2: dv ~ iv + I(iv^2)
Model 3: dv ~ iv + U1.iv + psi1.iv

  Res.Df    RSS Df Sum of Sq        F    Pr(>F)    
1    310 57.377                                    
2    309  2.307  1    55.070 8338.356 < 2.2e-16 ***
3    308  2.034  1     0.273   41.334 4.883e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> 
> # Parsimony criteria: AIC and BIC 
> AIC(mod_linear, mod_polynomial, mod_segmented)

               df       AIC
mod_linear      3  363.0917
mod_polynomial  4 -637.5626
mod_segmented   5 -674.8518

> BIC (mod_linear, mod_polynomial, mod_segmented)

               df       BIC
mod_linear      3  374.3207
mod_polynomial  4 -622.5906
mod_segmented   5 -656.1368
```
$\endgroup$
2
$\begingroup$
  • A polynomial regression model can not be represented as a segmented model, and

  • a segmented regression model cannot be represented as a polynomial model,

so, this models are not nested (but both are, separately, nested within your first simple regression model). There do exist tests for non-nested models, see this posts: Generalized log likelihood ratio test for non-nested models, Non-nested model selection and for problems with the usual tests Why can't likelihood ratio tests be used for non-nested models?.

On this site it is often recomended to use information criteria (AIC, BIC), see Comparing non nested models with AIC

$\endgroup$
1
  • $\begingroup$ Thank you for the helpful answer and relevant posts! $\endgroup$
    – Annalise
    Apr 23 '21 at 13:52

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.