Exploring REML edge cases

twitter_badREML.R

library(lme4)
library(nlme)
library(glmmTMB)
library(Matrix)

## what do various R packages do when confronted with an
## undefined REML problem (fixed effects *and* random effects
## for every subject?

## Twitter thread: https://twitter.com/ten_photos/status/1657399290166222850


N <- 100
K <- 10
set.seed(101)
Y  <-  rnorm(N*K) # iid sim is OK to make the point here
Subj  <-  rep(1:N, each = K)
X  <-  factor(Subj)
dd <- data.frame(Y, X, Subj)
summary(lmer(Y~  (1|Subj), dd)) # This should be fine
m2 <- lmer(Y~X+(1|Subj), dd)
## Warnings
## unable to evaluate scaled gradient
##  Hessian is numerically singular: parameters are not uniquely determined

VarCorr(m2)
m2@optinfo$derivs$Hessian  ## 0!

m3 <- lme(Y~X, random = ~ 1|Subj, method = "REML")
intervals(m3)
## Error in intervals.lme(m3) : 
##   cannot get confidence intervals on var-cov components: Non-positive definite a

## eq. 42 of vignette
Ldet <- det(getME(m2, "L"))  # 2e16
RX <- getME(m2, "RX")
RXdet <- det(RX) # 4e33
##image(Matrix(getME(m2,"RX")))

2*(log(Ldet) + log(RXdet))  ## log-det component of REML criterion

m4 <- glmmTMB(Y~X + (1|Subj), REML = TRUE, dd)
summary(m4)
VarCorr(m4)
confint(m4)  ## stand d

pp <- with(m4$obj$env, last.par.best[-random])
## hessian calc
H <- numDeriv::jacobian(m4$obj$gr, x = pp)
det(H) ## -7e-6