library("brglm2")      
     data("anorexia", package = "MASS")
     anorexML <- glm(Postwt ~ Prewt + Treat + offset(Prewt), family = gaussian, data = anorexia)
     update(anorexML, method = "brglmFit", type = "correction", control = brglm_control(trace = 1))
#> Error in brglmControl(epsilon = 1e-06, maxit = 100, trace = 1, response_adjustment = NULL, : unused arguments (Trans = expression(dispersion), inverseTrans = expression(transformed_dispersion))

But

library("brglm2")      
     data("anorexia", package = "MASS")
     anorexML <- glm(Postwt ~ Prewt + Treat + offset(Prewt), family = gaussian, data = anorexia)
     update(anorexML, method = "brglmFit", type = "correction", trace = 1)
#> Coefficients update: Outer/Inner iteration:  001/001
#> Dispersion update:   Outer iteration:     001 
#> max |step|: 0.141944      max |gradient|: 0.002168 
#> Coefficients update: Outer/Inner iteration:  001/001
#> Dispersion update:   Outer iteration:     001 
#> max |step|: 0.012124      max |gradient|: 0.000107
#> 
#> Call:  glm(formula = Postwt ~ Prewt + Treat + offset(Prewt), family = gaussian, 
#>     data = anorexia, method = "brglmFit", type = "correction", 
#>     trace = 1)
#> 
#> Coefficients:
#> (Intercept)        Prewt    TreatCont      TreatFT  
#>     49.7711      -0.5655      -4.0971       4.5631  
#> 
#> Degrees of Freedom: 71 Total (i.e. Null);  68 Residual
#> Null Deviance:       4525 
#> Residual Deviance: 3311  AIC: 490

Hi. I am trying unsuccessfully to add a random effect to a model in brmultinom. I am not sure if the package does not support random effects, or if I am just coding it incorrectly. I am using the following command:

Multi15_250 <- brmultinom(Behavior ~ PrevBehav + (1 | RandWal), data = BfrAftr15_25, ref = "Swm", type = "AS_mixed")

I get the following error:

Error in contrasts<-(*tmp*, value = contr.funs[1 + isOF[nn]]) :
contrasts can be applied only to factors with 2 or more levels
In addition: Warning message:
In Ops.factor(1, RandWal) : ‘|’ not meaningful for factors

Thank you!

Rebecca

Hi. Could you please tell me what is meant by the "linear.predictor" slot in the brmultinom object and its summary. When I reconstruct the linear predictors by hand and take the anti-logit, I get the probabilities I used to construct the data. I cannot tell what the linear predictors in the brmultinom object are. thanks.

Rebecca

Hello Dr. kosmidis,

Is the algorithm of "brglm2 AS_mean" the same as "brglm"? (pseudo-data method)
I ask the question because when I used "brglm" the results always finite. However, I got warning message said "brglmFit did not converge" when I use "brglm2 AS_mean".
I hope to understand what happened. Thank you!

I am using brmultinom with a small data set (120 cases). Is there a function or command to find the confidence intervals with the brmultinom function? I tried using the confint function but it gives no real output or error message.

All predictors are binary categorical data, except for "Number of Births per Year". The dependent outcome variable "Financial Sustainability" has 3 categories (best, good, poor)
###Set working directory
setwd("C:/Users/Melanie/OneDrive")

###install and library brglm2 package for analysis
install.packages("brglm2")
library(brglm2)

###install and library readxl package to read in data
install.packages("readxl")
library(readxl)

###Read in the data, make the outcome a factor
birthcenter<-read_excel("birthcenterdata")
birthcenter$Financial Sustainability<-as.factor(birthcenter$Financial Sustainability)
bccomplete<-birthcenter[complete.cases(birthcenter),]

Model D: uses "1" as baseline to compare "good" to "poor" and "excellent" to "poor."

###Model D shows each estimated coefficient and associated standard errors
modelD <- brmultinom(Financial Sustainability ~ Number of Births per Year +
Liability Insurance+ Accreditation +Urban/Suburban +
For Profit + Over 3 years +
Licensure Availability,data=bccomplete,ref=1)

summarizes the model

summary(modelD)

Odds ratios

exp(coefficients(modelD))

confidence intervals

confint(modelD,level = 0.95)

confint(modelD,level = 0.95)
2.5 % 97.5 %

Thank you for any / all thoughts + suggestions.

Hello,

I fit a model with the glm function (setting method="brglmFit"). When I implemented the step function from the stats package, the trace showed that the model with the lowest AIC value is not chosen. I do not want to share the data but I can share the trace. Please let me know if there is more information you need.

This is to avoid the need for an extra QR decomposition of the model matrix prior to fitting

Based on Wald and score statistics

After fitting two models with brmultinom using type = "AS_mixed" to handle quasi-complete separation, is it appropriate for me to test for differences in the same manner in which I would construct a LRT?

For example, I fit models "BRM0" and "BRM", used the anova command and compared the results to a chi-squared distribution as follows:

anova(BRM0,BRM)

Resid. Df Resid. Dev Df Deviance

1 1748 592.95

2 1746 593.67 2 -0.72255

pchisq(-0.72255,df=2,lower.tail = FALSE)

Could you kindly know if this is appropriate, and if so, what to call it?

Thank you!!

Rebecca

Is there a chance you forgot to add the parameter penalty when calculating AIC? For example, in a model I fit called Multi15_25.1, that I summarized in the object SmryMulti15_25.1, when I call SmryMulti15_25.1$AIC , it returns 543.0035. When I call SmryMulti15_25.1$logLik, it returns 'log Lik.' -271.5017 (df=6). And -2*(-271.5017) = 543.0035.

Also, can you please tell me what is the number the package calls "aic"? It is different than what the package calls "AIC."

Thank you.

Rebecca

devtools::install_github("ikosmidis/brglm2")

Gives:

Error in namespaceExport(ns, exports) : undefined exports: adjustedz.glm
Error: loading failed
Execution halted
ERROR: loading failed for 'i386', 'x64'

Hello,

I was excited to find this package to solve my complete separation problem. My count data is overdispersed so need to use the negative binomial family, which I understand is run using brnb(). I keep getting this error: NA/NaN/Inf in foreign function call (arg 1), even though I am specifying all the necessary terms:

salmonella_fm <- production ~ Treat_comboSource_Salinity +Rarity+ Region+Treat_comboTurion_length_num+replanted.x

N_T_6_2020_nbn <- brglm2::brnb(salmonella_fm,
data=turion_harvest_2021_all_gen_initial_turion,link = "log",transformation = "inverse", type = "AS_mean")

I've noticed that the error disappears if I use type="ML" but all other types trigger the error.

Could you please advise?

Thanks!

-Carrie

I want to test the separation, so I used this code :

glm(MODEL_FORM_GLM, data=train_std, family = binomial("logit"),method = "detect_separation")

As you can see below this does not show the variables responsible for the separation.

Implementation: | Linear program: | Purpose:
Separation: TRUE
Warning message:
'detect_separation' will be removed from 'brglm2' at version 0.8. A new version of 'detect_separation' is now maintained in the 'detectseparation' package.

glm(Histology ~ Cigarette + med1 * med2 * med3 * med4 + Age + Gender + Stage + Race, family = binomial("logit"), data = data2, method = "detect_separation")

gives：

Error in separator(x = x, y = y, linear_program = control$linear_program, : unexpected result from lpSolveAPI for primal test

data2.zip

PZIET-brglm2_bug.txt

When using brglmFit() for a logistic regression with one coefficient, NA is returned. However using brglm() from the brglm package returns a real valued coefficient. The code snippet below illustrates this issue,

# Author: Patrick Zietkiewicz
# File purpose: Investigate possible bug in brglm2 package
library(brglm)
library(brglm2)
dat <- MASS::Pima.tr
X <- as.matrix(dat[,1:7])
Y <- as.matrix(dat[,'type'] == 'Yes', 1, 0)
# fit models
res_brglm <- brglm(Y ~ -1 + X[,'bmi'],
                 family = binomial(link='logit'),
                 pl = FALSE)
res_brglm2 <- brglmFit(x = X[,'bmi'], y = Y, 
        family = binomial(link='logit'),
        intercept = FALSE)
cat('brglm:', res_brglm$coefficients,'\n')
cat('brglm2:', res_brglm2$coefficients,'\n')
cat('done')

Please see the attached txt file for a copy of the code snippet above

Hi Dr. kosmidis,

I got two warning message when I fit a logistic regression with one covariate, no intercept.
1: brglmFit: algorithm did not converge
2: brglmFit: fitted probabilities numerically 0 or 1 occurred
I got no problem when fitting standard logistic regression but got this problem when I use type = "MPL_Jeffreys".

Just a dataset with 50 samples and it look like this

I hope to understand what might happen here. Thank you.

First of all, thank you for your great work! I have tried your package. It seems to work for glm with family = "binomial". For the family = "gaussian", almost the same coefficients are generated when comparing standard glm and a glm with method = "brglmFit", type="correction". The difference in the coefficients is less than 10^-13. I guess there is something wrong with the link function and the calculation of the coefficients?

I am attempting to run a multinomial regression that has 3 levels of the categorical DV (reference category = 1) with a correction for data separation. I continue to receive an error that package ‘brglmfit’ is not available (for R version 3.6.1) and have attempted to run the analysis with package brglm2 instead. However, now with the below pasted syntax, I am getting the following error: Error in eval(extras, data, env) : object 'Freq' not found.

When I remove the weights = Freq syntax, the analysis simply doesn't run and does not give an error message.

brmultinomModel <- brmultinom(OCDPTSD3 ~ PCLintrusion+PCLavoidance+PCLcognitions+PCLarousal+DOCScontam+DOCSresponsibility+
DOCSunacceptablethts+DOCSjustright+LECSexAssault+LECCausedHarm, weights = Freq,
data = OCD_PTSD_valid_n_595_1, type = "AS_mean", ref = 1)

OCD-PTSD valid n = 595_1.xlsx

Thank you in advance for your help!

Hi, thanks for your work on this package. I am by no means an expert but I tried to read up on penalized logistic regression and encountered several times the claim that in principle Firth regression has better properties when it comes to convergence compared to standard GLM with logit link. I was therefore surprised to see the exact reverse on my data. I was hoping you could help me understand why this happens.

I narrowed down the problem in my regressions to a single variable. I have a longitudinal dataset where I use the year of the observation as a categorical variable, so I get intercepts for each year.

The data is available here: https://github.com/michalovadek/misc/blob/main/data.RDS

glm(dv ~ year, data = data, family = binomial()) # converges
glm(dv ~ year, data = data, family = binomial(), method = "brglmFit") # does not converge

Any advice would be much appreciated

Does brglm2 currently support ordered logistic regression? I had no luck trying to pass a cumulative() family object to brglmFit.

Thank you for the great package!

Fin

I was just wondering if it would also be possible to support quasipoisson & quasibinomial? Or is there any currently supported method to take into account overdispersion in binomial GLMs?