Errors are caused for certain models when running `fit_mgx_model()`

[svteichman]

Specifically, there are examples (see some cases here) in which the call to raoBust::gee_test calls geeasy::geelm() which creates the error "Error in .solve.checkCond(a, tol) : 'a' is computationally singular, rcond(a)=0".

@svteichman to try to figure out why these errors are happening.

[svteichman]

Notes on this issue:

• when the model does converge for the examples in the linked issue, the beta values are not unusually large
• when the model fails to converge, it runs into case in geeasy:::update_beta where the magnitude of the beta values gets very large, and eventually the hessian has a very large magnitude and cannot be inverted, leading to the error message Error in .solve.checkCond(a, tol): 'a' is computationally singular, ...

[zbartolek]

Hi Sarah and Amy, I think this thread corresponds to the issue we talked about on Monday. Here are a few examples where I've observed it happen with the all the new updates to the model and data that were done more recently. Right now this happens in ~5% of all instances that we run. From what I remember previously, it might still be perfect or quasi-perfect separation issue in the data that I think previously we addressed by skipping over those instances. I attached a few examples here, and the exact error message I get from these. Thanks a lot for all the help!

For example G1_K00052.csv

Error in .solve.checkCond(a, tol) : 
  'a' is computationally singular, rcond(a)=1.04759e-16
In addition: Warning message:
glm.fit: fitted rates numerically 0 occurred 

For example G1_K00477.csv

Error in .solve.checkCond(a, tol) : 
  'a' is computationally singular, rcond(a)=9.65544e-17
In addition: Warning messages:
1: glm.fit: algorithm did not converge 
2: glm.fit: fitted rates numerically 0 occurred 

For example G1_K00515.csv

Error in .solve.checkCond(a, tol) : 
  'a' is computationally singular, rcond(a)=0

For example G1_K00548.csv

Error in .solve.checkCond(a, tol) : 
  'a' is computationally singular, rcond(a)=1.95429e-16
In addition: Warning messages:
1: glm.fit: algorithm did not converge 
2: glm.fit: fitted rates numerically 0 occurred

For example G1_K00588.csv

Error in .solve.checkCond(a, tol) : 
  'a' is computationally singular, rcond(a)=0

Here is the code I use to test these:

library(plyr)
library(dplyr)

library(tidyr)
library(stringr)
library(data.table)
library(ggplot2)
library(splitstackshape)
library(gtools)
library(lsr)
library(tictoc)

library(here) # for file paths
library(R.utils)

remotes::install_github("statdivlab/raoBust")
remotes::install_github("statdivlab/enviromtx")
library(enviromtx)
library(geeasy)

# Load examples to test:
pair_data_funct <- read.csv("~/Desktop/G1_K00052.csv")
pair_data_funct <- read.csv("~/Desktop/G1_K00477.csv")
pair_data_funct <- read.csv("~/Desktop/G1_K00515.csv")
pair_data_funct <- read.csv("~/Desktop/G1_K00548.csv")
pair_data_funct <- read.csv("~/Desktop/G1_K00588.csv")

temp <- pair_data_funct$temp
NO3_NO2_uM <- pair_data_funct$NO3_NO2_uM
Mn_nM <- pair_data_funct$Mn_nM

##### Run the Model:
try({out_i_k <- fit_mgx_model(
  yy = pair_data_funct$counts_sum,
  xstar = pair_data_funct$xstar,
  xx = pair_data_funct$xx,
  formula = ~ temp, # ~temperature + Fe_nM + NO3_NO2_uM + PO4_uM + SiO4_uM + Zn_nM + Cd_pM + Ni_nM + Cu_nM + Pb_pM + Mn_nM + Co_pM #Environemntal co-variates
  enviro_df = cbind(temp, NO3_NO2_uM, Mn_nM),
  replicates = pair_data_funct$sample_name,
  wts = NULL,
  replace_zeros = "minimum")})

[svteichman]

Hi @zbartolek I think this is still a trickier issue, and I unfortunately don't have the bandwidth in the next few weeks to take it on. However, if you remind me again after Nov 5, I'm happy to take another look at this!

[zbartolek]

Hi Sarah, congrats on your defense!!!
Just following up on this issue.

[svteichman]

Thanks for the reminder @zbartolek, I'll see what I can do about these tests that run into this error!

[svteichman]

Hi @zbartolek, I think I've solved this problem in most cases! The problem right now is that in all 5 of your examples, the gee solver is failing due to computational reasons. What I've done in these cases is replaced the gee solver with a glm solver (this is a simpler algorithm and doesn't account for clustering between replicates, so we have to deal with this manually). When I do this, 4 out of the 5 examples work, and the 5th still fails. I think that this is going to probably be the best we can do, but hopefully it's only a small set of examples that still fail in this case. Before re-running anything, please re-install both raoBust and enviromtx (and you likely will need to close your R session and restart another one after the updated install of these packages to make sure that you have the new versions).

However, there is a bit of a complication with the fix that I added. In order to run the robust score tests (which I believe that Amy has recommended to use for inference), you're going to need to input an estimated within-cluster correlation coefficient for all of the cases in which the gee solver failed and the glm solver needs to be used. Amy and I think that a good estimate for this would be the average estimate from the datasets in which it can be estimated, this would be the average correlation:alpha value across your results from running fit_mgx_model(). Based on all of this, here is my recommended workflow:

1. Run all tests as you have in your current code
2. Remove results for all tests that fail (fail as in try({out_i_k <- fit_mgx_model...}) gives an error and doesn't produce results, there are tests that will not fail but will still print error messages to the console), and let us know if this is a still a large amount (I'm hoping that it'll be far less than before).
3. Loop through your remaining results and compute the average value of correlation:alpha for all of the results out_i_k for which it is not NA. You could do this with something like the following:

alpha_vec <- rep(NA, n_test) # n_test is the number of tests you're running
for (i in 1:n_test) {
    alpha_vec <- result[3, 1] # where result is the result of fit_mgx_model for test data i 
}
estimated_alpha <- mean(alpha_vec, na.rm = TRUE)

4. Identify all the results for which the column for Robust Score p in the results is set to NA (these are the results in which the gee solver failed). Re-run these tests, but now add an additional argument to fit_mgx_model():

out_i_k <- fit_mgx_model(
  yy = pair_data_funct$counts_sum,
  xstar = pair_data_funct$xstar,
  xx = pair_data_funct$xx,
  formula = ~ temp, 
  enviro_df = cbind(temp, NO3_NO2_uM, Mn_nM),
  replicates = pair_data_funct$sample_name,
  wts = NULL,
  replace_zeros = "minimum",
  cluster_corr_coef = estimated_alpha) # where you use the estimated_alpha from above

Please let me know if any of this is unclear or if you run into any issues when doing it! I'd also be curious about the results of running summary(alpha_vec) after you generate that vector as described above, and what proportion of your tests still fail.

[zbartolek]

Hi Sarah, thanks a lot for this! I will try it out tomorrow and let you know if I run into any issues!

[zbartolek]

Hi Sarah, I have a few clarifying questions for you.

I re-downloaded the packages and opened a new R session. I'm still having trouble running the new code on test cases. This is what I do:

1. Original code: (run on the last example case (G1_K00588.csv)

library(plyr)
library(dplyr)

library(tidyr)
library(stringr)
library(data.table)
library(ggplot2)
library(splitstackshape)
library(gtools)
library(lsr)
library(tictoc)

library(here) # for file paths
library(R.utils)

remotes::install_github("statdivlab/raoBust")
remotes::install_github("statdivlab/enviromtx")
library(enviromtx)
library(raoBust)
library(geeasy)

# Load examples to test:
pair_data_funct <- read.csv("~/Desktop/G1_K00052.csv")
pair_data_funct <- read.csv("~/Desktop/G1_K00477.csv")
pair_data_funct <- read.csv("~/Desktop/G1_K00515.csv")
pair_data_funct <- read.csv("~/Desktop/G1_K00548.csv")
pair_data_funct <- read.csv("~/Desktop/G1_K00588.csv")

temp <- pair_data_funct$temp
NO3_NO2_uM <- pair_data_funct$NO3_NO2_uM
Mn_nM <- pair_data_funct$Mn_nM

##### Run the Model:
try({out_i_k <- fit_mgx_model(
  yy = pair_data_funct$counts_sum,
  xstar = pair_data_funct$xstar,
  xx = pair_data_funct$xx,
  formula = ~ temp, # ~temperature + Fe_nM + NO3_NO2_uM + PO4_uM + SiO4_uM + Zn_nM + Cd_pM + Ni_nM + Cu_nM + Pb_pM + Mn_nM + Co_pM #Environemntal co-variates
  enviro_df = cbind(temp, NO3_NO2_uM, Mn_nM),
  replicates = pair_data_funct$sample_name,
  wts = NULL,
  replace_zeros = "minimum")})

Error message:

Error in .solve.checkCond(a, tol) : 
  'a' is computationally singular, rcond(a)=0
There were 50 or more warnings (use warnings() to see the first 50)

2. Trying out the new cluster_corr_coef = estimated_alpha with a fake estimated_alpha value (also on G1_K00588.csv):

###### If it fails, try to solve with an estimated alpha:
estimated_alpha <- 0.2# Fake value for testing only. 
try({out_i_k <- fit_mgx_model(
  yy = pair_data_funct$counts_sum,
  xstar = pair_data_funct$xstar,
  xx = pair_data_funct$xx,
  formula = ~ temp, # ~temperature + Fe_nM + NO3_NO2_uM + PO4_uM + SiO4_uM + Zn_nM + Cd_pM + Ni_nM + Cu_nM + Pb_pM + Mn_nM + Co_pM #Environemntal co-variates
  enviro_df = cbind(temp, NO3_NO2_uM, Mn_nM),
  replicates = pair_data_funct$sample_name,
  wts = NULL,
  replace_zeros = "minimum",
  cluster_corr_coef = estimated_alpha)})

I still get this error message:

Error in .solve.checkCond(a, tol) : 
  'a' is computationally singular, rcond(a)=0
There were 50 or more warnings (use warnings() to see the first 50)

I get this with the other examples too... hopefully the installation is fine this time, is there a way to check that?
I might be missing something very obvious here!

Before trying this on a full dataset, I also want to make sure I understand everything first.

So for the second step, when you say remove results for all tests that fail. From what I understand, we still expect the instances that have failed previously with gee to fail until this point too, is that correct?

Next, I calculate the average correlation:alpha for those that did not fail and use that to re-run the failed examples with the glm and estimated alpha. So just to be clear, we are still using gee for most cases and then redoing the failed ones with glm and the calculated average alpha. The reason I'm asking is that I could use the output I already have from gee (where I know which cases failed and didn't fail) and search for failed cases, calculate the average correlation:alpha from the ones that ran well, and re-run the failed ones only. Is this a valid thing to do?

Last question, when calculating average correlation:alpha, can I calculate this from the entire dataset of a particular cruise, or should this be specific to each pair that I test?

Thanks a lot!
Zinka

[zbartolek]

Just to check that the enviromtx updated on my computer, I ran this
enviromtx::fit_mgx_model

> enviromtx::fit_mgx_model
function (yy, xstar, xx, replicates = NULL, formula = NULL, enviro_df = NULL, 
    wts = NULL, replace_zeros = "minimum", use_jack_se = FALSE, 
    cluster_corr_coef = NULL) 
{
    if (replace_zeros == "minimum") {
        xx[xx == 0] <- min(xx[xx > 0])
    }
    else if (is.numeric(replace_zeros)) {
        if (replace_zeros == 0) 
            stop("You've told me to replace zeroes with zero!")
        xx[xx == 0] <- replace_zeros
    }
    if (any(xx == 0)) 
        stop("There are zeros in xx and you haven't told me what to do with them!")
    if (replace_zeros == "minimum") {
        xstar[xstar == 0] <- min(xstar[xstar > 0])
    }
    else if (is.numeric(replace_zeros)) {
        if (replace_zeros == 0) 
            stop("You've told me to replace zeroes with zero!")
        xstar[xstar == 0] <- replace_zeros
    }
    if (any(xstar == 0)) 
        stop("There are zeros in xstar and you haven't told me what to do with them!")
    predictor <- log(xstar/xx)
    if (!is.null(formula) & !is.null(enviro_df)) {
        my_formula <- lazyeval::f_new(lhs = quote(yy), rhs = lazyeval::as_call(paste("predictor +", 
            as.character(formula)[2])))
        my_df <- dplyr::bind_cols(tibble::tibble(yy, xstar, xx, 
            predictor), enviro_df)
    }
    else if (xor(is.null(formula), is.null(enviro_df))) {
        stop("One of `formula` or `enviro_df` is missing. Please correct this and try again.")
    }
    else if (is.null(formula) & is.null(enviro_df)) {
        my_df <- tibble::tibble(yy, xstar, xx, predictor)
        my_formula <- yy ~ predictor
    }
    else {
        stop("Ergh, there is an error and it is Amy's fault.")
    }
    if (any(apply(my_df, 2, is.infinite))) {
        print(my_df, n = nrow(my_df))
        stop("Infinities in my_df? Amy's fault.")
    }
    if (is.null(wts)) {
        wts <- rep(1L, nrow(data.frame(my_df)))
    }
    my_df$wts <- wts
    if (is.null(replicates)) {
        raoBust_out <- raoBust::glm_test(formula = my_formula, 
            offset = log(xx), family = stats::poisson(link = "log"), 
            data = my_df, weights = wts)
    }
    else {
        if (!all(wts == 1L)) {
            warning("Run this by Amy; not sure what this is doing off-the-cuff")
        }
        my_df$id <- replicates
        raoBust_out <- raoBust::gee_test(formula = my_formula, 
            offset = log(xx), family = stats::poisson(link = "log"), 
            id = id, data = my_df, use_jack_se = use_jack_se, 
            cluster_corr_coef = cluster_corr_coef)
    }
    raoBust_out[-1, ]
}
<bytecode: 0x1682c3190>
<environment: namespace:enviromtx>

This seems to be the same as your new update on github, unless I'm missing something!

[svteichman]

Great, thanks for this update!

I think one thing that would help (if you aren't doing this already in your code), is saving each result to a different element in a list. Something like the following:

files <- list.files(file path where csvs are stored) # this should be a vector of all the csv's you're working with
res <- vector(mode = "list", length = length(files)

for (f in 1:length(files)) {
   pair_data_funct <- read.csv(files[f])
   temp <- pair_data_funct$temp
   NO3_NO2_uM <- pair_data_funct$NO3_NO2_uM
   Mn_nM <- pair_data_funct$Mn_nM

   res[[f]] <- try({out_i_k <- fit_mgx_model(
                             yy = pair_data_funct$counts_sum,
                             xstar = pair_data_funct$xstar,
                             xx = pair_data_funct$xx,
                             formula = ~ temp, 
                             enviro_df = cbind(temp, NO3_NO2_uM, Mn_nM),
                             replicates = pair_data_funct$sample_name,
                             wts = NULL,
                             replace_zeros = "minimum")})
}

This way, when you're done you can loop through the list res to see which tests generated results and which ones failed (if they failed, the element res[[f]] for that dataset will be NULL. Something like:

passed <- !unlist(lapply(res, is.null))

Then, passed will be a vector of the length of your number of tests, with entry TRUE for that test if it has estimate results. Then you'll also want a vector that tells you whether or not the results include robust score test results (because using the glm solver won't give you robust score test results unless we give it the correlation coefficient).

robust_score_run <- unlist(lapply(res, function(x) {!(is.null(x)) && !(is.na(x[2, 4]))}))

Based on these results, any test with index i that has passed[i] = TRUE and robust_score_run[i] = TRUE has passed and has successfully run the robust score test, a test that has passed[i] = TRUE and robust_score_run[i] = FALSE has estimation results but will need to be rerun to get robust score test results, and a test that has passed[i] = FALSE can't be optimized with this method currently.

Using these results, you could then loop through the results res list for all tests for which robust_score_run = TRUE, and calculate the average estimated correlation coefficient for those tests. Then, you can rerun just the tests for which passed = TRUE and robust_score_run = FALSE, and rerun those adding in the cluster_corr_coef argument to fit_mgx_model.

It's possible that a few of these robust score tests could then fail due to something in the robust score test function, but I'm hoping that will happen pretty rarely if so.

---

This gets me to your first point in your comment above. My guess is that when you run the code in step 1), you are getting these error messages (which are still being thrown from the function when it tries to run the gee solver), but I'm hoping that the glm solver will work, so results will still be generated. So for most of the tests that were previously failing, I'm hoping that you'll have passed = TRUE for these tests because they are generating a result table.

"Next, I calculate the average correlation:alpha for those that did not fail and use that to re-run the failed examples with the glm and estimated alpha. So just to be clear, we are still using gee for most cases and then redoing the failed ones with glm and the calculated average alpha. The reason I'm asking is that I could use the output I already have from gee (where I know which cases failed and didn't fail) and search for failed cases, calculate the average correlation:alpha from the ones that ran well, and re-run the failed ones only. Is this a valid thing to do?"

Yes, totally valid!

"Last question, when calculating average correlation:alpha, can I calculate this from the entire dataset of a particular cruise, or should this be specific to each pair that I test?"

I think doing it for the entire dataset of a particular cruise sounds like the right call to me.

Basically what's happening here is that we can't optimize some of these models for certain datasets with gee solvers because they run into computational issues. A gee solver will give us parameter estimates and estimated within-cluster correlations for that dataset, which we then give to the robust score test function to get a robust score test p-value. So, when the gee fails, we instead use a glm which will give us parameter estimates (which won't be quite as accurate for a given test but are still on average accurate), but doesn't estimate the within-cluster correlation for that dataset. So for those settings, we are hoping that the average estimated correlation coefficient for all datasets for that cruise is a good guess for this specific dataset, so we can give that as our guess to the robust score test function.

[adw96]

@zbartolek related to my question from #16 -- re: the exact code you were running and if it included adjustment covariates -- running

fit_mgx_model(yy = pair_data_funct$counts_sum,
                             xstar = pair_data_funct$xstar,
                             xx = pair_data_funct$xx,
                             formula = ~ temp, 
                             enviro_df = cbind(temp, NO3_NO2_uM, Mn_nM),
                             replicates = pair_data_funct$sample_name,
                             wts = NULL,
                             replace_zeros = "minimum")

won't adjust for manganese and nitrate+nitrite. You need to also add them to the formula, eg

fit_mgx_model(yy = pair_data_funct$counts_sum,
                             xstar = pair_data_funct$xstar,
                             xx = pair_data_funct$xx,
                             formula = ~ temp + NO3_NO2_uM + Mn_nM, 
                             enviro_df = cbind(temp, NO3_NO2_uM, Mn_nM),
                             replicates = pair_data_funct$sample_name,
                             wts = NULL,
                             replace_zeros = "minimum")

Happy to discuss! I'll flag this in your text too so we are sure that the code and manuscript align.