lme4: constrain sigma to 0

Took me a while, but I have managed to write up some specifics of the 
problem that I had in mind.  The gremlins seem to have changed things 
since the last time I played around with this stuff, and what I am now
getting differs from what I recollect.  However there is still a bit of
a problem.

I communicated with both Doug Bates and Ben Bolker on this issue, but 
cannot now for the life of me find any record of the emails that went 
back and forth.  And I keep ***everything***.  (Always the way; 
everything but what you want is available.)

Anyhow --- I have attached what I hope is a clear write-up in pdf format 
and a script to demonstrate what goes on using simulated data.

I hope these get through ...

cheers,

Rolf Turner

Any chance someone could write out the specifics of such a model? If I

can wrap my head around it and its not too hard, I could try to throw it
into blme.

On Feb 11, 2014, at 6:18 PM, Rolf Turner wrote:

Several millennia ago I put in a feature-request that the facility for
constraining sigma to 0 be added. So far, as far as I can tell, it
hasn't been.

Apparently it is (very) difficult to add such a constraint due to
the way that lmer approaches the maximization of the likelihood.
(This, rather than any intrinsic mathematical block to such a
constraint.)

I wanted to be able to fit a fairly simple repeated measures model
with the covariance over time being an "arbitrary" n-x-n positive
definite matrix (where "n" is the number of time points).  This
can't be done using lmer() unless one can constrain the "overall"
variance to be zero, otherwise the diagonal of the covariance
matrix is un-identifiable.

Since one cannot constrain sigma to be 0, one can't fit this model
with lmer().  Bummer.

-------------- next part --------------
A non-text attachment was scrubbed...
Name: specifics.pdf
Type: application/pdf
Size: 33767 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20140213/a166695c/attachment-0001.pdf>
-------------- next part --------------
#
# Script demoScript
#

library(MASS)
library(lme4)

# Generate simulated data.  No attempt is made to make these data
# realistic.

N <- 2000
MU <- c(0,-1,2)
Sigma <- 2*matrix(c(1.00,0.50,0.25,
                    0.50,1.00,0.50,
                    0.25,0.50,1.00),ncol=3)
set.seed(42)
Y <- mvrnorm(N,MU,Sigma)
y <- as.vector(t(Y))
dat.sim <- data.frame(y=y,time=factor(rep(1:3,N)),
                          student=factor(rep(1:N,each=3)))

# Fit with lmer():
fit <- lmer(y ~ time + (time|student),data=dat.sim,REML=TRUE)
V   <- VarCorr(fit)
CM  <- V[[1]]
attributes(CM) <- attributes(CM)["dim"]
sig.hat <- attr(V,"sc")

# Fixed effect estimates:
m <- matrix(c(1,-1,-1,0,1,0,0,0,1),3,3)
print(m%*%apply(Y,2,mean)) # Aggrees with "fit".

# Standard errors:
print(sqrt(diag(CM))) # Agrees with the "Std.dev" column under
                      # "Random effects" from "fit".
# The following agree with each other and with the standard
# errors of the fixed effects from "fit".
print(sqrt(diag(m%*%var(Y)%*%t(m)))/sqrt(N))
print(sqrt(diag(CM + sig.hat^2 * m%*%t(m)))/sqrt(N))