1
$\begingroup$

Thanks for taking the time to read my question. I am using R to run a GLMM because I have taken repeated measurements on 100 individuals. The explanatory variable it is a proportion so I am using binomial family. The final model after a model selection is:

  1. treat: categorical variable with 4 levels
  2. place: categorical variable with 6 levels
  3. fdate: categorical variable with 13 levels
  4. fid: categorical variable with 100 levels

    M <- glmer(y ~ treatment + place + fdate + (1|fid), data = table, family = binomial(logit), glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5)))
    
    The single term deletions tell me that are all significant 
    Model:
    y ~ treatment + place + fdate6 + (1 | fid)
             Df    AIC    LRT   Pr(Chi)    
    <none>      828.25                     
    treatment 848.87 26.620 7.071e-06 ***
    place 830.27 12.014    0.0346 *  
    fdate6   12 851.13 46.881 4.887e-06 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    

Summary output

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: y ~ treatment + place + fdate + (1 | fid)
   Data: table
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))

     AIC      BIC   logLik deviance df.resid 
   828.3    942.0   -392.1    784.3     1278 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9511 -0.1789  0.0000  0.0000  5.0429 

Random effects:
 Groups Name        Variance Std.Dev.
 fid (Intercept) 0.1935   0.4399  
Number of obs: 1300, groups:  fgpsid, 100

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -1.17329    0.31499  -3.725 0.000195 ***
treatme02        -0.45551    0.26314  -1.731 0.083438 .  
treatme03        -1.63531    0.31361  -5.215 1.84e-07 ***
treatme04        -0.64832    0.30000  -2.161 0.030690 *  
place02          -0.57500    0.24866  -2.312 0.020757 *  
place03          -0.24367    0.38401  -0.635 0.525731    
place04          -0.48586    0.36544  -1.330 0.183675    
place05          -0.92812    0.41933  -2.213 0.026874 *  
place06           0.37485    0.49016   0.765 0.444412    
fdate62           1.23421    0.26453   4.666 3.08e-06 ***
fdate63           0.97002    0.27035   3.588 0.000333 ***
fdate64           1.12241    0.30254   3.710 0.000207 ***
fdate65           0.71227    0.43558   1.635 0.102005    
fdate66          -0.12241    0.66440  -0.184 0.853819    
fdate67           0.09972    0.80556   0.124 0.901484    
fdate68           2.19890    0.59646   3.687 0.000227 ***
fdate69           0.50455    0.56520   0.893 0.372024    
fdate610          0.92899    0.42959   2.162 0.030581 *  
fdate611          0.51314    0.56166   0.914 0.360920    
fdate612         -0.16656    0.43301  -0.385 0.700492    
fdate613          0.93091    0.31163   2.987 0.002815 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

My next question is regarding the levels of the fixed categorical variables. When I use the multcomp function my output is not significant. Is there something here that am I doing wrong? Which of these results should be reported in my manuscript? I do not understand why my variables are reported as significant in the summary, but not significant with I conduct pairwise comparisons.

> test02 <- glht(M, linfct = mcp(place= "Tukey"))
> summary(test02)
 Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: glmer(formula = y ~ treatment + place + fdate6 + (1 | fgpsid), 
    data = table, family = binomial(logit), control = glmerControl(optimizer = "bobyqa", 
        optCtrl = list(maxfun = 2e+05)))

Linear Hypotheses:
                             Estimate Std. Error z value Pr(>|z|)
place01 - place02 == 0      -0.57500    0.24866  -2.312    0.176
place03 - place02 == 0      -0.24367    0.38401  -0.635    0.987
place03 - place02 == 0      -0.48586    0.36544  -1.330    0.753
place04 - place02 == 0      -0.92812    0.41933  -2.213    0.217
place05 - place02 == 0       0.37485    0.49016   0.765    0.971
place06 - place01 == 0       0.33133    0.36982   0.896    0.943
place03 - place01 == 0       0.08914    0.32802   0.272    1.000
place04 - place01 == 0      -0.35313    0.40192  -0.879    0.947
place05 - place01 == 0       0.94985    0.43239   2.197    0.224
place03 - place06 == 0      -0.24219    0.41090  -0.589    0.991
place04 - place06 == 0      -0.68446    0.47250  -1.449    0.679
place05 - place06 == 0       0.61852    0.55486   1.115    0.865
place04 - place03 == 0      -0.44227    0.45707  -0.968    0.922
place05 - place03 == 0       0.86071    0.49542   1.737    0.486
place05 - place04 == 0       1.30298    0.57876   2.251    0.201
(Adjusted p values reported -- single-step method)
$\endgroup$
2
  • $\begingroup$ Can't comment on the multcomp thing, but you should generally report all your results (possibly in a supplementary material) - note what is significant but don't hide the non-significant results. Also it is not clear why are you doing model selection in the first place (inference? interpretation? decision making?), knowing your final goal might help people answer the question. $\endgroup$ Commented Aug 3, 2018 at 9:20
  • $\begingroup$ Thanks you for you answers @Martin Modrák. My goal is to know the effect of the treatment on the response variable and if there are differences between the treatment. During the data explotarion I saw differences on places on the response variable but I would like to know if there differences between them. The reason of model selection was to select the variable which could explain what happen on my system. $\endgroup$
    – R. Alcalá
    Commented Aug 3, 2018 at 9:39

1 Answer 1

2
$\begingroup$

Let's break this down into two parts:

1) What should you report: The model selection process you have run tells you that place should be important, but you didn't observe a significant contribution of individual places. That's IMHO totally OK and probably just means you have too little measurements per place. I wouldn't worry about it - it is important that you controlled for place as that might have an impact on your estimates of treatment - which are (I guess) of primary interest.

I would report the data you have shown as something like "All treatments have reduced the proportion of interest compared to treatment 01. The evidence is strongest for treatment 03 (<estimate, CI, p-value >), less so for treatment 04 (<estimate, CI, p-value >). For treatment 02 we cannot reject no effect at 0.05 level (<estimate, CI, p-value >)". Or you may report the contrasts between the treatments in a similar way.

You should generally report all your results (possibly in a supplementary material) - note what is significant but don't hide the non-significant results. Also make it clear you have run model selection.

2) Model selection When your goal is to make inferences about a treatment, model selection might not be a great idea - since you are choosing the "best" among multiple models, you are likely to get a model that is overly confident (really depends on the number of variables you've tried). It might be sensible to just use all your variables and combine that with some regularization or priors to avoid overfitting. You may want to give rstanarm a try - it has a drop-in replacement for glmer and will let you specify priors that will avoid overfitting even when using all the variables (all the information) you have i your data.

IMHO model selection has use only when simplifying the model is an important goal in itself (e.g. when you later want to make decisions with less measured variables) and even in that case approaches like the projpred package might be preferrable.

$\endgroup$

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.