Performing GLMM using binomial data


#1

Hi,
I am running an R code using glmer function

m38 <- glmer(smear ~ sex + age + hiv_status + treatmenthistory + month + (1|serialnumber) + (1|region), data = phd_dat, family = binomial("logit"), control = glmerControl(optimizer = "optimx", calc.derivs = FALSE, optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE))
             ,nAGQ = 1)

and get the error:

Error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,  : 
  (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate
In addition: Warning message:
In optwrap(optimizer, devfun, start, rho$lower, control = control,  :
  convergence code 9999 from optimx

Any help is highly appreciated.


#2

Can you please provide a minimal reprex (reproducible example)? The goal of a reprex is to make it as easy as possible for me to recreate your problem so that I can fix it: please help me help you!

If you've never heard of a reprex before, start by reading "What is a reprex", and follow the advice further down that page.


#3

Hi Max,

Many thanks for looking at the error code and the will to assist. Please see below a reprex of my problem ( Its my first time using reprex, hope this is exactly what youve asked for)

Looking forward to your response

m36 <- glmer(smear ~ sex + age + hiv_status + treatmenthistory + month + (1|serialnumber)
             , data = phd_dat, family = binomial("logit"), control = glmerControl(optimizer = "bobyqa",, optCtrl = list(maxfun = 100000)),
             nAGQ = 25)
#> Error in glmer(smear ~ sex + age + hiv_status + treatmenthistory + month + : could not find function "glmer"

Created on 2018-10-11 by the reprex package (v0.2.1)


#4

Hi,

the code you show uses the reprex package (good job! :+1:), but it's not really a reprex in the sense of a reproducible example. The goal of a reprex is to allow us to reproduce, on our machines, the same error (or an error as close as possible), to the error you're getting:

Error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,  : 
  (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate
In addition: Warning message:
In optwrap(optimizer, devfun, start, rho$lower, control = control,  :
  convergence code 9999 from optimx

Thus, your expression, even if generated with the reprex package, is not a reprex, because it generates a different error:

#> Error in glmer(smear ~ sex + age + hiv_status + treatmenthistory + month + : could not find function "glmer"

I'll show you quickly how to modify your reprex stub into an actual reprex, and then give you a link to learn more about building reprexes yourself: first of all, in a reprex you need to load all and only the packages which are strictly necessary to reproduce your problem. In your case, the only package needed is lme4, thus start by modifying your code to:

library(lme4)
m36 <- glmer(smear ~ sex + age + hiv_status + treatmenthistory + month + (1|serialnumber),
             data = phd_dat, family = binomial("logit"), 
             control = glmerControl(optimizer = "bobyqa",, optCtrl = list(maxfun = 100000)),
             nAGQ = 25)

Now, try creating a reprex with it:

library(lme4)
#> Loading required package: Matrix
m36 <- glmer(smear ~ sex + age + hiv_status + treatmenthistory + month + (1|serialnumber),
             data = phd_dat, family = binomial("logit"), 
             control = glmerControl(optimizer = "bobyqa",, optCtrl = list(maxfun = 100000)),
             nAGQ = 25)
#> Error: 'data' not found, and some variables missing from formula environment

Created on 2018-10-11 by the reprex package (v0.2.1)

Session info
sessionInfo()
#> R version 3.5.0 (2018-04-23)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 16.04.5 LTS
#> 
#> Matrix products: default
#> BLAS: /usr/lib/atlas-base/atlas/libblas.so.3.0
#> LAPACK: /usr/lib/atlas-base/atlas/liblapack.so.3.0
#> 
#> locale:
#>  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
#>  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
#>  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
#> [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] lme4_1.1-18-1 Matrix_1.2-14
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_0.12.19    lattice_0.20-35 digest_0.6.18   rprojroot_1.3-2
#>  [5] MASS_7.3-49     grid_3.5.0      nlme_3.1-137    backports_1.1.2
#>  [9] magrittr_1.5    evaluate_0.12   stringi_1.2.4   minqa_1.2.4    
#> [13] nloptr_1.2.1    rmarkdown_1.10  splines_3.5.0   tools_3.5.0    
#> [17] stringr_1.3.1   yaml_2.2.0      compiler_3.5.0  htmltools_0.3.6
#> [21] knitr_1.20

See? You're getting a new error, and again an error which is not the same you're getting. Why? Because now your code includes all and only the necessary packages, but it doesn't include the data necessary to make it run (which I guess reside in the phd_dat dataframe). Now, probably you don't want to share publicly data containing information on smear, sex, age, hiv_status, treatmenthistory, etc. for the subjects you're following. Or maybe you want, but it's a big data set and thus it would be inconvenient for you to share it with us (and for us to test it). In this case, you have three options:

  • reproduce your error with one of the default datasets included in base R or in the package lme4. Thisi s not always possible.
  • try creating a small, fake data set, which reproduces the error you're getting. This is a great option.
  • extract a random subset from your actual data set, and include it into your code using the function dput. This works if you don't have sensitive data.

For more info about reprexes, please see here .

PS your issue seems to be optimization-related. Once we manage to have a reprex up and running, it could be a good idea for you to ask on Stack Overflow or Cross Validated...I don't want to jinx you, but frequentist estimation of GLMMs can be...not straightforward. Maybe you'd have better luck with a Bayesian approach:


#5

You might also find the the GLMM FAQ helpful! Here’s what it has to say about your error message:

  • PIRLS step-halvings failed to reduce deviance in pwrssUpdate
    • When using lme4 to fit GLMMs with link functions that do not automatically constrain the response to the allowable range of the distributional family (e.g. binomial models with a log link, where the estimated probability can be >1, or inverse-Gamma models, where the estimated mean can be negative), it is not unusual to get this error. This occurs because lme4 doesn’t do anything to constrain the predicted values, so NaN values pop up, which aren’t handled gracefully. If possible, switch to a link function to one that constrains the response (e.g. logit link for binomial or log link for Gamma).

#6

Andrea,

Thanks for taking time to look at my nightmare, moreso, giving an elaborate guide.

My efforts are giving this


library(lme4)
#> Loading required package: Matrix

test <- read_dta("test_dat.dta")
#> Error in read_dta("test_dat.dta"): could not find function "read_dta"

dput(test)
#> Error in dput(test): object 'test' not found

#dput(test_dat)
testm <- glmer(smear ~ sex + age + hiv_status + treatmenthistory + typeoftb + month + (1|serialnumber), data = test, family = binomial("logit"), control = glmerControl(optimizer = "bobyqa",optCtrl = list(maxfun = 100000)),  nAGQ = 1)
#> Error: 'data' not found, and some variables missing from formula environment
summary(testm)
#> Error in summary(testm): object 'testm' not found

Created on 2018-10-13 by the reprex package (v0.2.1)

Kindly advice,

Thanks


#7

Hi Jcblum,

The code already has specified the link function as "logit" but the error is still there

Kindly assist

Thanks


#8

Dear Richard,

until you provide us with a working reprex, I'm afraid I can't help you. As you can see, your reprex doesn't reproduce the error you're having:

#> Error in read_dta("test_dat.dta"): could not find function "read_dta"

I need a reprex which gives the same error you're getting:

Error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,  : 
  (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate
In addition: Warning message:
In optwrap(optimizer, devfun, start, rho$lower, control = control,  :
  convergence code 9999 from optimx

Absent that, I can only suggest you to post on Stack Overflow or Cross Validated, but I'm not sure you'll get the support you need there, without a working reprex.

Best,

Andrea