0
$\begingroup$

I am fitting a linear mixed model with lme4 in R. The model has a single factor (des_days) with 4 levels (-1,1,14,48), and I am using random intercept and slopes. The model fit looks like this:

> lCtr <- lmeControl(maxIter = 5000, niterEM = 500, msMaxIter=1000, msMaxEval=1000, msVerbose = FALSE, opt = 'optim')
> lev.lm <- lme(data ~ des_days, random = ~des_days|ratID, data=data_red_trf, na.action=na.omit, method = "ML", control=lCtr)
> summary(lev.lm)

Linear mixed-effects model fit by maximum likelihood
 Data: data_red_trf 
        AIC       BIC   logLik
  -3289.891 -3216.717 1659.946

Random effects:
 Formula: ~des_days | ratID
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev      Corr                
(Intercept) 0.018205256 (Intr) ds_dy1 ds_d14
des_days1   0.014267902 -0.116              
des_days14  0.012291117  0.289  0.916       
des_days48  0.009514837 -0.866 -0.393 -0.725
Residual    0.043158947                     

Fixed effects: data ~ des_days 
                 Value   Std.Error  DF   t-value p-value
(Intercept)  0.8274313 0.007937938 962 104.23757  0.0000
des_days1   -0.0026322 0.007443294 962  -0.35363  0.7237
des_days14  -0.0011319 0.006635512 962  -0.17058  0.8646
des_days48   0.0112579 0.005452614 962   2.06469  0.0392
 Correlation: 
           (Intr) ds_dy1 ds_d14
des_days1  -0.213              
des_days14  0.068  0.694       
des_days48 -0.738 -0.078 -0.208

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-4.692631026 -0.524900364  0.006026586  0.514647345  3.766144781 

Number of Observations: 971
Number of Groups: 6 

I can clearly use the previous results to compare the estimations of each "des_day" to the intercept, using the provided t-statistics. This would suggest a slight difference between des_days48 and intercept (p=0.0392).

However, if I use post-hoc tests (z-statistics), this result changes (suggesting that des_days48 is non-different than intercept, p=0.101):

> ph_conditional <- c("des_days1  = 0",
                      "des_days14  = 0",
                      "des_days48 = 0");
> lev.ph <- glht(lev.lm, linfct = ph_conditional);
> summary(lev.ph)

Simultaneous Tests for General Linear Hypotheses

Fit: lme.formula(fixed = data ~ des_days, data = data_red_trf, random = ~des_days | 
    ratID, method = "ML", na.action = na.omit, control = lCtr)

Linear Hypotheses:
                 Estimate Std. Error z value Pr(>|z|)
des_days1 == 0  -0.002632   0.007428  -0.354    0.971
des_days14 == 0 -0.001132   0.006622  -0.171    0.996
des_days48 == 0  0.011258   0.005441   2.069    0.101
(Adjusted p values reported -- single-step method)

What result should be trusted and why? Thanks a lot in advance!

EDIT: We figured out that the p-values of the coefficient estimates and those of the post-hoc tests differ because the latter are adjusted with Bonferroni correction. I wonder, than, which results one should use. Does anybody have an answer?

$\endgroup$
  • $\begingroup$ "Adjusted p values reported" means that they were adjusted for multiple comparisons. $\endgroup$ – amoeba Mar 12 '18 at 20:43
  • $\begingroup$ ah , you are right, I overlooked it! Thanks! Shall I then use the adjusted p-values for z-scores, or the non-adjusted ones for t-statistics? And why? $\endgroup$ – Cristiano Mar 12 '18 at 20:47

Your Answer

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

Browse other questions tagged or ask your own question.