Issue with multinomial model in pygformula when using categorical covariates with more than 6 levels

[jgonzhijon]

I'm encountering an issue when using categorical covariates with more than 6 levels in a multinomial model with pygformula. The model fitting process produces NaN values, which then prevents bootstrap results from being computed — due to a NoneType error in the results.

This does not happen if I:

• Reduce the number of levels in the categorical variable to 6 or fewer.
• Use one-hot encoding and fit a one-vs-rest model instead of multinomial.

I would prefer to use the multinomial option directly with categorical covariates if possible.

My questions are:
What could be causing the multinomial model to fail when categorical variables have many levels?
Is there a known limitation or workaround in pygformula or statsmodels?
If this can't be fixed cleanly, would using one-vs-rest with one-hot encoding be a valid alternative?

Here is the relevant part of the error output, the warnings repeat across all the execution :

 
/opt/anaconda3/lib/python3.12/site-packages/statsmodels/discrete/discrete_model.py:3028: RuntimeWarning: invalid value encountered in divide
  return eXB/eXB.sum(1)[:,None]
Optimization terminated successfully.
Current function value: nan
Iterations 4
Optimization terminated successfully.
Current function value: 0.214036
Iterations 12
Optimization terminated successfully.
Current function value: 0.458476
Iterations 8
/opt/anaconda3/lib/python3.12/site-packages/statsmodels/discrete/discrete_model.py:3027: RuntimeWarning: overflow encountered in exp
  eXB = np.column_stack((np.ones(len(X)), np.exp(X)))
/opt/anaconda3/lib/python3.12/site-packages/statsmodels/discrete/discrete_model.py:3028: RuntimeWarning: invalid value encountered in divide
  return eXB/eXB.sum(1)[:,None]
Optimization terminated successfully.
Current function value: nan
Iterations 5
/opt/anaconda3/lib/python3.12/site-packages/statsmodels/discrete/discrete_model.py:3027: RuntimeWarning: overflow encountered in exp
  eXB = np.column_stack((np.ones(len(X)), np.exp(X)))
/opt/anaconda3/lib/python3.12/site-packages/statsmodels/discrete/discrete_model.py:3028: RuntimeWarning: invalid value encountered in divide
  return eXB/eXB.sum(1)[:,None]
Optimization terminated successfully.
Current function value: 0.471235
Iterations 8
Optimization terminated successfully.
Current function value: nan
Iterations 4
Optimization terminated successfully.
Current function value: 0.205321
Iterations 10
Optimization terminated successfully.
Current function value: 0.458697
Iterations 8
/opt/anaconda3/lib/python3.12/site-packages/statsmodels/discrete/discrete_model.py:3027: RuntimeWarning: overflow encountered in exp
  eXB = np.column_stack((np.ones(len(X)), np.exp(X)))
/opt/anaconda3/lib/python3.12/site-packages/statsmodels/discrete/discrete_model.py:3028: RuntimeWarning: invalid value encountered in divide
  return eXB/eXB.sum(1)[:,None]
Optimization terminated successfully.
Current function value: nan
Iterations 4
/opt/anaconda3/lib/python3.12/site-packages/statsmodels/discrete/discrete_model.py:3027: RuntimeWarning: overflow encountered in exp
  eXB = np.column_stack((np.ones(len(X)), np.exp(X)))
/opt/anaconda3/lib/python3.12/site-packages/statsmodels/discrete/discrete_model.py:3028: RuntimeWarning: invalid value encountered in divide
  return eXB/eXB.sum(1)[:,None]
Optimization terminated successfully.
Current function value: nan
Iterations 4
Traceback (most recent call last):
  File "/Users/juanito/Documents/Paper1-HSU/HSU-PSID-Data/Main_Analysis/prueba_error.py", line 147, in 
    g.fit()
  File "/opt/anaconda3/lib/python3.12/site-packages/pygformula/parametric_gformula/parametric_gformula.py", line 742, in fit
    boot_results_dicts[i]['boot_results'][j] for j in range(len(self.int_descript)))

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^
TypeError: 'NoneType' object is not subscriptable

 

Thanks for all the work on this package!

[LilJing]

Hi, thanks for posting this issue. The statsmodels used to fit a multinomial model in the pygformula does support categorical covariates with more than 6 levels. There is no limitation on the number of levels. However, the model may fail if some categories are rare. If that's the case, switching to fitting a one-vs-rest model with one-hot encoding is a reasonable alternative, as it is more robust to class imbalance. Could you check whether your categorical variable includes any categories with very few observations?

Best,
Jing

[jgonzhijon]

Hi Jing,

Thank you very much for your response and thank you for confirming the feasibility of adjusting a one-vs-rest.

Regarding having categories with low observations, it's something I considered. I merged the categories "Entertainment and Recreation Services," "Agriculture, Forestry, Fisheries and Mining," and "Information Services" just to test if it worked, and the error persisted.

Below, I attach the variable stratified by time point.

                                                               0    1    2    3    4    5    6    7    8    9
  Accommodation and Food Services                              1299 1212  908  687  483  357  247  158   99   52
  Agriculture, Forestry, Fisheries and Mining                   669  539  356  309  247  194  153  113   77   41
  Business and Repair Services                                 2131 1782 1424 1127  936  793  614  457  283  147
  Construction                                                 2081 1666 1199  897  720  578  425  304  183   73
  Entertainment and Recreation Services                         343  279  214  156  126   96   64   48   31   15
  Finance, Insurance, and Real Estate                          1264 1126  910  736  594  495  379  275  195  114
  Health Care and Social Assistance                            1709 1948 1599 1349 1058  835  609  459  285  138
  Information Services                                          323  502  382  321  257  203  154  118   77   41
  Manufacturing                                                4213 3339 2496 2039 1684 1362 1030  760  500  280
  Professional, Scientific, Technical and Educational Services 2406 1209  969  837  669  541  424  319  203  119
  Public Administration                                        1371 1187  970  867  696  570  456  366  255  144
  Transportation, Communications, and Other Public Utilities   2093 1611 1252  983  820  654  505  383  264  141
  Wholesale and Retail Trade                                   3848 2897 2236 1694 1310  994  740  509  353  189

All the best,
Juan

[LilJing]

Hi Juan,

Thanks for providing your data. Looking at the values for the categorical variable, I didn’t identify a clear issue that might be causing the error. To further debug this issue, could you try running your analysis using the R version of the gformula (https://github.com/CausalInference/gfoRmula) to see if a similar error occurs? Alternatively, if it’s convenient, feel free to send the full dataset to my email - I’ll do a cross-check to identify the potential issues.

Best,
Jing

[jgonzhijon]

Hi Jing,

Sorry for not being able to contact you earlier. Unfortunately, I can't share the original dataset. However, I've been running some tests with gfoRmula, and this issue doesn't appear—the program runs smoothly.
But the bootstrap doesn't seem to work correctly when I include the categorical variable. Below is the gfoRmula code I used (the variable names have been anonymized).

library(data.table)
library(gfoRmula)
library(nnet)
set.seed(1234)

df <- fread("./path/to/input.csv")
df[, `:=`(
  bs1 = as.integer(gender == "Female"),
  bs2 = as.integer(race_white == "White"),
  tv3 = as.factor(industry)
)]
df[, id := ID][, c("ID", "gender", "race_white", "industry") := NULL]
setorder(df, id, register_num)
df[, time := seq_len(.N) - 1L, by = id]

never_if_employed <- function(df, pool, intvar, intvals, time_name, t) {
  idx <- df[[time_name]] == t & df$tv1 == 1
  df[idx, (intvar) := 0]
  invisible(NULL)
}
always_if_employed <- function(df, pool, intvar, intvals, time_name, t) {
  idx <- df[[time_name]] == t & df$tv1 == 1
  df[idx, (intvar) := 1]
  invisible(NULL)
}

covnames <- c("tv1", "tv2", "tv3")
covtypes <- c("binary", "binary", "categorical")
covmodels <- list(
  tv1 ~ lag1_tv1 + tv2 + lag1_tv2 + tv3 + lag1_tv3,
  tv2 ~ lag1_tv2 + tv1 + lag1_tv1 + tv3 + lag1_tv3,
  tv3 ~ lag1_tv3 + tv1 + lag1_tv1 + tv2 + lag1_tv2
)
covparams <- list(covmodels = covmodels)

histories <- list(lagged)
histvars  <- list(covnames)

ymodel <- outcome ~
  tv2 + lag1_tv2 +
  tv1 + lag1_tv1 +
  tv3 + lag1_tv3 +
  bs1 + bs2

time_points  <- 10L
nsimul       <- 24047L
nsamples     <- 10L
int_descript <- c("Never if employed", "Always if employed")

fit <- gformula(
  obs_data      = df,
  id            = "id",
  time_name     = "time",
  time_points   = time_points,

  covnames      = covnames,
  covtypes      = covtypes,
  covparams     = covparams,
  histories     = histories,
  histvars      = histvars,

  basecovs      = c("bs1", "bs2"),
  ymodel        = ymodel,
  outcome_name  = "outcome",
  outcome_type  = "continuous_eof",

  seed        = 1234,
  nsimul      = nsimul,
  nsamples    = nsamples,
  boot        = TRUE,
  boot_diag   = TRUE,
  return_sims = TRUE,

  intervention1.tv2 = list(never_if_employed),
  intervention2.tv2 = list(always_if_employed),

  int_descript = int_descript,

  parallel   = TRUE,
  ncores     = 10,
  model_fits = TRUE
)

print(fit)

Here are my results. As you can see, the bootstrap standard errors are not estimated. However, this issue does not occur when a categorical variable is not included.

PREDICTED RISK UNDER MULTIPLE INTERVENTIONS

Intervention 	 Description
0 		 Natural course
1		 Never if employed
2		 Always if employed

Sample size = 24047, Monte Carlo sample size = 24047
Number of bootstrap samples = 10
Reference intervention = natural course (0)
 

     k Interv.  NP mean g-form mean Mean SE Mean lower 95% CI
 <num>   <num>    <num>       <num>   <num>             <num>
     9       0 2.788385    2.815759      NA                NA
     9       1       NA    2.821295      NA                NA
     9       2       NA    2.838735      NA                NA
 Mean upper 95% CI Mean ratio MR SE MR lower 95% CI
             <num>      <num> <num>           <num>
                NA   1.000000    NA              NA
                NA   1.001966    NA              NA
                NA   1.008160    NA              NA
 MR upper 95% CI Mean difference MD SE MD lower 95% CI
           <num>           <num> <num>           <num>
              NA     0.000000000    NA              NA
              NA     0.005535982    NA              NA
              NA     0.022976063    NA              NA
 MD upper 95% CI % Intervened On Aver % Intervened On
           <num>           <num>                <num>
              NA         0.00000             0.000000
              NA        40.32520             5.110409
              NA        96.72308            28.313303

I also don't know how to retrieve the simulated datasets in R to check if everything is working correctly — I just get an empty object.

Sorry for such a long message and for the delay. And thank you so much for everything!!

All the best,
Juan

[LilJing]

Hi Juan,

Thanks for posting! The simulated dataset can be retrieved by setting the argument "sim_data_b" to TRUE. Note that this only applies when bootstrapping is not used. The NA error might occur if the model fails to converge when fitting the bootstrap samples. Could you check whether any coefficients have NA estimates in the bootstrap replicates? You can set the argument "boot_diag" to TRUE and retrieve the bootstrap coefficients via "fit$bootcoeffs".

Best,
Jing