SuSiE can crash when run in small samples

[william-denault]

Hello,

As suggested by @stephens999 , I am opening an issue about that.

If you run the script below for a couple of minutes it will crash and output

Error in init_finalize(s) : 
  Residual variance sigma2 must be positive (is your var(Y) zero)

library(susieR
        )

N <- 20
set.seed(123)

res <- list( )
data(N3finemapping)
h=1
for ( i in (length(res)+1):10000){

  L <- sample (1:10, size=1)

  X <- N3finemapping$X
  X <- N3finemapping$X[sample (1:nrow(X), size=N, replace=FALSE),]
  true_pos <- sample( 1:ncol(X), L)
  if( L==1){
    if (var( X[,true_pos])==0) next
    y <- X[, true_pos]
  }else{
    y <- apply( X[, true_pos],1, sum)
  }

  y <- y + rnorm( N, sd=  1*( sd(y)))

  out <-  susie(X,y, L=1 )  
  out$sets
    if(!is.null(out$sets$cs)){

    n_true_cs <-   Reduce("+",sapply(1:length(out$sets$cs), function(k)
        ifelse( length(which(true_pos%in%out$sets$cs[[k]] ))==0, 0,1)
      ))
    n_cs   <-    length(  out$sets$cs  )
    n_effect <- length(true_pos)

      res [[h]] <- c( n_true_cs ,   n_cs,n_effect   )

        print( res[[h]])

      h=h+1
    }

} 

[pcarbo]

@william-denault The error message says to check that the variance of y is zero, and indeed it is zero in your example:

> i
[1] 429
> var(y)
[1] 0

One issue perhaps is that this step has no effect when y has no variance:

y <- y + rnorm( N, sd= 1*( sd(y)))