One option since you are considering nonlinear regression is a generalized additive model (GAM). I have provided an example below. First I loaded the mgcv package for running GAMs and MASS for the mcycle dataset, which is referenced in my link at the end of this answer.
#### Libraries ####
library(mgcv)
library(MASS)
I then fit a model. Its fairly similar to a normal regression fit using lm with some minor tweaks. First, I added a simple spline term here (basically a function which creates a smooth function for nonlinear modeling), then set the number of knots to 10, basically stating that the maximum amount of curving points should be around this number. Thereafter I just add the dataset and set the method to either REML or ML (these are the two suggested methods to use).
#### Fit Model ####
fit <- gam(
accel ~ s(times,
k=10),
data = mcycle,
method = "REML"
)
From there you can just summarize the model like an lm object.
#### Summarize and Check ####
summary(fit)
Which gives you these terms:
Family: gaussian
Link function: identity
Formula:
accel ~ s(times, k = 10)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -25.546 1.951 -13.09 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(times) 8.625 8.958 53.4 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.783 Deviance explained = 79.7%
-REML = 616.14 Scale est. = 506.35 n = 133
You'll see the summary looks different. Here, if you chose linear terms, they would be listed with the intercept under the parametric section. We only have one slope term, the smooth fit function for times, so we skip to that in the Approximate Smoothness section. Here we see the term is significant and the edf is around 8.625, indicating a fair amount of curvature. We can visualize it by simply plotting it directly. Here I've added the residual argument to show the data points and changed the point shapes.
#### Plot Model ####
plot(fit,
residuals = T,
pch=21)
You can see that the data clearly has a nonlinear function, so fitting a linear model would be fairly bad. We can also check to make sure our model doesn't need to be adjusted (this is simply a heuristic check and if p values are flagged you can consider adjusting your splines to fit better). It automatically spits out a plot of the fitted/response values as well.
gam.check(fit)
The summary and plot are seen below:
Method: REML Optimizer: outer newton
full convergence after 7 iterations.
Gradient range [-1.436342e-06,1.179834e-06]
(score 616.142 & scale 506.3529).
Hessian positive definite, eigenvalue range [3.329028,65.73378].
Model rank = 10 / 10
Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.
k' edf k-index p-value
s(times) 9.00 8.62 1.15 0.95
You should get 4 plots, and the qq plot this time looks normal.
There appear to be no issues with our check, so we are free to move on.
Useful Resources
For learning all about GAMs, the DataCamp course is great, but this is also an excellent free resource:
https://noamross.github.io/gams-in-r-course/
