6
$\begingroup$

I have data on which I would like to compare the mean of the value variable between day 28 and day 83. Because the experience involve pseudo-replication (culture), I was thinking to use a mixed model with an random effect on culture I have two questions. The first one is about this dataset:

  library(lme4)
#> Loading required package: Matrix
  library(lmerTest)
#> 
#> Attaching package: 'lmerTest'
#> The following object is masked from 'package:lme4':
#> 
#>     lmer
#> The following object is masked from 'package:stats':
#> 
#>     step

df <- structure(list(
  day_true = c(28, 83, 28, 83, 28, 83), value = c(
    758453.333333333,
    575133.333333333, 684160, 656933.333333333, 816840, 734700
  ),
  culture = c(1L, 1L, 2L, 2L, 3L, 3L)
), row.names = c(NA, -6L), class = c("data.frame"))

If I fit it as follows, I do not have warnings.

lme4::lmer(value ~ factor(day_true) + (1 | culture), data = df)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: value ~ factor(day_true) + (1 | culture)
#>    Data: df
#> REML criterion at convergence: 102.7974
#> Random effects:
#>  Groups   Name        Std.Dev.
#>  culture  (Intercept) 47535   
#>  Residual             55990   
#> Number of obs: 6, groups:  culture, 3
#> Fixed Effects:
#>        (Intercept)  factor(day_true)83  
#>             753151              -97562

However, if I fit it with lmerTest I have a warning. Should I care about it?

lmerTest::lmer(value ~ factor(day_true) + (1 | culture), data = df)
#> Warning in as_lmerModLT(model, devfun): Model may not have converged with 1
#> eigenvalue close to zero: 2.6e-09
#> Linear mixed model fit by REML ['lmerModLmerTest']
#> Formula: value ~ factor(day_true) + (1 | culture)
#>    Data: df
#> REML criterion at convergence: 102.7974
#> Random effects:
#>  Groups   Name        Std.Dev.
#>  culture  (Intercept) 47535   
#>  Residual             55990   
#> Number of obs: 6, groups:  culture, 3
#> Fixed Effects:
#>        (Intercept)  factor(day_true)83  
#>             753151              -97562

According to the resulting model, there is no significant difference. But visually it seems to have one. I just want to make sure to understand the convergence issue.

boxplot(value ~ factor(day_true), data = df)

My second question is about this data, for which I have a singular fit message. I can’t find the reason. Is it because I have very few points (n = 3 per group)? Alternatively, is there a better analysis to use to compare the mean of repeated measures between these two groups?

df <- structure(list(day_true = c(0, 28, 0, 28, 0, 28), value = c(
  34.6732447526395,
  31.5635584852635, 34.2763264775584, 32.1719125021771, 35.0747566289866,
  31.7318622838194
), culture = c(1L, 1L, 2L, 2L, 3L, 3L)), row.names = c(
  NA,
  -6L
), class = c("data.frame"))

  lme4::lmer(value ~ factor(day_true) + (1 | culture), data = df)
#> singular fit
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: value ~ factor(day_true) + (1 | culture)
#>    Data: df
#> REML criterion at convergence: 5.3578
#> Random effects:
#>  Groups   Name        Std.Dev.
#>  culture  (Intercept) 0.0000  
#>  Residual             0.3592  
#> Number of obs: 6, groups:  culture, 3
#> Fixed Effects:
#>        (Intercept)  factor(day_true)28  
#>             34.675              -2.852  
#> convergence code 0; 1 optimizer warnings; 0 lme4 warnings

Created on 2019-02-06 by the reprex package (v0.2.1)

Session info

devtools::session_info()
#> ─ Session info ──────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 3.5.2 (2018-12-20)
#>  os       Linux Mint 19.1             
#>  system   x86_64, linux-gnu           
#>  ui       X11                         
#>  language en_CA:en                    
#>  collate  en_CA.UTF-8                 
#>  ctype    en_CA.UTF-8                 
#>  tz       America/Toronto             
#>  date     2019-02-06                  
#> 
#> ─ Packages ──────────────────────────────────────────────────────────────
#>  package     * version    date       lib source                          
#>  assertthat    0.2.0      2017-04-11 [1] CRAN (R 3.5.1)                  
#>  backports     1.1.3      2018-12-14 [1] CRAN (R 3.5.1)                  
#>  callr         3.1.1      2018-12-21 [1] CRAN (R 3.5.2)                  
#>  cli           1.0.1      2018-09-25 [1] CRAN (R 3.5.1)                  
#>  colorspace    1.4-0      2019-01-13 [1] CRAN (R 3.5.2)                  
#>  crayon        1.3.4      2017-09-16 [1] CRAN (R 3.5.1)                  
#>  curl          3.3        2019-01-10 [1] CRAN (R 3.5.2)                  
#>  desc          1.2.0      2018-05-01 [1] CRAN (R 3.5.1)                  
#>  devtools      2.0.1      2018-10-26 [1] CRAN (R 3.5.1)                  
#>  digest        0.6.18     2018-10-10 [1] CRAN (R 3.5.1)                  
#>  dplyr         0.8.0      2019-02-01 [1] Github (tidyverse/dplyr@e9902f1)
#>  evaluate      0.12       2018-10-09 [1] CRAN (R 3.5.1)                  
#>  fs            1.2.6      2018-08-23 [1] CRAN (R 3.5.1)                  
#>  ggplot2       3.1.0      2018-10-25 [1] CRAN (R 3.5.1)                  
#>  glue          1.3.0      2018-07-17 [1] CRAN (R 3.5.1)                  
#>  gtable        0.2.0      2016-02-26 [1] CRAN (R 3.5.1)                  
#>  highr         0.7        2018-06-09 [1] CRAN (R 3.5.1)                  
#>  htmltools     0.3.6      2017-04-28 [1] CRAN (R 3.5.1)                  
#>  httr          1.4.0      2018-12-11 [1] CRAN (R 3.5.1)                  
#>  knitr         1.21       2018-12-10 [1] CRAN (R 3.5.1)                  
#>  lattice       0.20-38    2018-11-04 [1] CRAN (R 3.5.1)                  
#>  lazyeval      0.2.1      2017-10-29 [1] CRAN (R 3.5.1)                  
#>  lme4        * 1.1-20     2019-02-04 [1] CRAN (R 3.5.2)                  
#>  lmerTest    * 3.0-1      2018-04-23 [1] CRAN (R 3.5.2)                  
#>  magrittr      1.5        2014-11-22 [1] CRAN (R 3.5.1)                  
#>  MASS          7.3-51.1   2018-11-01 [1] CRAN (R 3.5.1)                  
#>  Matrix      * 1.2-15     2018-11-01 [1] CRAN (R 3.5.1)                  
#>  memoise       1.1.0      2017-04-21 [1] CRAN (R 3.5.1)                  
#>  mime          0.6        2018-10-05 [1] CRAN (R 3.5.1)                  
#>  minqa         1.2.4      2014-10-09 [1] CRAN (R 3.5.1)                  
#>  munsell       0.5.0      2018-06-12 [1] CRAN (R 3.5.1)                  
#>  nlme          3.1-137    2018-04-07 [1] CRAN (R 3.5.2)                  
#>  nloptr        1.2.1      2018-10-03 [1] CRAN (R 3.5.1)                  
#>  numDeriv      2016.8-1   2016-08-27 [1] CRAN (R 3.5.1)                  
#>  pillar        1.3.1.9000 2019-02-01 [1] Github (r-lib/pillar@3a54b8d)   
#>  pkgbuild      1.0.2      2018-10-16 [1] CRAN (R 3.5.1)                  
#>  pkgconfig     2.0.2      2018-08-16 [1] CRAN (R 3.5.1)                  
#>  pkgload       1.0.2      2018-10-29 [1] CRAN (R 3.5.1)                  
#>  plyr          1.8.4      2016-06-08 [1] CRAN (R 3.5.1)                  
#>  prettyunits   1.0.2      2015-07-13 [1] CRAN (R 3.5.1)                  
#>  processx      3.2.1      2018-12-05 [1] CRAN (R 3.5.1)                  
#>  ps            1.3.0      2018-12-21 [1] CRAN (R 3.5.2)                  
#>  purrr         0.3.0      2019-01-27 [1] CRAN (R 3.5.2)                  
#>  R6            2.3.0      2018-10-04 [1] CRAN (R 3.5.1)                  
#>  Rcpp          1.0.0      2018-11-07 [1] CRAN (R 3.5.1)                  
#>  remotes       2.0.2      2018-10-30 [1] CRAN (R 3.5.1)                  
#>  rlang         0.3.1.9000 2019-02-01 [1] Github (tidyverse/rlang@7243c6d)
#>  rmarkdown     1.11       2018-12-08 [1] CRAN (R 3.5.1)                  
#>  rprojroot     1.3-2      2018-01-03 [1] CRAN (R 3.5.1)                  
#>  scales        1.0.0      2018-08-09 [1] CRAN (R 3.5.1)                  
#>  sessioninfo   1.1.1      2018-11-05 [1] CRAN (R 3.5.1)                  
#>  stringi       1.2.4      2018-07-20 [1] CRAN (R 3.5.1)                  
#>  stringr       1.3.1      2018-05-10 [1] CRAN (R 3.5.1)                  
#>  testthat      2.0.1      2018-10-13 [1] CRAN (R 3.5.1)                  
#>  tibble        2.0.1      2019-01-12 [1] CRAN (R 3.5.2)                  
#>  tidyselect    0.2.5      2018-10-11 [1] CRAN (R 3.5.1)                  
#>  usethis       1.4.0      2018-08-14 [1] CRAN (R 3.5.1)                  
#>  withr         2.1.2      2018-03-15 [1] CRAN (R 3.5.1)                  
#>  xfun          0.4        2018-10-23 [1] CRAN (R 3.5.1)                  
#>  xml2          1.2.0      2018-01-24 [1] CRAN (R 3.5.1)                  
#>  yaml          2.2.0      2018-07-25 [1] CRAN (R 3.5.1)                  
#> 
#> [1] /home/pmassicotte/R/x86_64-pc-linux-gnu-library/3.5
#> [2] /usr/local/lib/R/site-library
#> [3] /usr/lib/R/site-library
#> [4] /usr/lib/R/library

Improve this question
$\endgroup$
7
  • $\begingroup$ I cannot replicate this problem, it works without errors for me. Did you update your packages? $\endgroup$
    – user2974951
    Feb 6, 2019 at 13:21
  • $\begingroup$ Seems I have everything updated. I will paste my session info. $\endgroup$
    – Philippe Massicotte
    Feb 6, 2019 at 13:23
  • $\begingroup$ I think it is from lmer in lmerTest. lmer from lme4 gives me the same results, but without the warning. $\endgroup$
    – Philippe Massicotte
    Feb 6, 2019 at 13:28
  • 1
    $\begingroup$ You have six data points in total, or is this a reduced example for posting? $\endgroup$
    – Wayne
    Feb 6, 2019 at 14:01
  • 1
    $\begingroup$ This is real data. So maybe there is a better analysis to use. $\endgroup$
    – Philippe Massicotte
    Feb 6, 2019 at 14:07

2 Answers 2

Reset to default
14
$\begingroup$

This is, in all likelihood, not a warning that you need to worry about. As you can see, the parameter estimates are the same in both cases. The version of lmer in lmertest apparently has a more conservative check for convergence than the current lme4 version.

The problem in lmertest::lmer is caused by the variables being on vastly different scales, which can make some of the optimisation routines unstable, and thus generate the warning. If you re-scale, the problem vanishes:

> df$value <- df$value / 10000

> lmerTest::lmer(value ~ factor(day_true) + (1 | culture), data = df)

Linear mixed model fit by REML ['lmerModLmerTest']
Formula: value ~ factor(day_true) + (1 | culture)
   Data: df
REML criterion at convergence: 29.1147
Random effects:
 Groups   Name        Std.Dev.
 culture  (Intercept) 4.753   
 Residual             5.599   
Number of obs: 6, groups:  culture, 3
Fixed Effects:
       (Intercept)  factor(day_true)83  
            75.315              -9.756

Note that, since you have only 3 groups, it may be better to model culture as fixed:

> lm(value ~ factor(day_true) + culture, data = df)
Coefficients:
   (Intercept)  factor(day_true)83             culture  
        64.417              -9.756               5.449

Clearly the estimate for the fixed effect of day_true is the same in both analyses.

The reason for not finding a statistically significant estimate, this is because the sample size is so small. It is highly preferable to run a "power analysis" prior to collecting data and fitting the model.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thank you Robert for this good explanation. Using culture as fixed might be a solution. Based on the results, however, it seems there is no difference between day 28 and day 83. Is that correct to say that? If you see the plot I made above, I seems to have a difference. $\endgroup$
    – Philippe Massicotte
    Feb 6, 2019 at 14:51
  • 2
    $\begingroup$ @PhilippeMassicotte No, it isn't saying there is no difference. The difference is estimated as -9.756 in both models, so they are both saying the same thing. If you are referring to the estimate's p-value or confidence interval, then don't pay much attention to this - with such a small sample you it could be difficult to get a "significant" result and you won't have enough statistical power to detect the effect size (this should have been checked before collecting the data). $\endgroup$
    – Robert Long
    Feb 6, 2019 at 15:02
  • $\begingroup$ Do you think it's OK to fit a mixed model for n=6? $\endgroup$
    – amoeba
    Feb 6, 2019 at 15:17
  • 2
    $\begingroup$ @amoeba no I don't that's why I suggested the linear model with fixed effects for the grouping variable, in my answer. $\endgroup$
    – Robert Long
    Feb 6, 2019 at 15:22
11
$\begingroup$

I can’t speak to the calculation issue, but six data points is not very much at all, much less if you want to fit random effects — of which you only have 2 and only 3 examples of each.

The not-statistically-significant result makes a lot of sense here: you have a tiny amount of data that hints at something, but not enough to calculate anything well enough to draw any general conclusions.

Your box plots are misleading you: it’s calculating means, quartiles, etc with only 3 numbers for each box. There are more visualized features than data in this case.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thank you for your helpful comment. As I understand, given the little number of observations it is difficult to get clear conclusions. Do you see any other venues to tackle this? $\endgroup$
    – Philippe Massicotte
    Feb 6, 2019 at 15:42
  • 5
    $\begingroup$ @PhilippeMassicotte Get rid of the box plot and simply plot three dots for actual values and connect them by lines. In terms of test, it's a paired t-test, but there is little sense in running tests for n=3. Just state that the value decreased for all 3 cultures. $\endgroup$
    – amoeba
    Feb 6, 2019 at 15:45
  • 3
    $\begingroup$ @PhilippeMassicotte: What amoeba says. Basically you don't have enough data to generalize and make an overall statement about the world, but you can speak about the data you have and what happened in those cases. If all three lines go down, that begins to suggest that it's possible that there is an effect there. Which might justify doing a new experiment with enough data to determine if you've discovered something about the way things actually work. $\endgroup$
    – Wayne
    Feb 6, 2019 at 16:25
  • 2
    $\begingroup$ it would be nice for these comments to be preserved in an answer (either edits to @Wayne's answer, or as a separate answer). There can never be enough sensible explanations about what to do with tiny data sets. $\endgroup$
    – Ben Bolker
    Feb 6, 2019 at 21:46
  • $\begingroup$ I can add them, or @amoeba could add an answer as well, or we can leave them as comments. Either way. $\endgroup$
    – Wayne
    Feb 7, 2019 at 1:59

Your Answer

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

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