Orange: Linear v. Additive v. Nonlinear Single and Multilevel Models

Orange Tree Models

###################################################
### linear v additive v nonlinear models
###################################################
# necessary packages
library(ggplot2)
library(lme4)
library(mgcv)

# rough idea of relationship
cor(Orange$age, Orange$circumference)

# plot model
ggplot(data = Orange, mapping = aes(x = circumference, y = age, color = Tree)) + geom_point() + geom_line()

## single level
# linear model
slm1 <- lm(formula = circumference ~ 1 + age, data = Orange)
summary(slm1)
plot(density(resid(object = slm1, type = "pearson")))

# linear model without intercept (when age is approx. 0 tree circumference is also approx. zero so constraining the intercepts to 0 might improve the model by adding a DF)
slm2 <- lm(formula = circumference ~ 0 + age, data = Orange)
summary(slm2)
plot(density(resid(object = slm2, type = "pearson")))

# additive model (single dimension, max 7 knots)
sgam <- gam(formula = circumference ~ 1 + s(age, bs = "tp", k = 7), data = Orange)
summary(sgam)
plot(density(resid(object = sgam, type = "pearson")))

# non-linear model
snls <- nls(formula = circumference ~ SSlogis(input = age, Asym = Asym, xmid = xmid, scal = scal), data = Orange)
summary(snls)
plot(density(resid(object = snls, type = "pearson")))



## multilevel model
# linear model
mlmer <- lmer(formula = circumference ~ 1 + age + (1 | Tree), REML = F, data = Orange)
summary(mlmer)
plot(density(resid(object = mlmer, type = "pearson")))

# linear model with covariance model
lme(fixed = circumference ~ 1 + age, random = ~ 1 | Tree, method = "ML", data = Orange)
mlme <- lme(fixed = circumference ~ 1 + age, random = ~ 1 | Tree, correlation = corLin(), method = "ML", data = Orange)
summary(mlme)
plot(density(resid(object = mlme, type = "pearson")))

# additive model (single dimension, max 7 knots)
mgam <- gam(formula = circumference ~ 1 + s(age, bs = "tp", k = 7) + s(Tree, bs = "re"), data = Orange)
summary(mgam)
plot(density(resid(object = mgam, type = "pearson")))

# non-linear model (starting values taken from single level model)
mnlmer <- nlmer(circumference ~ SSlogis(input = age, Asym, xmid, scal) ~ Asym | Tree, start = c(Asym = 193, xmid = 729, scal = 354), nAGQ = 0, data = Orange)
summary(mnlmer)
plot(density(resid(object = mnlmer, type = "pearson")))


data.frame(row.names = c("linear1", "linear2", "additive", "non-linear"), slevelBIC = sapply(X = list(slm1, slm2, sgam, snls), FUN = BIC), mlevelBIC = sapply(X = list(mlmer, mlme, mgam, mnlmer), FUN = BIC), slevelLL = sapply(X = list(slm1, slm2, sgam, snls), FUN = logLik), mlevelLL = sapply(X = list(mlmer, mlme, mgam, mnlmer), FUN = logLik))