1
$\begingroup$

I am trying to assess the effect of a treatment. $Y$ is the response variable with two levels - 0 and 1. The predictor variable (denoted with column name $grp$ in the data given below) also has two levels -- groups A and B (with and without treatment).

There are also some random effects that I have to consider. The samples were collected from three sources (denoted with column name $src$ in the data given below), two sources for group A and one for group B. Each source had a different environmental condition, which was quantified as one number (called $T$) for each source.

The data is as given below.

v = data.frame(src = c(rep("1", 773),
                      rep("2", 896),
                      rep("3", 2032)),
           T = c(rep(0.1452560, 773),
                      rep(0.1407996, 896),
                      rep(0.1730382, 2032)),
           Y = c(rep(1, 40), rep(0, 733),
                 rep(1, 50), rep(0, 846),
                 rep(1, 110), rep(0, 1922)),
           grp = c(rep("A", 773+896),
                   rep("B", 2032)))
v$Y = factor(v$Y)

I used lme4 package in R to fit generalized mixed model to the above data, once with the fixed effect ($grp$) and a second time without it.

library(lme4)

alt = glmer(Y ~ grp + (1|T) + (1|src), data = v, family = "binomial")

null = glmer(Y ~ (1|T) + (1|src), data = v, family = "binomial")

test = anova(null, alt)

Here is the summary of alt and null models.

> summary(alt)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [glmerMod]
Family: binomial  ( logit )
Formula: Y ~ grp + (1 | T) + (1 | src)
Data: v

 AIC      BIC   logLik deviance df.resid 
1564.2   1589.1   -778.1   1556.2     3697 

Scaled residuals: 
   Min      1Q  Median      3Q     Max 
-0.2392 -0.2392 -0.2392 -0.2387  4.1886 

Random effects:
Groups Name        Variance Std.Dev.
 T      (Intercept) 4e-14    2e-07   
 src    (Intercept) 4e-14    2e-07   
Number of obs: 3701, groups:  T, 3; src, 3

Fixed effects:
         Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.864737   0.108343 -26.441   <2e-16 ***
grpB         0.004096   0.146090   0.028    0.978    
---  
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
     (Intr)
grpB -0.741

> summary(null)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [glmerMod]
Family: binomial  ( logit )
Formula: Y ~ (1 | T) + (1 | src)
Data: v

   AIC      BIC   logLik deviance df.resid 
1562.2   1580.9   -778.1   1556.2     3698 

Scaled residuals: 
Min     1Q     Median     3Q    Max 
-0.239 -0.239 -0.239 -0.239  4.184 

Random effects:
Groups Name        Variance Std.Dev.
T      (Intercept) 0        0       
src    (Intercept) 0        0       
Number of obs: 3701, groups:  T, 3; src, 3

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -2.8625     0.0727  -39.37   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Here are the test details.

Data: v
Models:
null: Y ~ (1 | T) + (1 | src)
alt: Y ~ grp + (1 | T) + (1 | src)
     Df    AIC    BIC logLik deviance Chisq Chi Df Pr(>Chisq)
null  3 1562.2 1580.9 -778.1   1556.2                        
alt   4 1564.2 1589.1 -778.1   1556.2 8e-04      1     0.9776

I have following questions:

  1. Is it ok to use random effects here? I am concerned because each random effect variable has only three levels. I am wondering if there is a requirement for random effects to be continuous variables or have large number of levels.

  2. Because the variables $T$ and $src$ are correlated, I am wondering if it is ok to drop one of the two. Alternatively, is it better to use $T$ as fixed effect and $src$ as random effect?

  3. How do I find the effect size associated with the treatment? Also, I am unable to understand how to interpret negative intercept values in the summary for fixed effects. On a related note, why does the summary of null model also include fixed effects?

  4. From the test details should I concern with anything other than the p-value? Is there a command to get it from the test variable directly?

If there are any flaws in my analysis, please do point out. Thanks for looking into this!

$\endgroup$
  • $\begingroup$ Can you tell us a bit more about src and T? Are the 3 levels of src the only ones you are interested in or just a set supposed to be representative of a larger set of sources? Is T a number (e.g., temperature) describing environmental conditions for each source? $\endgroup$ – Isabella Ghement Nov 24 '18 at 4:37
  • $\begingroup$ @IsabellaGhement Yes, there are only 3 levels. You could think of src as three days on which the samples were collected and each day had a different temperature (T). $\endgroup$ – Mr K Nov 24 '18 at 5:11
  • $\begingroup$ It's hard to estimate the variance of a random intercept effect associated with the grouping variable src when src has only 3 levels. As a rule of thumb, you would want a minimum of 5 levels. Since T is a continuous variable, you shouldn't treat it as a grouping variable in your model - including T as a fixed effect in the model seems appropriate. $\endgroup$ – Isabella Ghement Nov 24 '18 at 5:41
2
$\begingroup$
  1. ok to use random effects here, but your results indicate that model without random effect would be good. There is no a requirement for random effects to be continuous variables or have large number of levels.

    Groups Name        Variance Std.Dev.
        T      (Intercept) 4e-14    2e-07   
        src    (Intercept) 4e-14    2e-07   
    
  2. You can try to use both $T$ and $src$ as fixed effect, but it seems they have no contribution to the model, so it is possible both of them will not be in teh model.

  3. About effect size, you need find which one you want. Intercept: exp(-2.864737)/(1+exp(-2.864737)) is the probability of response =1 in grpA. Null model still has intercept, and intercept is fixed effect.

  4. Your p-value from fixed effect in the first model is the same as that from test. No more test is needed.

    grpB         0.004096   0.146090   0.028    0.978    
    
    alt   4 1564.2 1589.1 -778.1   1556.2 8e-04      1     0.9776
    
$\endgroup$
  • $\begingroup$ I just need to quantify the effect of treatment. Is it 0.004096 from the summary of my alt model? $\endgroup$ – Mr K Nov 24 '18 at 23:03
  • $\begingroup$ Yes, odds ratio of grpB vs grpA is exp(0.004095). But obviously, it is not significant with p = 0.978. $\endgroup$ – user158565 Nov 24 '18 at 23:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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