kollerma / robustlmm

This is an R-package for fitting linear mixed effects models in a robust manner. The method is based on the robustification of the scoring equations and an application of the Design Adaptive Scale approach.
28 stars 9 forks source link

missing value where TRUE/FALSE needed in updateSigma() #27

Closed jay-sf closed 1 year ago

jay-sf commented 1 year ago

Hi kollerma,

first of all thanks for the robustlmm package!

I face an issue in an rlmer call while using weights.

In updateSigma function, line 376 scale <- fun2(scale0, object@resp$wtres) can become NaN which leads to a NA value in converged which in turn leads, while trying to evaluate !converged && (it <- it + 1) < max.iter to an error. Couldn't find out what goes wrong.

Self-contained example:

n <- 200
set.seed(42)
data <- matrix(sample(0:1, n*8, replace=TRUE, prob=c(.9, .1)), n, 8)
data <- data.frame(id=seq_len(n), y=cbind(1, data, runif(n)) %*% rnorm(ncol(data) + 2), 
                   data, runif(n)) |>
  setNames(c("id", "y", "cond", "treat", "d1", "d2", "d3", "d4", "d5", "x", "w"))
data <- merge(data, expand.grid(id=1:n, trial=1:10))
data$y <- data$y*rnorm(nrow(data))

> rlmer(y ~ cond + treat + d1 + d2 + d3 + d4 + d5 + (1 | id) + (1 | trial) + x +
+         cond:treat + cond:d1 + cond:d2 + cond:d3 + cond:d4 + cond:d5, 
+       data, weights=w)
boundary (singular) fit: see help('isSingular')
Error in while (!conv && (it <- it + 1) < max.it && tau[i] <= 1) { : 
  missing value where TRUE/FALSE needed
In addition: Warning message:
In sqrt(tau2) : NaNs produced

The boundary message and sqrt(tau2) warning do not appear with my real data.

It would be great if you could look at this.

Cheers!

kollerma commented 1 year ago

Thanks for reporting this. I suspect this model is too complex for the data and results in a degenerated fit. Use the verbose argument to find out what happens during the fitting process.

Your example crashes R, which shouldn't happen. I don't have time to investigate right now, but I expect to be able to investigate it in a few weeks at least.

kollerma commented 1 year ago

I finally found some time to investigate this. Some of the weights are very small. The weights are so small that the formula used to approximate the variance of the individual residual's variance produces negative values. Negative values don't make sense and are not allowed.

If you make sure the weights are not lower than a certain value, e.g., 0.1, then your example runs through just fine.

I will eventually release a new version of the package that warns about weights that small.

kollerma commented 1 year ago

New version of the package is on CRAN.