How to create a B-spline basis without intercept and linear trend included? I want to fit the following model using splines:
\begin{align} Y(t) = \beta_0 + \beta_1t +  \sum_{j=2}^{d} \beta_jB(t)_j  \end{align}
where $B_j$ are the basis functions. However, when I run the following code:
fit <- lm(y ~ bs(t, df=5, intercept = FALSE), data = df)

The output only contains a coefficient for $\beta_0$ and the splines, not for the linear trend of $t$. When I add $t$ manually:
fit <- lm(y ~ t + bs(t, df=5, intercept = FALSE), data = df)

I do get an estimate for $\beta_1$, but the estimate for $\beta_d$, the coefficient for the last spline, is NA. Can I just remove the last column of the B-basis matrix and trust the results?
B <- bs(t, df=5, intercept = FALSE)
B <- B[,-ncol(B)]
fit <- lm(y ~ t + B, data = df)

The reason why I want to have the linear trend separately, is because in a next step I would like to extend the model by including a random intercept and a random slope:
\begin{align} Y(t) = (\beta_0 + b_0) + (\beta_1 + b_1)t +  \sum_{j=2}^{d} \beta_jB(t)_j  \end{align}
I know you can add a random intercept and random slope in R with the following code:
fit = lme(y ~ bs(t, df=5, intercept=FALSE),
                         random = ~ t | ID,
                         method = 'REML',
                         data = df)

And then you get indeed random intercepts, random slopes and a fixed coefficient for the intercept. But again not a fixed coefficient for the linear trend of $t$. Is this a mathematically correct model, since you have a random effect for $t$ but no "explicit" fixed effect? And is there a way to model the fixed effect for $t$ explicitly? For instance, by removing the last spline function from the design matrix of the B-splines and adding a linear trend manually?
If this is not possible with B-splines, is there a way to do it with natural cubic splines?

Your instinct is right--you can just drop one of the spline bases; and yes, this can be done with any kind of spline or indeed any basis whatsoever.  All you have to do is remove the $t$ component from all the other variables and then proceed.
The procedure can be coded efficiently by constructing the model matrix -- throw in everything you can think of--; using that as the response in a multivariate regression; and extracting its residuals:
library(splines) # Exports `bs`
k <- 5
X <- model.matrix(y ~ t + bs(t, df = k), df)
X. <- residuals(lm(X ~ df$t))

That's it: the last line is all the extra work that might be needed.  The rest is just organizing your results to make them easy to use.
Some of these columns might be essentially zeros.  You can detect them from the norms of the columns and remove them:
q <- colSums(X.^2)
X. <- X.[, zapsmall(q) > 0]

The adjusted model matrix will likely have a superfluous column--you may select one and remove it at this juncture if you wish.
Adjoin a column for $t$ and another one for your response vector $y$ and perform the regression:
colnames(X.) <- paste0("S.", seq_len(ncol(X.)) + 1)
Y <- cbind(subset(df, select = c("t", "y")), as.data.frame(X.))
obj <- lm(y ~ ., Y)

Let's examine a small example.  Here's a time series of 30 random variables with a slight linear trend:
n <- 30
set.seed(17)
df <- data.frame(t = seq_len(n), y = rnorm(n, seq_len(n)/n))

I wrote a function to plot the contribution of each variable to $y$ as a function of time.  (These are line plots of the $(t_i, b_j x_{ij})$ pairs where $i$ indexes the observations, $j$ indexes the variables, and $b_j$ is the estimated coefficient of $x_j.$)
Here's the original fit from lm(y ~ t + bs(t, df = k), data = df).

You can see there is a time term (upper left) and that the last spline term had been automatically dropped (bottom right).
Here's the fit based on the adjusted model matrix:

The spline terms look a little different because the contribution of $t$ has been removed from each of them: they have been "untilted."
If instead we had fit a model without $t$ originally (lm(y ~ bs(t, df = k), data = df)), here is what the fit would look like:

Each of these plots has some net nonzero slope: in this sense, $t$ is a component of all of them.  That's why you might have had some trouble originally in isolating the effect of $t$.
All three models give identical predictions.  If you do your significance tests correctly (by "chunking" all the spline components related to $t,$ including $t$ itself) they will also give identical p-values.  By plotting the models in this form you are giving yourself an additional way to understand their similarities and differences, perhaps enabling you to choose one you like and can interpret best.

I hope it is evident that (1) the nature and dimension of the splines doesn't affect this approach and (2) you can play this game not only with a time variable, but with any variable (or, indeed, set of variables) you wish to isolate and emphasize in your modeling.
