I was looking for some detailed proof for the simple linear model. By the factorisation theorem we have the following solution
Considering the classical model
$$y_i=\beta_0+\beta_1 x_i+\varepsilon_i,\ \mbox{ em que}\ \varepsilon_i\sim \mathrm{N}(0;\sigma^2),\ \hbox{Cov}(\varepsilon_i,\varepsilon_j)=0, i\neq j,\ i,j=1,\ldots, n$$
We start from $\mathbf{y}=(y_1,y_2,\ldots,y_n)$ a simple random sample and $x_1,x_2,\ldots, x_n$ not random. $\beta_0$ and $\beta_1$ are the classical parameters for the model (both maximum likehood or OLS). Consequently
$$y_i \sim \mathrm{N}(\beta_0+\beta_1 x_i;\sigma^2),\ i=1,\ldots, n.$$
Density function for $y_i$
$$f(y_i)= \frac{1}{\sigma \sqrt{2\pi}} \exp\left\{ -\frac{1}{2\sigma^2}\left[y_i-(\beta_0+\beta_1 x_i)\right]^2 \right\}\ -\infty <x_i<\infty $$
Likehood function
\begin{align*}
\hbox{L}(\mathbf{y}) &= \prod_{i=1}^{n} \frac{1}{\sigma \sqrt{2\pi}} \exp\left\{ -\frac{1}{2\sigma^2}\left[(y_i-\overline{y})-(\beta_0+\beta_1 x_i-\overline{y})\right]^2 \right\},\ -\infty <x_i<\infty \\
&= (2\pi\sigma^2)^{-n/2}\\
&= \exp\left\{ -\frac{1}{2\sigma^2}\left[\underbrace{\sum_{i=1}^{n}(y_i-\overline{y})^2-2\sum_{i=1}^{n}(y_i-\overline{y})(\beta_0+\beta_1 x_i-\overline{y})+\sum_{i=1}^{n}(\beta_0+\beta_1 x_i-\overline{y})^2}_{\delta}\right] \right\}
\end{align*}
Be $s_{yy}=\sum_{i=1}^{n}(y_i-\overline{y})^2=\sum_{i=1}^{n}y_i^2-n\overline{y}$ e $s_{xx}=\sum_{i=1}^{n}x_i^2-n\overline{x}$, $y$ random variable e $x$ constant.
\begin{align*}
\delta&=\sum_{i=1}^{n}(y_i-\overline{y})^2-2\sum_{i=1}^{n}(y_i-\overline{y})(\beta_0+\beta_1 x_i-\overline{y})+\sum_{i=1}^{n}(\beta_0+\beta_1 x_i-\overline{y})^2\\
\delta&=s_{yy}-2\sum_{i=1}^{n}(y_i\beta_0+\beta_1 x_iy_i-\overline{y}^2) +\sum_{i=1}^{n}\left[(\beta_0+\beta_1 x_i)^2-2(\beta_0+\beta_1 x_i)\overline{y}+\overline{y}^2\right]\\
\delta&=s_{yy}-2\alpha\sum_{i=1}^{n}y_i-2\beta\sum_{i=1}^{n}x_iy_i +2n\overline{y}^2+n\beta_0^2+2\beta_0\beta_1\sum_{i=1}^{n}x_i+\beta^2_1\sum_{i=1}^{n}x_i^2-2n\beta_0\overline{y}-2\beta_1\overline{y}\sum_{i=1}^{n}x_i +n\overline{y}^2\\
\delta&=\sum_{i=1}^{n}y_i^2-4n\beta_0\overline{y}+n\beta_0^2-2\beta_1\sum_{i=1}^{n}x_iy_i+2n\overline{y}^2+2\beta_0\beta_1\sum_{i=1}^{n}x_i+\beta^2_1\sum_{i=1}^{n}x_i^2-2n\beta_1\overline{x}\ \overline{y}\\
\delta&=\sum_{i=1}^{n}y_i^2+n\beta_0^2-4n\beta_0\overline{y}+2\beta_0\beta_1\sum_{i=1}^{n}x_i-2\beta_1\sum_{i=1}^{n}x_iy_i+2n\overline{y}^2+\beta^2_1\sum_{i=1}^{n}x_i^2-2n\beta_1\overline{x}\ \overline{y}\\
\delta&=\sum_{i=1}^{n}y_i^2-2n\beta_0 -2n\beta_0\left(\overline{y}-\beta_1\overline{x}\right)-2\beta_1\left(\sum_{i=1}^{n}x_iy_i-n\overline{x}\ \overline{y}\right)+2n\overline{y}^2+\beta^2_1\sum_{i=1}^{n}x_i^2\\
\delta&=\sum_{i=1}^{n}y_i^2+2n\overline{y}^2-2n\beta_0 -2n\beta_0\left(\overline{y}-\beta_1\overline{x}\right)+\beta^2_1\sum_{i=1}^{n}x_i^2-2\beta_1\left(\sum_{i=1}^{n}x_iy_i-n\overline{x}\ \overline{y}\right)\\
\delta &= \underbrace{\sum_{i=1}^{n}y_i^2+2n\overline{y}^2}_{g(\mathbf{y})}-\underbrace{2n\beta_0 -2n\beta_0\left(\overline{y}-\beta_1\overline{x}\right)}_{h(\mathbf{y},\beta_0)}+\underbrace{\beta^2_1 \sum_{i=1}^{n}x_i^2-2 s_{xx}\beta_1 \left(\frac{\sum_{i=1}^{n}x_iy_i-n\overline{x}\ \overline{y}}{s_{xx}}\right)}_{s(\mathbf{y},\beta_1)}
\end{align*}
Since $\hat{\beta_1}=\frac{\sum_{i=1}^{n}x_iy_i-n\overline{x}\ \overline{y}}{s_{xx}}$ and $\hat{\beta}_0=\overline{y}-\hat{\beta}\overline{x}$ are OLS and ML estimators
\begin{eqnarray*}
\hbox{L}(y_i) &=&(2\pi\sigma^2)^{-n/2} \exp\left[g(\mathbf{y})\right]\times \exp\left[h(\mathbf{y},\hat{\beta_0})+s(\mathbf{y},\hat{\beta_1})\right]
\end{eqnarray*}
Finally by the factorisation theorem we prove that for a simple linear regression model:
$$y_i=\beta_0+\beta_1x_i+\varepsilon_i,$$
the conditional distribution
$$Y|\hat{\beta}_0,\hat{\beta}_1$$
do not depends on $\beta_0$ and $\beta_1$, the original parameters for $Y$.
