How to combine trend component and seasonality component from time series model?

I saw some solutions for people who wanted to jump straight into an R package to do the work, but I am trying to do some of the process piece-by-piece. The model that I refer to is

$$Y_t = m_t + s_t + X_t,$$

where $$Y_t$$ are the observations, $$m_t$$ is the trend component, $$s_t$$ is the seasonal component, and $$X_t$$ is the rough component.

I'm using R and the dataset is AirPassengers. What I can understand is fitting a Loess or polynomial trend separately using either loess() or lm(...poly()). Also, I understand creating a series of indicator variables in a matrix and using lm() to find the coefficients for the seasonality component. This would in essence give me values for $$s_t$$ of the seasonality and $$m_t$$ of the trend separately, but I am not sure that I can simply add these two values together to find the rough. I tried looking at the residuals for this setup and they don't appear like white noise.

In the case of polynomial degree 3 trend, $$m_t = \beta_0 + \beta_1t+\beta_2t^2+\beta_3t^3.$$ Also, with a seasonal component using 12 months, there is $$s_t = \beta_4I_{t,1} + \beta_5I_{t,2} + \cdots + \beta_{14}I_{t,11},$$ where $$I_{t,1} \begin{cases} 1 & \text{if time t is Jan} \\ -1 & \text{if time t is Dec}\\ 0 & \text{otherwise} \end{cases}, I_{t,2} \begin{cases} 1 & \text{if time t is Feb} \\ -1 & \text{if time t is Dec}\\ 0 & \text{otherwise} \end{cases} ,\cdots, I_{t,11} \begin{cases} 1 & \text{if time t is Nov} \\ -1 & \text{if time t is Dec}\\ 0 & \text{otherwise} \end{cases}.$$

So my question is how do I find these components at the same time? The reason is that I would like to ultimately test different transformations of Box-Cox using $$R^2$$, but the improper looking rough makes me hesitant togo on further.