Documenting: interesting behaviour susie in small cases

[william-denault]

@pcarbo @stephens999

To reproduce this example your first need to have the susie package from the commit of Aug 2025, the data used in this example is now modified.

remotes::install_github("stephenslab/susieR@b5f5f7bddce71c1dfd2f401a24985338292e4f31")

Now we can look at the problematic example from microglia

library(susieR)
data(data_small)
y <- data_small$y
X <- data_small$X
dim(X)



res_susie <- susie(X  ,y,L = 10,verbose = FALSE)
res_susie$sets$cs
 
susie_plot(res_susie,y = "PIP")


ypred <- predict(res_susie,X)
pve <- 1 - drop(res_susie$sigma2/var(y))
round(100*pve,digits = 1)

This returns 10 CS and explains 99.9% of the variance

res_susie$sets$cs
$L1
[1] 4919 4920

$L2
[1] 455 456 457 462
 .....

Something interesting happen if one remove SNP in CS directly from the data set

res_susie1 <- susie(X[,-c(450:470)]  ,y,L = 10,verbose = FALSE)
 res_susie1$sets$cs
$L1
[1] 4898 4899

This is quite strange.

[pcarbo]

@william-denault I was able to reproduce this result.

I agree this result is quite strange; it shows that susie is quite unpredictable at smaller sample sizes.

The X matrix is not the original genotype matrix, so I can't check this, but an additional issue could be that the MAFs of some of the SNPs are quite small.

It is reassuring though that when you rerun with small = TRUE, the results are very similar!

[pcarbo]

I'll also note that when you do standardize = FALSE and small = FALSE, the results look much better. Regardless, we would expect it to work well with standardize = TRUE also.

[gaow]

Using unmappable_effects="ash" in version 0.15.3, no CS is found after a long time to converge, although given the small sample size n=47 it looks reasonable.

res_susie <- susie(X  ,y,L = 10, unmappable_effects="ash", verbose = TRUE)

unmappable_effects="inf" is messier, see below,