Fatal model failure

[hyjforesight]

Hello CellRank,
I'm running the CytoTRACE kernel. Everything works fine until plotting the gene trends.
Appreciate if you could help me with this issue.
Thanks!
Best,
YJ

model = cr.ul.models.GAM(adata)
adata.varm['terminal_lineage_drivers'].sort_values(by="Pit_corr", ascending=False)
	Pit_corr	Pit_pval	Pit_qval	Pit_ci_low	Pit_ci_high
Pclaf	0.509537	3.796250e-136	5.135187e-132	0.475956	0.541634
Gsta1	0.403946	7.050554e-80	4.768642e-76	0.366171	0.440392
Esco2	0.401800	6.000396e-79	2.705578e-75	0.363951	0.438325
Pigr	0.387159	8.471574e-73	2.864875e-69	0.348813	0.424210
Gkn3	0.380823	3.060182e-70	8.279016e-67	0.342269	0.418097
...	...	...	...	...	...
2410006H16Rik	-0.214614	5.875677e-22	7.033653e-20	-0.256506	-0.171917
Ly6c1	-0.282888	8.739462e-38	2.955468e-35	-0.323167	-0.241585
Rhox5	-0.302470	2.729246e-43	1.190920e-40	-0.342214	-0.261646
Ly6a	-0.304436	7.212288e-44	3.484308e-41	-0.344124	-0.263662
Ly6d	-0.333813	4.396291e-53	4.247759e-50	-0.372631	-0.293829

cr.pl.gene_trends(adata, model, data_key="X", genes=['Pclaf', 'Gsta1'], time_key="ct_pseudotime", same_plot=True, hide_cells=True, n_test_points=200, backward=False, n_jobs=8)
ValueError                                Traceback (most recent call last)
~\anaconda3\envs\HYJ_py38\lib\site-packages\cellrank\ul\models\_pygam_model.py in fit()
    207         try:
--> 208             self.model.fit(self.x, self.y, weights=self.w, **kwargs)
    209             return self

~\anaconda3\envs\HYJ_py38\lib\site-packages\pygam\pygam.py in fit()
    894         # validate data
--> 895         y = check_y(y, self.link, self.distribution, verbose=self.verbose)
    896         X = check_X(X, verbose=self.verbose)

~\anaconda3\envs\HYJ_py38\lib\site-packages\pygam\utils.py in check_y()
    224         if np.any(np.isnan(link.link(y, dist))):
--> 225             raise ValueError('y data is not in domain of {} link function. ' \
    226                              'Expected domain: {}, but found {}' \

ValueError: y data is not in domain of log link function. Expected domain: [0.0, inf], but found [-0.26, 2.31]

The above exception was the direct cause of the following exception:

RuntimeError                              Traceback (most recent call last)
~\anaconda3\envs\HYJ_py38\lib\site-packages\cellrank\ul\models\_base_model.py in wrapper()
     64                     instance.reraise()
---> 65                 return wrapped(*args, **kwargs)
     66             except Exception as e:  # noqa: B902

~\anaconda3\envs\HYJ_py38\lib\site-packages\cellrank\ul\models\_pygam_model.py in fit()
    210         except Exception as e:  # noqa: B902
--> 211             raise RuntimeError(
    212                 f"Unable to fit `{type(self).__name__}` for gene "

RuntimeError: Unable to fit `GAM` for gene `'Pclaf'` in lineage `'Pit'`. Reason: `y data is not in domain of log link function. Expected domain: [0.0, inf], but found [-0.26, 2.31]`

The above exception was the direct cause of the following exception:

RuntimeError                              Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_19192/2586172927.py in <module>

---->  cr.pl.gene_trends(adata, model, data_key="X", genes=['Pclaf', 'Gsta1'], time_key="ct_pseudotime", same_plot=True, hide_cells=True, n_test_points=200, backward=False, n_jobs=8)

~\anaconda3\envs\HYJ_py38\lib\site-packages\cellrank\ul\_utils.py in _genesymbols(wrapped, instance, args, kwargs)
    284         use_raw=kwargs.get("use_raw", False),
    285     ):
--> 286         return wrapped(*args, **kwargs)

~\anaconda3\envs\HYJ_py38\lib\site-packages\cellrank\pl\_gene_trend.py in gene_trends(adata, model, genes, lineages, backward, data_key, time_key, time_range, transpose, callback, conf_int, same_plot, hide_cells, perc, lineage_cmap, abs_prob_cmap, cell_color, cell_alpha, lineage_alpha, size, lw, cbar, margins, sharex, sharey, gene_as_title, legend_loc, obs_legend_loc, ncols, suptitle, return_models, n_jobs, backend, show_progress_bar, figsize, dpi, save, plot_kwargs, **kwargs)
    197     models = _create_models(model, genes, lineages)
    198 
--> 199     all_models, models, genes, lineages = _fit_bulk(
    200         models,
    201         _create_callbacks(adata, callback, genes, lineages, **kwargs),

~\anaconda3\envs\HYJ_py38\lib\site-packages\cellrank\pl\_utils.py in _fit_bulk(models, callbacks, genes, lineages, time_range, parallel_kwargs, return_models, filter_all_failed, **kwargs)
    483     logg.info("    Finish", time=start)
    484 
--> 485     return _filter_models(
    486         models, return_models=return_models, filter_all_failed=filter_all_failed
    487     )

~\anaconda3\envs\HYJ_py38\lib\site-packages\cellrank\pl\_utils.py in _filter_models(models, return_models, filter_all_failed)
    532                     model, FailedModel
    533                 ), f"Expected `FailedModel`, found `{type(model).__name__!r}`."
--> 534                 model.reraise()
    535 
    536     if not np.all(modelmask.values):

~\anaconda3\envs\HYJ_py38\lib\site-packages\cellrank\ul\models\_base_model.py in reraise(self)
   1349         """Raise the original exception with additional model information."""
   1350         # retain the exception type and also the original exception
-> 1351         raise type(self._exc)(f"Fatal model failure `{self}`.") from self._exc
   1352 
   1353     def _get_colors(

RuntimeError: Fatal model failure `<FailedModel[origin=GAM[gene='Pclaf', lineage='Pit', model=GammaGAM(callbacks=[Deviance(), Diffs()], fit_intercept=True, max_iter=2000, scale=None, terms=s(0), tol=0.0001, verbose=False)]]>`.

RuntimeError: Fatal model failure `<FailedModel[origin=GAM[gene='Pclaf', lineage='Pit', model=GammaGAM(callbacks=[Deviance(), Diffs()], fit_intercept=True, max_iter=2000, scale=None, terms=s(0), tol=0.0001, verbose=False)]]>`.

Versions:

scanpy==1.8.2 anndata==0.7.8 umap==0.5.2 numpy==1.21.5 scipy==1.7.3 pandas==1.3.5 scikit-learn==1.0.2 statsmodels==0.13.1 python-igraph==0.9.9 pynndescent==0.5.5
scvelo==0.2.4 scanpy==1.8.2 anndata==0.7.8 loompy==3.0.6 numpy==1.21.5 scipy==1.7.3 matplotlib==3.5.1 sklearn==1.0.2 pandas==1.3.5
cellrank==1.5.1 scanpy==1.8.2 anndata==0.7.8 numpy==1.21.5 numba==0.55.0 scipy==1.7.3 pandas==1.3.5 pygpcca==1.0.3 scikit-learn==1.0.2 statsmodels==0.13.1 python-igraph==0.9.9 scvelo==0.2.4 pygam==0.8.0 matplotlib==3.5.1 seaborn==0.11.2

...

[michalk8]

Hi @hyjforesight ,

the error shows that ValueError: y data is not in domain of log link function. Expected domain: [0.0, inf], but found [-0.26, 2.31] and we use by default log link function, which has to have [0.0, inf].
I would assume that the data in adata.X is somehow normalized (therefore the negative values), so I'd recommend either to change the link function that can support this domain (e.g., identity + gaussian distribution when initialized the model) or working on normalized/scaled data (e.g. in adata.raw, if present).

[hyjforesight]

Hello @michalk8
Thanks for the response. As you recommended, use model = cr.ul.models.GAM(adata, distribution='gaussian', link='identity') can solve it. This error was caused by I running the regressed and scaled data. I used the raw data and error was also solved. Thanks!

[nynkekatinka]

Hi @michalk8 , I just wanted to let you know that I got the same error from running the "getting started with CellRank" tutorial. I saw the data is also normalized at the beginnen, so using model = cr.models.GAM(adata, distribution='gaussian', link='identity') allowed me to plot a similar ish version of the plots in the tutorial. Just thought you might like to know this happens and maybe make a note of it in the tutorial :)