How to avoid singularity error in fitting piece-wise linear mixed effect regression model in R?

Question

I am trying to fit a linear mixed effect model to longitudinal haemoglobin (Hb) count data from a two-arm clinical trial. A subset of the data in long form is shown below:

          pid  site  arm pday Hb  pday2 I
1 MOB2004_001 Boane  SP    0  8.5     0 0
2 MOB2004_001 Boane  SP    3  8.0     9 0
3 MOB2004_001 Boane  SP    7  8.9    49 1
4 MOB2004_001 Boane  SP   14  8.9   196 1
5 MOB2004_002 Boane  SP    0 11.5     0 0
6 MOB2004_002 Boane  SP    3  9.1     9 0

The figure attached shows the mean response curves for the observations partitioned by study arm. The shape of the curve is such that I would like to model a quadratic relationship between haemoglobin and time, with different parameters for the curve before a knot at day 3. To this extent, I created an indicator variable 'I' which is zero when day < 3 and one otherwise. The objective is to assess differences between the study arms and thus a studyarm-day interaction is necessary too.

I have thus specified this model as:

lme(data=dat.mod2, fixed=Hb~ arm*pday + pday*I + pday2 + 
    pday2*I , random=~1|site/pid, na.action = na.omit)

The problem occurs when trying to fit this model, I get the error:

Error in MEEM(object, conLin, control$niterEM) : Singularity in backsolve at level 0, block 1

I see that this is a multicollinearity error, but I don't understand why it is occurring in this case.

I am aware that there are other ways to model this but want to use linear mixed effect models first. Would really appreciate any ideas on how to specify this model correctly.

I have tried different model specifications and all seem to work unless both pday*I and pday2*I are included. The model fits when I have either pday*I or pday2*I as fixed effects, but not with both. This is not great because I would like to estimate quadratic terms for both before and after the knot, and do not want to smooth over these.

I tried to fit with lmer in lme4 but I get the rank deficient error.

Answer 1

Based on the data you've shown, I suspect you may only have four unique time points (0, 3, 7, 14). Ignoring all the (I think) extraneous components, we can boil the fixed-effect part of your model down to this (I don't think discarding arm will matter for this question):

dd <- data.frame(day = c(0, 3, 7, 14))
dd <- transform(dd, day2 = day^2, I = as.numeric(day>3))
X <- model.matrix(~ I * (day + day2), dd)

  (Intercept) I day day2 I:day I:day2
1           1 0   0    0     0      0
2           1 0   3    9     0      0
3           1 1   7   49     7     49
4           1 1  14  196    14    196

X has six columns, you only have four unique day values. If you drop I*day or I*day2 from the specification (e.g. ~I*day or ~I*day2) that brings it down to four, which is (just) identifiable, although you will get a perfect/zero-residual fit ...

I may have gotten some of the details wrong, but you should be able to diagnose the problem yourself this way.

If you build the model matrix and can't immediately see what's wrong you can apply caret::findLinearCombos() to it:

caret::findLinearCombos(X)
$linearCombos
$linearCombos[[1]]
[1] 5 2 3 4

$linearCombos[[2]]
[1] 6 2 3 4


$remove
[1] 5 6