[R-sig-ME] Split-plot Design

[Douglas Bates]

Thanks for your response John.  I just have one quick comment about
the Kenward-Roger degrees of freedom calculation.  It has been some
time since I looked at that paper but my impression at the time was
that the equations would not easily translate into the formulation
used in lmer.  The approach used in lmer is like that of Henderson's
mixed model equations (with many modifications).  That is, it is based
on a penalized least squares problem, not a generalized least squares
problem.  My recollection is that Kenward and Roger wrote their
equations in terms of the generalized least squares problem.

You are not the first person to suggest that incorporating the
Kenward-Roger calculation would enhance lmer.  The reason I haven't
done so is that I believe it is far from trivial to do so and I have
many, many other enhancements I would prefer to spend my time on.
However, this is open source software and any enterprising person who
wants to implement it is more than welcome to do so.  It may be
sufficiently involved to be a thesis topic - I don't know because I
haven't studied it carefully enough.

Any person considering doing that should read or reread Bill Venables'
"Exegeses on Linear Models" before embarking on it.  As he points out
in his discussion of modifications made in S-PLUS to emulate some
calculations in SAS, the "brute force" approach of taking a set of
equations and implementing them literally is rarely a good approach.
So I am not talking about a "pidgin R" implementation here where a
linear least squares calculation is written

XpX <- t(X) %*% X
XpXinv <- solve(XpX)
Xpy <- t(X) %*% y
betahat <- XpXinv %*% Xpy

An mer object includes slots L, RZX and RX that define the Cholesky
decomposition of the crossproduct matrix that is more-or-less like
that in the Henderson mixed model equations.  The L slot itself has a
Perm slot that gives the fill-reducing permutation P for the random
effects.  That may not be relevant - I'm not sure.  The original model
matrices are available as the X and Zt (transpose of Z) slots.  The
(transpose of the) derived model matrix for the orthogonal random
effects is available as A.  The terms attribute of the model matrix X
and the assign attribute of the model frame (in the slot named
"frame") should be used to associate terms with columns of the model
matrix.

It's possible that the calculation would be straightforward.  As i
said, I don't know.  My gut feeling is that it is not, which is why I
haven't embarked on it.