# Generalised Logistic Regression with mixed effects - Binomial model creation / selection

I have a dataset with 72 individuals from 5 separate groups, with repeated measures - each individual was sampled 4 times. All data is binomial, and in most cases there are more 0s than 1s.

I have a lot of potential predictor variables BacteriaA,BacteriaB etc, (9+ depending on how I group the data), and a number of outcome variables I want to test DiseaseA, Disease B, etc.

data <- read.csv("https://pastebin.com/raw/gwEFqh79 ")
cols <-   colnames(data)
data[cols] <- lapply(data[cols], factor)
str(data)

'data.frame':   500 obs. of  18 variables:
$$Group : Factor w/ 5 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...$$ ID       : Factor w/ 73 levels "E1","E10","E11",..: 1 1 1 1 2 2 2 2 2 2 ...
$$Time : Factor w/ 4 levels "1","2","3","4": 2 4 2 4 1 2 3 1 2 3 ...$$ DiseaseA : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$$DiseaseB : Factor w/ 2 levels "0","1": 2 1 2 1 1 1 1 1 1 1 ...$$ DiseaseC : Factor w/ 2 levels "0","1": 1 2 1 2 1 1 1 1 1 1 ...
$$DiseaseD : Factor w/ 2 levels "0","1": 2 2 2 2 1 1 1 1 1 1 ...$$ DiseaseE : Factor w/ 2 levels "0","1": 2 2 2 2 1 1 1 1 1 1 ...
$$DiseaseF : Factor w/ 2 levels "0","1": 2 2 2 2 1 1 1 1 1 1 ...$$ BacteriaA: Factor w/ 2 levels "0","1": 2 1 1 1 1 2 1 2 1 1 ...
$$BacteriaB: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...$$ BacteriaC: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 2 1 1 ...
$$BacteriaD: Factor w/ 2 levels "0","1": 2 1 2 1 2 2 1 1 2 1 ...$$ BacteriaE: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$$BacteriaF: Factor w/ 2 levels "0","1": 1 1 1 2 1 1 1 1 1 1 ...$$ BacteriaG: Factor w/ 2 levels "0","1": 1 2 2 2 1 2 2 1 2 2 ...
$$BacteriaH: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...$$ BacteriaI: Factor w/ 2 levels "0","1": 2 1 1 1 1 1 1 1 1 1 ...



I have run a model using glmer because I believe I need a mixed effects model, and I have binomial data.

model1 <- glmer(DiseaseA ~ (1|Group) + (1|ID) + Time + BacteriaA + BacteriaB + BacteriaC + BacteriaD + BacteriaE + BacteriaF + BacteriaG + BacteriaH + BacteriaI, data = data, family = 'binomial'(link = "logit"))

summary(model1)

Model failed to converge with max|grad| = 0.0337228 (tol = 0.001, component 1)Model failed to converge with max|grad| = 0.0337228 (tol = 0.001, component 1)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial  ( logit )
Formula: DiseaseA ~ (1 | Group) + (1 | ID) + Time + BacteriaA + BacteriaB +
BacteriaC + BacteriaD + BacteriaE + BacteriaF + BacteriaG +      BacteriaH + BacteriaI
Data: data

AIC      BIC   logLik deviance df.resid
426.2    489.4   -198.1    396.2      485

Scaled residuals:
Min      1Q  Median      3Q     Max
-1.2452 -0.3965 -0.2829 -0.1610  4.8949

Random effects:
Groups Name        Variance  Std.Dev.
ID     (Intercept) 1.298e+00 1.139325
Group  (Intercept) 3.841e-05 0.006198
Number of obs: 500, groups:  ID, 73; Group, 5

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -1.5414     0.4423  -3.485 0.000493 ***
Time2         0.6009     0.4518   1.330 0.183535
Time3         0.4617     0.4782   0.966 0.334287
Time4         0.4136     0.4577   0.904 0.366223
BacteriaA1   -0.5733     0.4593  -1.248 0.211939
BacteriaB1   -0.6688     0.5057  -1.323 0.185963
BacteriaC1   -0.6867     0.4371  -1.571 0.116228
BacteriaD1   -0.7494     0.3458  -2.167 0.030236 *
BacteriaE1   -0.1796     0.7715  -0.233 0.815964
BacteriaF1   -1.3041     0.6412  -2.034 0.041972 *
BacteriaG1   -0.7576     0.3260  -2.324 0.020133 *
BacteriaH1  -10.1035   124.1785  -0.081 0.935154
BacteriaI1   -0.1340     0.7987  -0.168 0.866712
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation matrix not shown by default, as p = 13 > 12.
Use print(x, correlation=TRUE)  or
vcov(x)        if you need it

convergence code: 0
Model failed to converge with max|grad| = 0.0337228 (tol = 0.001, component 1)


I am super new to modelling, never mind mixed models, and I'm just really stuck on what I should do next.

Is this the right model? How do I work out the odds ratios from this? What does Model failed to converge with max|grad| mean?

And most of all, how do I choose a model from this? Is there a stepwise way of dealing with so many variables?

• It's worth pausing to see if the sample size is adequate. The simplest thing you can do with binary outcome data is to estimate the probability that Y=1. This takes a sample size of n=96 independent subjects, to get a margin of error of +/- 0.1 in estimating the unknown probability. You are trying to do something that will require a much larger sample size than that, although having repeated measures does boost the effective sample size a bit. Commented May 9, 2023 at 11:32

I think your questions are a little too broad, and it's likely you should consider other approaches before jumping into a generalized linear mixed-effects model (GLMM), especially if you're super new to modelling.

• There is no "right model" -- remember, "all models are wrong, but some are useful". To determine if yours is even useful, you need to perform model validation (look at residuals, homoskedasticity, etc.) and consider what standard you'd want to hold your model to.
• Because you're fitting a model with the logit link, your regression parameters are already transformed into log-odds space (NOTE: if you have many more zeroes than ones, you may wan to consider the cloglog link instead). Assume you were to fit a null model (no predictors, only an intercept term) through model2 <- glm(DiseaseA ~ 1, data = data, family = 'binomial'(link = "logit")), and let's say the intercept parameter was estimated to be -2. Your model would be written as follows:

$$\log(\frac{p}{1-p}) = \beta_0 = -2$$

This means the log odds are -2 -- you would exponentiate to get the odds ($$\frac{p}{1-p} = e^{\beta_0} = e^{-2} = 0.1353$$) and back-transform to get the probability ($$p = \frac{0.1353}{0.1353+1} = 0.1192$$). Your log-odds will follow this process, although you have many more regression parameters.
• See this post.
• To your final question -- again, it is a bit too broad to really answer. There are many ways to perform model selection, and stepwise selection as you suggest is typically frowned upon. If you are new to modelling, how did you know to fit this model in the first place? Mixed-effects GLMMs are not easy and you should not just jump into one without a lot of appreciation for the statistics behind them. I suggest you read Ben Bolker's FAQ on GLMMs as a helpful resource to start (here).
• I wouldn't choose a link function based on the distribution of the outcome. Commented May 9, 2023 at 11:31