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?
 A: When we use a spline in the context of an additive model as soon as the spline basis is created, fitting reverts to standard GLM modelling basis coefficients for each separate basis function.
The issue is that by construction, the B-spline is designed to be a representation of $x$ itself. Practically this means that given B-spline basis expansion $B$ if we append $x$ itself to that basis expansion then this is an over-determined system - that matrix is not full rank, ergo non-invertible, ergo the NA.
set.seed(123)
N = 1000
x = sort(runif(N, min = -10, max = 10)) 
library(splines)
Bmatrix <- bs(x = x,  df=5,  intercept = FALSE)
colnames(Bmatrix) <- paste0("BF", 1:5)
matplot( x, Bmatrix, type = "l", lwd=2); 
grid()
legend("topleft", legend = colnames(Bmatrix), lwd=2, col=1:5)


BmatrixPlus = cbind(Bmatrix, x)
# Do we have machine 0 eigenvalues?
any(eigen(cov(BmatrixPlus))$values <= .Machine$double.eps) # Yes... :(

The solution you propose is "good". In a sense, removing any of the basis functions and then adding $x$ ensure that we go around the rank deficiency. But that creates the question, "which basis function to remove?". Realistically, either we start looking at which one we "like the least"/our expert knowledge suggests it is a very unlikely pattern, or we pick one at random. The first option is not "that bad" but we maybe we can do better...
Enter Legendre polynomials; Legendre polynomials (LPs) form a basis expansion on $[-1,1]$ and they can easily be redefined on a shifted space, LPs are considered "classical orthogonal polynomials" and most importantly for us they also define a representation of $x$ that uses $x$ directly (or a directly scaled combination of it) as a building block. That means we can immediately use LPs instead of an intercept and slope in our design matrix:
library(orthopolynom)
LP5coef = legendre.polynomials(5, normalized = FALSE) # Read these
LPmatrix <- as.matrix(as.data.frame(polynomial.values(polynomials=LP5coef, x=x/max(x))))
colnames(LPmatrix) <- paste0("LP", 0:5)
matplot(x, LPmatrix, type = "l", lwd=2); 
grid()
legend("topleft", legend = colnames(LPmatrix), lwd=2, col=1:6)


We can rename our first two LPs to "intercept" and "slope" if that makes it easier too. Please note that it makes sense to look at the Legendre polynomial coefficient to ensure they are at the scale you expect them to be.
Finally, let us note that there is also the "unloved option" of using a direct polynomial basic expansion of $x$ too, i.e. something in like poly(x, 5). That is not as catastrophic as people make it to be but it can be inefficient to capture some of the finer variances our $x$ might entail compared to a $B$-spline basis or a Legendre polynomial basis.
A: 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.
