covaruber / sommer

36 stars 23 forks source link

Can't get mmec to run #41

Closed Ax3man closed 2 years ago

Ax3man commented 2 years ago

I'm trying to see if mmec would be useful for a project I am working on.

If I run this (slightly modified) example:

data(DT_gryphon)
DT <- DT_gryphon
Ag <- A_gryphon

## mmec uses the inverse of the relationship matrix
Agi <- solve(Ag + diag(1e-4,ncol(Ag),ncol(Ag)))
mix1b <- mmec(
  BWT~1,
  random=~vsc(isc(ANIMAL),Gu=Agi),
  rcov=~units,
  data=DT
)

I get the following error:

Error in if (class(Gu) %!in% c("dgCMatrix")) { : 
the condition has length > 1

This may be because of R 4.2.0 any if condition of length > 1 is now an error, and nowadays matrices have both class matrix and array. To avoid this error, check the class of an object using inherits. So in this case I think this would solve the issue:

if (!inherits(Gu, "dgCMatrix")) { ...... }

Looking at the source I can see that I need to pass a sparse matrix, but the above generates its error just before your stop.

covaruber commented 2 years ago

Dear @Ax3man since the mmec relies in sparse operations you need to make sure that Agi is also of class sparse matrix. Please see the line of code I added:

library(sommer)
data(DT_gryphon)
DT <- DT_gryphon
Ag <- A_gryphon

## mmec uses the inverse of the relationship matrix
Agi <- solve(Ag + diag(1e-4,ncol(Ag),ncol(Ag)))
Agi[1:4,1:4]

Agi <- as(Agi, Class="sparseMatrix") ## this is the missing part

mix1b <- mmec(
  BWT~1,
  random=~vsc(isc(ANIMAL),Gu=Agi),
  rcov=~units,  tolParConvNorm = .001,
  data=DT
)

I modify the tolParConvNorm parameter because the mmec function is based on 2 stopping criteria, the norm of VC change (explained in the documentation) and the log-likelihood. Looking at the example seems like is better to use LL in this case so I set the tolParConvNorm to a smaller value than the default.

Cheers, Eduardo

Ax3man commented 2 years ago

Hi Eduardo,

Thanks! Yes, I figured out the sparse matrix issue. But this line:

if(class(Gu) %!in% c("dgCMatrix")){
      stop("Gu matrix is not of class dgCMatrix. Please correct \n", call. = TRUE )
}

which is supposed to tell the user that, will not work on newer R versions when the user passes a matrix instead of a sparse matrix. Instead, it will give the uninformative error the condition has length > 1.

Thanks so much for the very useful package, Wouter.

covaruber commented 2 years ago

Oh OK, I didn't realize of that issue, I will adopt your suggestion and use inherit. I will leave the ticket open until I push the change tonight :) Thanks Wouter!

Cheers, Eduardo

covaruber commented 2 years ago

Change has been pushed to github, thanks for the report :)

Cheers, Eduardo