Estimating a conditional logistic model with hierarchical splines in mgcv leads leads to error "indefinite penalized likelihood in gam.fit5" One can estimate a conditional logistic model in the R package mgcv by using the cox.ph family, which I have done successfully. I also estimated a logistic regression model with hierarchical spline coefficients.
However, when I try to estimate a conditional logistic model with hierarchical spline coefficients as follows:
fit = gam(cbind(D,ID) ~ te(X, bs = "tp", m = 2) + t2(X, edu, bs = c("tp","re"), m = 2, full = TRUE), 
          data = my_data,
          family = cox.ph,
          weights = is.case)

i get the error

Error in gam.fit5(x, y, sp, Sl = Sl, weights = weights, offset = offset,  :
indefinite penalized likelihood in gam.fit5

I'd be grateful for any hint pointing to a potential solution or mistake I am making.
 A: An issue with models of this form, where we essentially have multiple smooths of a covariate – in this case X – is that our models might be unidentifiable or over parameterised given the data.
In this case model the OP is attempting to estimate some average effect (what we called a "global" effect in Pedersen et al (2019)) plus subject specific smooths. Decomposing the data into those two components is not always possible given the system and the sample of data being analysed.
The problem may be due to there being little variation among subjects, and only the "global" effect is important, although this tends to be less of a problem as the specific parameterisation from the subject-specific smooths includes includes an (random) intercept for each subject. This typically results in very similar subject-specific smooths that are offset from one another by the subject-specific intercepts. These terms are fully penalized too – their null space, the constants (subject-specific intercepts) and linear components can be shrunk by the penalty on the smooth, unlike most other smooths in {mgcv} – so if you don't need subject-specific effects, the penalty can shrink them effectively out of the model entirely.
The other situation is that the "global" smooth is not identifiable or estimable given the data. In that case switching to the simpler model with just the subject specific smooths is a reasonable next modelling step. One might also fit a model with only the global smooth and compare it with the model containing just subject-specific smooths and see which explains the data better. To isolate the comparison to be on whether one needs subject-specific smooths or not, the global model could include a random intercept term for the subjects such that we compare
y ~ s(x) + s(subject, bs = "re")

with
y ~ s(x, subject, bs = "fs")

because the random factor smooth (which is equivalent in the univariate case to the `t2()`` term used by the OP) includes a subject-specific intercept also.
Another possibility could be that the fitting code is having a hard time getting good starting values. I would look to see if the error happens in the early phase of the optimization or whether it happens later, by turning on some tracing elements by adding
gam.control(trace = TRUE)

to the call to gam().
It can help sometimes to focus the subject-specific smooths on differences from the global. One way to do that is to put a penalty on the first derivative for the subject-specific terms, by adding m = 1 to the smooth definition: s(x, subject, bs = "fs", m = 1). Because the model is now penalizing the global and subject-specific smooths in different ways, it can help with fitting. That said, it has less of an effect with these random factor smooth models as those smooths are fully penalised.
