Fix PLNnetwork crash issues from glasso

[rfriedman22]

The problem

Short version

When some sparsity parameters are too small, PLNnetwork() crashes completely and returns nothing. It can take quite a while to reach the crash point and there is not a straightforward way to select sparsities or the min_ratio to avoid this issue. All that you get back is an error message like the following:

Error in value[[3L]](cond) : 
  Error in `colnames<-`(`*tmp*`, value = c("ENSDARG00000044056", "ENSDARG00000012763", : attempt to set 'colnames' on an object with less than two dimensions

Similarly, even if the PLNnetworkfamily fits successfully, the same issue can occur when running stability_selection(), which gives a slightly different error:

Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'x' in selecting a method for function 'as.matrix': invalid 'nrow' value (too large or NA)

Long version

When the sparsity for a PLNnetworkfit is sufficiently small, glassoFast diverges on the first iteration of VEM and returns NAs. This causes a PLNnetworkfit to break out of the loop before Omega is set:

  while (!cond) {
    iter <- iter + 1
    if (config$trace > 1) cat("", iter)
    ## CALL TO GLASSO TO UPDATE Omega
    glasso_out <- glassoFast::glassoFast(private$Sigma, rho = self$penalty * self$penalty_weights)
    if (anyNA(glasso_out$wi)) break
    private$Omega <- args$params$Omega <- Matrix::symmpart(glasso_out$wi)
    ...
  }

This doesn't stop myPLN$optimize() from running for the rest of the PLNnetworkfamily. However, it fails during post-treatments because each PLNnetworkfit inherits its post-treatment method from PLNfit, which errors at the following line because Omega is NULL:

rownames(private$Omega) <- colnames(private$Omega) <- colnames(responses)

This causes the entire call to PLNnetwork to crash without returning anything, even if the rest of the network family is fine.

This issue can also come up when running myPLN$stability_selection(), even if it does not cause an issue when fitting the PLNnetworkfamily on the full data. Here, the error is caused when calculating the stability across samples. Specifically, as.matrix(model$latent_network("support")) errors because the latent network is NULL and thus cannot be converted to a matrix.

The solution

I decided to modify the PLNfamily$postTreatment() method so that the call to each model's post-treatment is wrapped in a tryCatch. If it catches an error for model i, PLNfamily$postTreatment() emits a warning and then truncates the family to only include models 1 to i - 1.

Then, I extended the postTreatment() method to PLNnetworkfamily so that after calling the super-method, it also truncates the list of penalties accordingly.

For the stability_selection() in PLNnetworkfamily and ZIPLNnetworkfamily, if Omega is NULL on a subsample, then I interpret the entire network to be unstable at that sparsity, so I return a matrix of zeros.

I'm open to other solutions -- what matters to me is that I am able to run PLNnetwork() and stability_selection() without a complete crash midway through the grid of sparsities.

Validation

I reinstalled a local build of the package and was able to call PLNnetwork() with the exact same data and parameters and get a result where it was previously crashing. I confirmed that the PLNnetworkfamily is correctly truncated and did not see other errors when Iprinted it. I was also able to call PLNnetwork() on a second dataset where I only had an issue with stability_selection() (and not fitting on the whole data) and was able to get a result when it was previously crashing.

I don't run PLNPCA() so I haven't checked or modified PLNPCAfamily which also inherits from PLNfamily.