I made a time series decomposition with tbats. There is weekly and yearly seasonality in the data (and maybe also monthly - not really important for the question)

x <- msts(data, start=c(2005,1,1), seasonal.period=c(7,30.4,365.25))
fit <- tbats(x, use.box.cox=FALSE)
plot(fit)

As far as i understand the tbats result is an "additive decomposition", right?

I can now access the different parts of the tbats decomposition:

level <- as.numeric(tbats.components(fit)[,'level'])
season1 <- as.numeric(tbats.components(fit)[,'season1'])
season2 <- as.numeric(tbats.components(fit)[,'season2'])
season3 <- as.numeric(tbats.components(fit)[,'season3'])

Since i suppressed 'Box-Cox transformation' earlier, level is roughly the same then trend according to a post from Rob J Hyndman in the comments section here: here

In order to get the remainder part of the decomposition is it a legit way to just subtract all the parts from the original data?

remainder <- data - level - season1 - season2 - season3

What do i get when i use this:

y <- resid(fit)

I am a bit confused right now about the right way to do it... Many thanks for all your input!

Update:

resid(fit)
fit$errors

Those two are both the same. I guess these values are related to the tbats method. I am doubting that i can take them as the "remainder" of the decomposition? Is there a way to extract the remainder out of the tbats method? In his paper 'Forecasting time series with complex seasonal patterns using exponential smoothing' Rob J Hyndman shows remainder graphs for the tbats method so that's why i think it is possible to get in R as well. Anyone any thoughts about that?