In a linear regression model $y= \beta_0+\beta_1X_1+\epsilon$, why is the model linear through the $\beta$'s and not $X_1$? In a linear regression model
$$
y_i= \beta_0+\beta_1X_{i1}+\epsilon_i
$$
I have seen it often stated that the model is linear through the $\beta$ 's  and not $X_1$. Specifically, if we have an example where $y_i$ is the sale price of home $i$ and $X_{i1}$ is square footage of home $i$, then I have seen that even if $\beta_{i2} = (\text{number of sinks})^2$, then the model is still linear in $\beta$. Why is that?
 A: The reason we focus on linearity with respect to the model parameters is that this is what determines the complexity of the mathematics for fitting the model and computing all the outputs.  In a "linear model" (i.e., linear with respect to the model parameters) we get an explicit form for the OLS estimator for the coefficients, and we likewise get nice explicit forms for all the various goodness-of-fit statistics in the model.  However, as soon as we have one or more parameters that enter the model in a nonlinear way (that is not "linearisable") we now have a "nonlinear model", and the OLS estimator and all subsequent quantities of interest no longer have simple closed forms; they have to be computed using iterative methods, and standard errors for our estimators are then approximated using the delta method.
Mathematically speaking ---in terms of fitting the model and computing all the outputs--- we really don't care about linearity with respect to the explanatory variables because we can easily compute the design matrix for these variables, with any nonlinear transformations we want on the variables.  However, if we have (non-reducible) nonlinearity in any of the parameters then the mathematics of the model becomes much more complicated.  Consequently, the distinction of interest for mathematicians/statisticians who build and analyse these models is between the nice simple "linear model" and the nasty complex "nonlinear model".

Example: Consider a periodic regression model of the form:
$$Y_t = \mu + \beta \sin(2 \pi \delta (\alpha + t)) + \varepsilon_t
\quad \quad \quad \varepsilon_t \sim \text{N}(0,\sigma^2).$$
In this model we have the following regression parameters:
$$\begin{align}
\text{Mean parameter} & & & & & \mu \\[6pt]
\text{Amplitude parameter} & & & & & \beta \\[6pt]
\text{Frequency parameter} & & & & & \delta \\[6pt]
\text{Phase-angle parameter} & & & & & \alpha \\[6pt]
\end{align}$$
Now, it turns out that we can write this model in a form that is linear with respect to all parameters except for the frequency parameter.  Specifically, using the rules for decomposing an out-of-phase sinusoidal wave we can rewrite the periodic regression model as:
$$Y_t = \mu + \beta_1 \sin(2 \pi \delta t) + \beta_2 \cos(2 \pi \delta t) + \varepsilon_t
\quad \quad \quad \varepsilon_t \sim \text{N}(0,\sigma^2),$$
using the transformed parameters $\beta_1 = \beta/\cos(2 \pi \delta \alpha)$ and $\beta_2 = \beta/\sin(2 \pi \delta \alpha)$.  This means that if the frequency is known (so that we are only estimating the other three parameters) then this is a "linear model" and we can fit it using the mathematics of multiple linear regression.  However, if the frequency parameter is unknown then this is a "nonlinear model" and we need to fit it using the mathematics of nonlinear regression.  This distinction is important because it makes a great deal of difference to the complexity of the model and the form of the resulting estimators, goodness-of-fit statistics, and standard errors of the various estimators.
