Best way to incorporate variable that is not linear among other linear variables in multiple regression? I am trying to predict a variable called epa, using 4 other variables. 3 of them are linear but this one, as you can see, is curved slightly downwards on each end (I plotted this using a Loess function). Based on my (limited) knowledge of multiple regression, it would not be useful to throw this into the equation as if it is linear like the others. Would there be a better way to incorporate it? I was thinking maybe a categorical variable that identifies whether the value is < 25, between 25 and 75, or > 75. I'm guessing this might be suboptimal though. I've tried Googling a bit but haven't had much luck, probably because I don't know what to call this non-linear variable form.

Also, I do realize this variable isn't going to be very predictive regardless. The question is more for learning purposes :)

Here is a way you can do this with splines.
In general, we can model a non-linear process with a linear model using splines.  In this example, I will generate some data with the true functional form
$$ y = 2x + z + \sin(z) $$
I will estimate this function using
$$ \hat{y} = \beta_0 + \beta_1x + f(z) $$
Where $f(z)$ is a linear combination of basis functions.  For this example, I will be using a natural spline, but there are other ways to estimate $f(z)$ with different basis functions.
Here is some code to do this in R.
library(tidyverse)
library(modelr)
library(splines)

# SImulate some data
N = 10000
x = rnorm(N)
z = rnorm(N)

# This is the true functional form
y = 2 * x + z + sin(z) + rnorm(N)

# Store data in a tibble
d = tibble(x, y, z)

# Basis functions require "knots".  Here, I am choosing them 
# relatively unintelligently.
k = quantile(z, seq(0.05, 0.95, length.out = 4))

# Estimate the model.  Here the spline is specified by ns
# ns <-> natural spline
model = lm(y ~ x + ns(z, knots = k), data = d)

# Predict on new data and plot
d %>%
  data_grid(x = c(-1,0, 1), z) %>%
  add_predictions(model) %>%
  ggplot(aes(z, pred, color = factor(x))) +
  geom_line()

Here is the result.  Note that $f(z)$ conditional on z is non-linear even though we estimated it through a linear model.  Additionally, as x changes, $f(z)$ moves up or down (simply because x is an additive effect).

