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?

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.
