You can incorporate this in the Bayesian framework by specifying a prior distribution for the age variable. And for the posterior, you have:
$$p(\theta|DI)\propto p(\theta|I)p(D|\theta I)$$
Now you simply take $D\equiv (18+)$ for example. This is no more difficult "in-principle" compared to when you actually do know the ages. The difference is that your likelihood function must be a cumulative distribution function instead of a density. As an example, suppose age is the only regressor you have (denoted $x_i$), and you are fitting a OLS model. This is for my benefit - but the generalisation is just details, rather than conceptual. If you have observed the ages exactly the likelihood function is:
$$p(y_1\dots y_N|x_1\dots x_N\alpha\beta\sigma I)=(2\pi\sigma^2)^{-\frac{N}{2}}\exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^{N}(y_i-\alpha-\beta x_i)^2\right)$$
But now suppose that the $(N+1)$th observation, you only observe that $L<x_{N+1}<U$. Lets call this piece of information $Z$. Now we can use the brilliant trick of introducing a nuisance parameter and then integrating it out again (via the sum rule). The nuisance parameter we introduce is $x_{N+1}$ (the actual unobserved age), and we have:
$$p(y_1\dots y_N y_{N+1}|x_1\dots x_N Z\alpha\beta\sigma I)=\int_{L}^{U} p(y_1\dots y_N y_{N+1}x_{N+1}|x_1\dots x_N Z\alpha\beta\sigma I)dx_{N+1}$$
Now we can split the integrand by using the product rule $P(AB|C)=P(A|C)P(B|AC)$ and we get:
$$p(x_{N+1}|x_1\dots x_N Z\alpha\beta\sigma I)p(y_1\dots y_N y_{N+1}|x_1\dots x_N x_{N+1}Z\alpha\beta\sigma I)$$
Note that in the second density, the information $Z\equiv L<x_{N+1}<U$ is redundant because we are already conditioning on the true value $x_{N+1}$. So we can remove it. Note that this second term could be called the "clean" data. The first term is basically a statement of how likely the unobserved age is given $L<x_{N+1}<U$, in addition to the position of the "true line" $(\alpha,\beta)$, the noise level $\sigma$, and the values of all other ages $(x_1\dots x_N)$. And so you have an integrated likelihood (sometimes called quasi-likelihood):
$$p(Y|XZ\alpha\beta\sigma I)=(2\pi\sigma^2)^{-\frac{N+1}{2}}\exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^{N}(y_i-\alpha-\beta x_i)^2\right)$$
$$\times\int_{L}^{U}p(x_{N+1}|X\alpha\beta\sigma I)\exp\left(-\frac{(y_{N+1}-\alpha-\beta x_{N+1})^2}{2\sigma^2}\right)dx_{N+1}$$
Now for every "messy" data, you will have a similar integral. You can take the above integral as multi-dimensional (with appropriate matrix sum of squares in the exponential).
I have heard something like this called the "Missing information Principle". You basically create a "nice" dataset from your "messy" one (i.e. the data set you wish you had), and then average out the "nice" inferences. You give more weight to certain nice data sets according to what your "messy" information is.
