2
$\begingroup$

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?

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$ Commented May 9, 2023 at 11:32

1 Answer 1

0
$\begingroup$

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.

To briefly answer your questions:

  • 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).
$\endgroup$
1
  • $\begingroup$ I wouldn't choose a link function based on the distribution of the outcome. $\endgroup$ Commented May 9, 2023 at 11:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.