[R-sig-ME] Model with interaction alone

[Ben Bolker]

Pierre Montpied <montpied at ...> writes:

> 
> Hi list,
> 
> I could not find any answer to my problem thus I post it on this list.
> 
> I have a design where  sites (Site) are nested in different regions 
> (Reg) and a pair of plots in each site
> is assigned two levels of a treatment (Trait).
> 
> I fit this model with lme (or lmer)  with  two crossed fixed factors 
> (Reg *Trait) and a random factor Site
> 
> All works fine when I fit the complete model "Response~Reg + Treat + 
> Reg:Treat"t or equivalently "Response~Reg*Trait"
> 
> But when I try the interaction alone by "Response~Reg:Trait" I get the 
> error message:
> "Error in MEEM(object, conLin, control$niterEM) :   Singularity in 
> backsolve at level 0, block 1".
> 
> When I replace "Reg:Trait" by "interaction(Reg,Trait)" or by 
> "factor(Reg:Trait)" all works fine
> whereas these should be equivalent to "Reg:Trait"
> 
> Any explanation ?
> 
> Below is a script simulating this


Thanks for the reproducible example.  I have to admit I don't fully
understand what's going on here, and in particular the difference
between "Reg:Trait" and "interaction(Reg,Trait)" and "factor(Reg:Trait)",
but the bottom line is that in the first case, R has decided that
it should include an intercept column in the model (fixed-effect
model matrix) as well as the 4 levels of the interaction, so the
fixed-effect model is overparameterized, which leads to the error.
In the second two cases R drops the intercept, so the model matrix
constructed by model.matrix() has only 4 rather than 5 columns
(see the end of this post for example).
The equivalent error from lme4::lmer is *slightly* more transparent
(and I'll try to improve it further).  If you specify Reg:Trait-1 or
Reg:Trait+0 (explicitly suppressing the intercept), everything works
again.

> 
library(nlme)
 
Nrep <- 10 # replicates = site number per region
Data <- expand.grid(Trait=c("T1","T2"),
                    Reg=rep(c("R1","R2"),each=Nrep),
                    Site=paste0("S",rep(1:(Nrep*2),each=2)))
Data <- transform(Data,RegTrait=factor(paste(Data$Reg,Data$Trait,sep=".")))
## with(Data,table(Reg,Trait,Site))
with(Data,table(Site,Trait,Reg))
 
Mu<-20    # Intercept
Alpha2<-1    # Fixed effect Trait
Beta2<-2    # Fixed effect Region
Gamma22<-2    # Interaction Trait:Region
Sigi<-1.5    # Random effect Site
Sige<-1    # Residual
 
## Simulation:
set.seed(101)
Data$Rep<- Mu +
         Alpha2*(Data$Trait=="T2") +
          Beta2*(Data$Reg=="R2") +
          Gamma22*(Data$Reg=="R2")*(Data$Trait=="T2") +
          rep(rnorm(Nrep*2,0,Sigi),each=2) +
          rnorm(Nrep*4,0,Sige)

library(nlme)

## the following two models should be _exactly_ identical
Modele.lme.1 <-lme(Rep~Trait*Reg,
     random=~1|Site,
     data=Data)# OK
 
Modele.lme.2 <-lme(Rep~Trait + Reg + Trait:Reg,
      random=~1|Site,
      data=Data)# OK
 
## Interaction alone:
 
Modele.lme.3 <- lme(Rep~Trait:Reg,
      random=~1|Site,
      data=Data)# Error:
## Error in MEEM(object, conLin, control$niterEM) :
##  Singularity in backsolve at level 0, block 1

Modele.lme.4 <- lme(Rep~interaction(Trait,Reg),
      random=~1|Site,
      data=Data)# OK
names(fixef(Modele.lme.4))
 
Modele.lme.5 <- lme(Rep~factor(Trait:Reg),
      random=~1|Site,
      data=Data)# OK
names(fixef(Modele.lme.5))

Modele.lme.6 <- lme(Rep~Trait:Reg-1,
      random=~1|Site,
      data=Data)# OK
names(fixef(Modele.lme.6))

library(lme4)
Modele.lmer.1 <- lmer(Rep~Trait*Reg+(1|Site),
                      data=Data)
Modele.lmer.2 <- lmer(Rep~Trait:Reg+(1|Site),
                      data=Data)  ## error: rank of X 
Modele.lmer.3 <- lmer(Rep~Trait:Reg-1+(1|Site),
                      data=Data)
fixef(Modele.lmer.3)

## check model matrices
X1 <- model.matrix(~Trait*Reg,Data)
X2 <- model.matrix(~Trait:Reg,Data)
X3 <- model.matrix(~Trait:Reg-1,Data)