fit_nullmodel() fails with “Not a matrix” when sparse GRM is supplied

[komais]

Hello,

 obj_nullmodel <- fit_nullmodel(pheno ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10,
    data = phenotype, kins = sgrm, use_sparse = TRUE, kins_cutoff = 0.022, id = "FID",maxiter=10,
    family = binomial(link = "logit"), verbose = TRUE, groups = NULL
) 

got an error：

Fixed-effect coefficients:
 (Intercept)          PC1          PC2          PC3          PC4          PC5
 -5.87561278 472.44855848 -93.49517238  -1.13570207   8.72329401  -6.69928427
         PC6          PC7          PC8          PC9         PC10
 -2.93669950  -0.06602734  -1.28533097   3.12009658   1.54106524

Iteration  1 :
Error in glmmkin.brent(fit0, kins[[1]], method = method, maxiter = maxiter,  :
  Not a matrix.
In addition: Warning messages:
1: glm.fit: algorithm did not converge
2: Average Information REML not converged, refitting model using Brent method...

i have sgrm looks like

str(sgrm)
Formal class 'dsCMatrix' [package "Matrix"] with 7 slots
  ..@ i       : int [1:347] 0 0 1 2 2 3 4 4 5 6 ...
  ..@ p       : int [1:344] 0 1 3 4 6 7 9 10 12 13 ...
  ..@ Dim     : int [1:2] 343 343
  ..@ Dimnames:List of 2
  .. ..$ : chr [1:343] "0_HG00475" "0_HG00542" "0_HG00581" "0_HG00607" ...
  .. ..$ : chr [1:343] "0_HG00475" "0_HG00542" "0_HG00581" "0_HG00607" ...
  ..@ x       : num [1:347] 0.501 0.0626 0.4098 0.4996 0.0525 ...
  ..@ uplo    : chr "U"
  ..@ factors : list()

how can i fix this error ?
thank you

[ArdaCet]

Hi @komais!,

I encountered a similar issue when running a STAAR null model for binary traits. I am wondering if you found a workaround, or if someone reached out to you with a solution?

I am providing some information below for the developers @xihaoli @zilinli1988 @LVMEHINOVIC:

Here is my model:

obj_nullmodel <- fit_nullmodel(Disease_Status~Sex+Cohort1+Cohort2+Cohort3+Cohort4+
C1+C2+C3+C4+C5+C6+C7+C8+C9+C10,
data=phenotype, kins=sgrm, use_sparse=TRUE, kins_cutoff=0.022, id="sample.id" ,family=binomial(link = "logit"), verbose=TRUE)

The sparse GRM matrix structure:

str(matrix)
Loading required package: Matrix
Formal class 'dsCMatrix' [package "Matrix"] with 7 slots
  ..@ i       : int [1:244976] 0 0 1 2 2 3 2 3 4 5 ...
  ..@ p       : int [1:18849] 0 1 3 4 6 9 10 12 15 19 ...
  ..@ Dim     : int [1:2] 18848 18848
  ..@ Dimnames:List of 2
  .. ..$ : chr [1:18848] "A" "B" "C" "D" ...
  .. ..$ : chr [1:18848] "A" "B" "C" "D" ...
  ..@ x       : num [1:244976] 0.495 0.233 0.505 0.5 0.237 ...
  ..@ uplo    : chr "U"
  ..@ factors : list()

and the error I received after several iterations:

Iteration  500 :
Variance component estimates:
[1] 1 0
Fixed-effect coefficients:
 [1] -27.9864971   0.1367875  28.9970812  29.0402314  30.8960334  28.0911089
 [7] -13.5624838 -12.0874552 -45.7784889 -13.0274918  -7.4289290   2.3963156
[13]   4.9003013   6.0348662  -7.1453006  -2.3972775
Fixed-effect coefficients:
(Intercept)         Sex  Cohort1       Cohort2       Cohort3       Cohort4
-21.6669975   0.1367875  22.6775823  22.7207317  24.5765340  21.7716083
         C1          C2          C3          C4          C5          C6
-13.5624838 -12.0874552 -45.7784889 -13.0274918  -7.4289290   2.3963156
         C7          C8          C9         C10
  4.9003013   6.0348662  -7.1453006  -2.3972775

Iteration  1 :
Error: Not a matrix.
In addition: Warning messages:
1: Variance estimate on the boundary of the parameter space observed, refitting model...
2: Average Information REML not converged, refitting model using Brent method...
Execution halted

[xihaoli]

Hi both,

Thanks for reaching out. STAAR share the statistical (GLMM) framework with GMMAT. I am including @hanchenphd who may provide additional insights in null model fitting.

Best,
Xihao

@

[hanchenphd]

Hi @komais and @ArdaCet,

Convergence issues in logistic mixed models occur more frequently when the sample size is small or when the kinship matrix is close to a diagonal matrix. One issue that I noticed was that in your initial step, you both had very large coefficients, which could imply issues in your data. It is worth noting that the initial step is a regular logistic regression model, and the coefficients are log odds ratios. For example, in @komais's results, the PC1 coefficient was 472, meaning that one unit increase in PC1 corresponds to an odds ratio of exp(472.45)=1.5x10^205. In @ArdaCet's results, the Cohort1, 2, 3, 4 coefficients were all large, Cohort1 effect estimate from a logistic regression model (without any random effects) shows an odds ratio of exp(22.6776)=7x10^9. So the first step to troubleshoot is to double check your cases and controls and see how you can interpret such large odds ratios (usually it means you only have cases in Cohort1, or only have controls in the reference group).

Also both kinship matrices were very sparse, with most off-diagonal numbers being 0, which increased the chance of observing convergence issues. In @komais's results, the sample size was small (N = 343), with 347 non-zero values in total in the kinship matrix, which means only 4 off-diagonal numbers were non-zero. After you figure out the issue with PC1 and PC2, I would still not recommend running a logistic mixed model, because most of your samples were uncorrelated. It would be better to exclude these a few related samples and run a logistic regression model instead, which is more parsimonious.

Best,
Han

[ArdaCet]

Dear @xihaoli and @hanchenphd ,

Thank you very much for your kindness, incredibly fast, and helpful responses. I will check for the imbalance in my cohorts and work on finding other covariates, which will not be biased for one group (case/control). Hopefully, I will find a workaround and correct these high coefficients in the logistic regression model.

Again, thank you for your consideration, time, and the support for the tools.

Best,
Arda

[xihaoli]

Dear @ArdaCet,

Thank you for your reply. Hope these recommendations can resolve your question about null model fitting. Please feel free to let us know if you have any other questions.

Best,
Xihao