# 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?

• 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.