If I use GAM, what is the formula for the estimated smoother $s$, i.e. if we have the model $$y \sim s(x)$$ using mgcv?
I can plot the term of $s$, but that isn't really all that informative. I want the exact formula, so that I can relay it in my report?
There isn't an easily conveyable formula or equation for this model. An alternative way of writing it is in terms of the basis expansion of $x$. If we gather each of the basis functions in the expansion of $x$ (with suitable identifiability constraints applied), plus a column vector of 1s for the model intercept, your model is just
$$ \boldsymbol{y} = \boldsymbol{X\beta} + \boldsymbol{\varepsilon}$$
which is just a multiple linear regression (assuming the family was Gaussian).
As such, you can write out an equation for the model taking the estimated coefficients for each of the terms in $\boldsymbol{X}$ just as you would for a linear model.
The real question is whether this is of any use? By default, what mgcv is doing to the thin plate spline means that the basis functions may well not be all that informative.
If you want something to use to predict from the model in software outside of R, then you can use the predict(model, newdata, type = "lpmatrix") option.
Beyond that there isn't much to offer you; if you want a formula in terms of the basis functions you can write it as:
$$\hat{\boldsymbol{y}} = \beta_0 + \beta_1 \boldsymbol{x}_{b_1} + \beta_2 \boldsymbol{x}_{b_2} + \cdots + \beta_9 \boldsymbol{x}_{b_9}$$
assuming 9 basis functions for the default k = 10 and where $\boldsymbol{x}_{b_j}$ is the $j$th basis function of $x$, replacing the $\beta_j$ with the estimated values from coef(model).
