Accounting for pretreatment differences (linear mixed effects model) using R/nlme

Question

I have an experiment that is as follows:

• 4 years of data (1 year pre-treatment, 3 years during treatment)
• 20 plots total
• 5 plots of each type (Control, Treatment A, Treatment B, Treatment A+B)
• Each plots is divided into 4 quadrants in which measurements are recorded.

Initial stats revealed pretreatment differences between plot types.

The data I would like to use in a linear mixed-effects model include:

1. my response variable;
2. BA - before/after (as "1" or "2");
3. Treatment A (as "yes" or "no");
4. Treatment B (as "yes" or "no");
5. The interaction of treatments and Year;
6. Plot and quadrant as nested random effects.

I've created a linear mixed effects model using the nlme package in R that currently looks like this:

model<-lme(response variable ~ BA * treatmentA * treatmentB * year ,
           method= "REML", data=data, random=~1|plot/quad, na.action=na.omit)

I get the following error message:

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

I can run the model successfully without BA or without year but am unable to get results when I run it with both of them.

Any advice or suggestions on how to account for pretreatment differences and get my LME to work happily?

Answer 1

Hello Daniel and welcome to CV. What you are witnessing are most probably problems with collinearity between the columns of your design matrix. Drawing analogy with the solution of a weighted least squares regression you fail to compute $(X^TWX)^{-1}$. In your case I suspect that you are defining a complex interaction matrix between BA * treatmentA * treatmentB * year that ends up being collinear with the intercept. So my immediate response would be to try fitting your model as:

model<-lme(response variable ~ -1 + BA * treatmentA * treatmentB * year ,  ... )

If this works take notice that you have done: There is no mean response in your model anymore; understand the implication this has to your modelling assumptions.

Assuming that this doesn't work the slightly nastier version will be to fit the standard OLS model and the feed the design matrix of that model to your LME model. (Dont' panic, I explain and provide code) In R-terms, you'll be doing this:

OLSmodel <- lm(response variable ~ -1 + BA * treatmentA * treatmentB * year,...) 
ModelMatrixOLS  <- model.matrix(OLSmodel)
NonApplicable   <- which(is.na(OLSmodel$coefficients))
X_Matrix  <- ModelMatrixOLS[,-NonApplicable]
model<-lme(response variable ~ -1 + X_Matrix , ...)

So what have you done here is the following:

• Line 1 : You instructed R to fit a matrix that you know/suspect is rank deficient. R's work-horse lm() fitted it by taking care automatically of the collinearities

• Line 2 : You get you model matrix X

• Line 3 : You find which of the original coefficient of the model matrix you tried to fit, R deemed as collinear with some other and set them to NA (practically realized it couldn't estimate them)

• Line 4 : You get a new model matrix X where you have omitted the columns that R found problematic.

• Line 5 : You fit your LME model, without an intercept, using the model matrix that is not rank-deficient anymore.

Elegant? No. Works? Most probably. Under this scenario MAYBE you can run the original OLS model with intercept also but it wouldn't be my first choice. Take notice that if you want to provide you own model matrix you need to set a no-intercept model (by the -1 argument). If you do not do that, the lm, lmer, lme etc. will append a column of 1 in your model matrix (and essentially nullify all the work you did building up you design matrix).

In the general case, assuming that even this fails... Transform your variables. It doesn't seem that you can meaningfully do it here but I am just saying; that might also be an option.

I am not really well-versed with lme, I work exclusively with it's sibling function lmer but in any case I think this will probably work. What you are experiencing has to do with problems of computing the model you define, not with pretreatment differences specifically. Good luck!