# Is it possible to specify a lmer model without any fixed effects?

In R, how do I specify lmer model without global fixed effect? For example, if I say something like

lmer(y ~ (1 | group) + (0 + x | group), data = my_df)


the fitted model will be

$y_{ij} = a + \alpha_{i} + \beta_{i} x_{ij}$

How do I fit model

$y_{ij} = \alpha_{i} + \beta_{i} x_{ij}$?

-
My kneejerk answer was that lmer(y~0+(1|group)+(0+x|group)) would work, but this yields an error. –  Mike Lawrence Nov 30 '11 at 7:47
I do not think it is possible to specify a model without a fixed effect with lmer because the lme4 package is dedicated to mixed models only (with at least one fixed effect and one random effect). I do not remember seeing random effects models in the documentation, in any case. It could be useful to ask for it in the R-sig-mixed-models list (stat.ethz.ch/mailman/listinfo/r-sig-mixed-models) –  maxTC Nov 30 '11 at 13:38
Just out of curiosity, why would you want to do so? Maybe there's another way to approach your underlying goal. –  jbowman Nov 30 '11 at 21:37
center y first :) –  Macro Jul 14 '12 at 4:44

As @Mike Lawrence mentioned the obvious thing to do when defining a model without fixed effects is something in the form of:

lmer(y ~ -1 + (1|GroupIndicator))


which is actually quite straightforward; one defines no intercept or an X matrix. The basic reason which this doesn't work out is that as @maxTC pointed out "lme4 package is dedicated to mixed models only".

In particular what lmer() fitting does is calculate the profiled deviance by solving the penalized least square regression between the $\hat{y}$ and ${y}$ as well as the spherical random effects $u$ and $0$ (Eq. (11), Ref.(2)). Computationally this optimization procedure computes the Cholesky decomposition of the corresponding system exploiting the system's block structure (Eq. (5), Ref.(1)). Setting no global fixed effects practically distorts that block structure in a way that the code of lmer() can't cope. Among other things the conditional expected value of $u$ is based on $\hat{\beta}$'s but solving for $\hat{\beta}$ asks the solution of a matrix system that never existed (the matrix $R_{XX}$ in Ref.(1), or $L_X$ in Ref.(2)). So you get an error like:

Error in mer_finalize(ans) :
Cholmod error 'invalid xtype' at file:../Cholesky/cholmod_solve.c, line 970


cause after all there was nothing to solve for in the first place.

Assuming you don't want to re-write lmer() profiled deviance cost function the easiest solution is based on the CS-101 axiom: garbage in, garbage out.

 N = length(y); Garbage <- rnorm(N);
lmer(y ~ -1 + Garbage + (1|GroupIndicator));


So what we do is define a variable $Garbage$ that is just noise; as before lmer() is instructed to use no fixed intercept but only the X matrix defined us (in this case the single column matrix Garbage). This extra Gaussian noise variable will be in expectation uncorrelated to our sample measurement errors as well as with your random effects variance. Needless to say the more structure your model has the smaller the probability of getting unwanted but statistically significant random correlations.

So lmer() has a placebo $X$ variable (matrix) to play with, you in expectation will get the associated $\beta$ to be zero and you didn't have to normalize your data in any way (centring them, whitening them etc.). Probably trying a couple random initialization of the placebo $X$ matrix won't hurt either. A final note for the "Garbage": using Gaussian noise wasn't "accidental"; it has the largest entropy amongst all random variables of equal variance so the least chance of providing an information gain.

Clearly, this is more a computational trick than a solution, but it allows the user to effectively specify an lmer model without global fixed effect. Apologies for hoping around the two references. In general I think Ref.(1) is the best bet for anyone for realize what lmer() is doing, but Ref.(2) is closer to the spirit of the actual code.

Here's a bit of code show-casing the idea above:

library(lme4)
N= 500;                 #Number of Samples
nlevA = 25;             #Number of levels in the random effect
set.seed(0)             #Set the seed
e = rnorm(N); e = 1*(e - mean(e) )/sd(e); #Some errors

GroupIndicator =  sample(nlevA, N, replace=T)   #Random Nvel Classes

Q = lmer( rnorm(N) ~ (1| GroupIndicator ));      #Dummy regression to get the matrix Zt easily
Z = t(Q@Zt);            #Z matrix

RA <-  rnorm(nlevA )                        #Random Normal Matrix
gammas =c(3*RA/sd(RA))                      #Colour this a bit

y = as.vector(  Z %*% gammas +  e )         #Our measurements are the sum of measurement error (e) and Group specific variance

lmer_native <- lmer(y ~ -1 +(1| GroupIndicator)) #No luck here.
Garbage <- rnorm(N)                              #Prepare the garbage

lmer_fooled <- lmer(y ~ -1 + Garbage+(1| GroupIndicator)) #OK...
summary(lmer_fooled)                             #Hey, it sort of works!


References:

1. Linear mixed models and penalized least squares by D.M. Bates and S. DebRoy, Journal of Multivariate Analysis, Volume 91 Issue 1, October 2004 (Link to preprint)
2. Computational methods for mixed models by Douglas Bates, June 2012. (Link to source )
-
Is there any reason to prefer the approach you've put forward over the one mentioned by Marco in the comments? –  rpierce Mar 3 at 22:51
Yes; you don't alter the data so no back transformations are needed. As a result all the standard diagnostics and goodies by lmer() eg. fitted variables, residuals, random effects levels etc. etc. are directly interpretable as they correspond to your "true" dataset and not an "altered one". –  user11852 Mar 3 at 23:00