# comparing models in R

Consider the following example:

x <- c(0.25,0.5,0.75,1,1.25,1.5,1.75,2,2.25,2.5,2.75,3,3.25)
y <- c(0.516,0.125,-0.0781,0,0.453,1.38,2.86,5,7.89,11.6,16.3,22,28.8)
d <- data.frame(x = x,
y = y)

if(!require("segmented")) {
install.packages("segmented")
require("segmented")
}
if(!require("mgcv")) {
install.packages("mgcv")
require("mgcv")
}

par(mfrow = c(2,2))

# first model
g1 <- lm(y ~ x,data = d)
g2 <- segmented(g1, seg.Z = ~ x,
psi = list(x = c(1.5)))
pdat <- data.frame(x = d$x, y = broken.line(g2, link = FALSE)[,1]) pdat <- pdat[with(pdat, order(x)), ] plot(y ~ x, data = d, pch = 21, bg = "white") lines(y ~ x, data = pdat, type = "l", col = "red") # second model g3 <- glm(y ~ x,data = d) g4 <- segmented(g3, seg.Z = ~ x, psi = list(x = c(1.5))) pdat <- data.frame(x = d$x,
y = broken.line(g4, link = FALSE)[,1])
pdat <- pdat[with(pdat, order(x)), ]
plot(y ~ x, data = d, pch = 21, bg = "white")
lines(y ~ x, data = pdat, type = "l", col = "red")

# third model
g5 <- lm(y ~ poly(x, 2), data = d)
pdat <- with(d, data.frame(x = exp(seq(min(x),
max(x), length = 100))))
tmp2 <- predict(g5, newdata = pdat, se.fit = TRUE)
pdat <- transform(pdat, pred = tmp2$fit, se = tmp2$se.fit)

plot(y ~ x, data = d)
lines(pred ~ x, data = pdat4, type = "l", col = "red")


Here, I am fitting the same data with three models.

Is there a correct way of comparing different model types?

Usually I use AIC

> AIC(g2,g4,g5)
df      AIC
g2  5 45.06724
g4  5 43.06962
g5  4 32.49106


But I wasn't sure if this was accurate seeing as the models are different i.e. quadratic vs lm vs glm.

g1 (fit with OLS) and g3 (fit as a GLM with Gaussian family) are identical models and indeed:

all.equal(AIC(g1), AIC(g3))
#[1] TRUE
logLik(g1)
#'log Lik.' -36.83844 (df=3)
logLik(g3)
#'log Lik.' -36.83844 (df=3)


Thus, we would expect the AIC values of the corresponding segmented models also to be identical. That they aren't seems like an oversight:

set.seed(1)
g2 <- segmented(g1, seg.Z = ~ x,
psi = list(x = c(1.5)),
control=seg.control(toll=1e-100, it.max=100))
set.seed(1)
g4 <- segmented(g3, seg.Z = ~ x,
psi = list(x = c(1.5)),
control=seg.control(toll=1e-100, it.max=100))

all.equal(coef(g2), coef(g4))
#[1] TRUE


So, the coefficients are equal.

AIC(g2)
#[1] 45.06724
AIC(g4)
#[1] 43.06724
logLik(g2)
#'log Lik.' -17.53362 (df=5)
logLik(g4)
#'log Lik.' -16.53362 (df=5)


But the log-likelihood values differ exactly by 1. In help("logLik") we read:

For a "glm" fit the family does not have to specify how to calculate the log-likelihood, so this is based on using the family's aic() function to compute the AIC. For the gaussian, Gamma and inverse.gaussian families it assumed that the dispersion of the GLM is estimated and has been counted as a parameter in the AIC value, and for all other families it is assumed that the dispersion is known.

segmented doesn't seem to consider that adequately when calculating g4\$aic.

In summary, you should be careful when using AIC with segmented models.

Edit:

I've contacted Vito Muggeo, the maintainer of the segmented package. He answered:

Currently segmented.glm is not including the breakpoint parameter in the df. Therefore there are differences equal to '2*n.of.breaks' in the AIC values (and equal to "1*n.breaks" in the logLik values) [...] I have corrected the code. The next version of segmented will include this correction.

How to compare models depend by the aim of your research. If you are more interested in predictive performance, a statistic like the Root Mean Squared Error of prediction RMSE (check the caret package for a quick implemenation) would probably be ok. If you are more interested in probabilistic (likelihood based) fit then AIC could be a reasonable choice