Why this simple mixed model fail to converge?
3
$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
r mixed-model lme4-nlme power
edited 10 hours ago
Robert Long
9,74622549
asked 12 hours ago
Philippe MassicottePhilippe Massicotte
284
$endgroup$
|
show 1 more comment
3
$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
r mixed-model lme4-nlme power
edited 10 hours ago
Robert Long
9,74622549
asked 12 hours ago
Philippe MassicottePhilippe Massicotte
284
$endgroup$
$begingroup$
I cannot replicate this problem, it works without errors for me. Did you update your packages?
$endgroup$
– user2974951
12 hours ago
$begingroup$
Seems I have everything updated. I will paste my session info.
$endgroup$
– Philippe Massicotte
12 hours ago
$begingroup$
I think it is from lmer in lmerTest. lmer from lme4 gives me the same results, but without the warning.
$endgroup$
– Philippe Massicotte
12 hours ago
1
$begingroup$
You have six data points in total, or is this a reduced example for posting?
$endgroup$
– Wayne
11 hours ago
1
$begingroup$
This is real data. So maybe there is a better analysis to use.
$endgroup$
– Philippe Massicotte
11 hours ago
|
show 1 more comment
3
3
3
1
$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
r mixed-model lme4-nlme power
edited 10 hours ago
Robert Long
9,74622549
asked 12 hours ago
Philippe MassicottePhilippe Massicotte
284
$endgroup$
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
r mixed-model lme4-nlme power
r mixed-model lme4-nlme power
edited 10 hours ago
Robert Long
9,74622549
asked 12 hours ago
Philippe MassicottePhilippe Massicotte
284
edited 10 hours ago
Robert Long
9,74622549
asked 12 hours ago
Philippe MassicottePhilippe Massicotte
284
edited 10 hours ago
Robert Long
9,74622549
edited 10 hours ago
Robert Long
9,74622549
edited 10 hours ago
Robert Long
9,74622549
9,74622549
asked 12 hours ago
Philippe MassicottePhilippe Massicotte
284
asked 12 hours ago
Philippe MassicottePhilippe Massicotte
284
asked 12 hours ago
Philippe MassicottePhilippe Massicotte
284
284
$begingroup$
I cannot replicate this problem, it works without errors for me. Did you update your packages?
$endgroup$
– user2974951
12 hours ago
$begingroup$
Seems I have everything updated. I will paste my session info.
$endgroup$
– Philippe Massicotte
12 hours ago
$begingroup$
I think it is from lmer in lmerTest. lmer from lme4 gives me the same results, but without the warning.
$endgroup$
– Philippe Massicotte
12 hours ago
1
$begingroup$
You have six data points in total, or is this a reduced example for posting?
$endgroup$
– Wayne
11 hours ago
1
$begingroup$
This is real data. So maybe there is a better analysis to use.
$endgroup$
– Philippe Massicotte
11 hours ago
|
show 1 more comment
$begingroup$
I cannot replicate this problem, it works without errors for me. Did you update your packages?
$endgroup$
– user2974951
12 hours ago
$begingroup$
Seems I have everything updated. I will paste my session info.
$endgroup$
– Philippe Massicotte
12 hours ago
$begingroup$
I think it is from lmer in lmerTest. lmer from lme4 gives me the same results, but without the warning.
$endgroup$
– Philippe Massicotte
12 hours ago
1
$begingroup$
You have six data points in total, or is this a reduced example for posting?
$endgroup$
– Wayne
11 hours ago
1
$begingroup$
This is real data. So maybe there is a better analysis to use.
$endgroup$
– Philippe Massicotte
11 hours ago
$begingroup$
I cannot replicate this problem, it works without errors for me. Did you update your packages?
$endgroup$
– user2974951
12 hours ago
$begingroup$
I cannot replicate this problem, it works without errors for me. Did you update your packages?
$endgroup$
– user2974951
12 hours ago
$begingroup$
Seems I have everything updated. I will paste my session info.
$endgroup$
– Philippe Massicotte
12 hours ago
$begingroup$
Seems I have everything updated. I will paste my session info.
$endgroup$
– Philippe Massicotte
12 hours ago
$begingroup$
I think it is from lmer in lmerTest. lmer from lme4 gives me the same results, but without the warning.
$endgroup$
– Philippe Massicotte
12 hours ago
$begingroup$
I think it is from lmer in lmerTest. lmer from lme4 gives me the same results, but without the warning.
$endgroup$
– Philippe Massicotte
12 hours ago
1
1
$begingroup$
You have six data points in total, or is this a reduced example for posting?
$endgroup$
– Wayne
11 hours ago
$begingroup$
You have six data points in total, or is this a reduced example for posting?
$endgroup$
– Wayne
11 hours ago
1
1
$begingroup$
This is real data. So maybe there is a better analysis to use.
$endgroup$
– Philippe Massicotte
11 hours ago
$begingroup$
This is real data. So maybe there is a better analysis to use.
$endgroup$
– Philippe Massicotte
11 hours ago
|
show 1 more comment
2 Answers
2
active
oldest
votes
6
$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 error. 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.
answered 10 hours ago
Robert LongRobert Long
9,74622549
$endgroup$
$begingroup$
Thank you Robert for this good explanation. Usingcultureas 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
10 hours ago
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
10 hours ago
$begingroup$
Do you think it's OK to fit a mixed model for n=6?
$endgroup$
– amoeba
10 hours ago
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
10 hours ago
add a comment |
8
$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.
answered 10 hours ago
WayneWayne
16.2k23976
$endgroup$
$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
9 hours ago
4
$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
9 hours ago
2
$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
9 hours ago
$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
3 hours ago
add a comment |
Your Answer
StackExchange.ifUsing("editor", function () {
return StackExchange.using("mathjaxEditing", function () {
StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix) {
StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
});
});
}, "mathjax-editing");
StackExchange.ready(function() {
var channelOptions = {
tags: "".split(" "),
id: "65"
};
initTagRenderer("".split(" "), "".split(" "), channelOptions);
StackExchange.using("externalEditor", function() {
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled) {
StackExchange.using("snippets", function() {
createEditor();
});
}
else {
createEditor();
}
});
function createEditor() {
StackExchange.prepareEditor({
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: false,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: null,
bindNavPrevention: true,
postfix: "",
imageUploader: {
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
},
onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
});
}
});
draft saved
draft discarded
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstats.stackexchange.com%2fquestions%2f391072%2fwhy-this-simple-mixed-model-fail-to-converge%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
2 Answers
2
active
oldest
votes
2 Answers
2
active
oldest
votes
active
oldest
votes
active
oldest
votes
6
$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 error. 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.
answered 10 hours ago
Robert LongRobert Long
9,74622549
$endgroup$
$begingroup$
Thank you Robert for this good explanation. Usingcultureas 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
10 hours ago
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
10 hours ago
$begingroup$
Do you think it's OK to fit a mixed model for n=6?
$endgroup$
– amoeba
10 hours ago
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
10 hours ago
add a comment |
6
$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 error. 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.
answered 10 hours ago
Robert LongRobert Long
9,74622549
$endgroup$
$begingroup$
Thank you Robert for this good explanation. Usingcultureas 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
10 hours ago
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
10 hours ago
$begingroup$
Do you think it's OK to fit a mixed model for n=6?
$endgroup$
– amoeba
10 hours ago
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
10 hours ago
add a comment |
6
6
6
$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 error. 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.
answered 10 hours ago
Robert LongRobert Long
9,74622549
$endgroup$
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 error. 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.
answered 10 hours ago
Robert LongRobert Long
9,74622549
edited 8 hours ago
answered 10 hours ago
Robert LongRobert Long
9,74622549
answered 10 hours ago
Robert LongRobert Long
9,74622549
answered 10 hours ago
Robert LongRobert Long
9,74622549
9,74622549
$begingroup$
Thank you Robert for this good explanation. Usingcultureas 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
10 hours ago
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
10 hours ago
$begingroup$
Do you think it's OK to fit a mixed model for n=6?
$endgroup$
– amoeba
10 hours ago
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
10 hours ago
add a comment |
$begingroup$
Thank you Robert for this good explanation. Usingcultureas 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
10 hours ago
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
10 hours ago
$begingroup$
Do you think it's OK to fit a mixed model for n=6?
$endgroup$
– amoeba
10 hours ago
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
10 hours ago
$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
10 hours ago
$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
10 hours ago
2
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
10 hours ago
$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
10 hours ago
$begingroup$
Do you think it's OK to fit a mixed model for n=6?
$endgroup$
– amoeba
10 hours ago
$begingroup$
Do you think it's OK to fit a mixed model for n=6?
$endgroup$
– amoeba
10 hours ago
2
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
10 hours ago
$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
10 hours ago
add a comment |
8
$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.
answered 10 hours ago
WayneWayne
16.2k23976
$endgroup$
$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
9 hours ago
4
$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
9 hours ago
2
$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
9 hours ago
$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
3 hours ago
add a comment |
8
$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.
answered 10 hours ago
WayneWayne
16.2k23976
$endgroup$
$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
9 hours ago
4
$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
9 hours ago
2
$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
9 hours ago
$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
3 hours ago
add a comment |
8
8
8
$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.
answered 10 hours ago
WayneWayne
16.2k23976
$endgroup$
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.
answered 10 hours ago
WayneWayne
16.2k23976
answered 10 hours ago
WayneWayne
16.2k23976
answered 10 hours ago
WayneWayne
16.2k23976
answered 10 hours ago
WayneWayne
16.2k23976
16.2k23976
$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
9 hours ago
4
$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
9 hours ago
2
$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
9 hours ago
$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
3 hours ago
add a comment |
$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
9 hours ago
4
$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
9 hours ago
2
$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
9 hours ago
$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
3 hours ago
$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
9 hours ago
$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
9 hours ago
4
4
$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
9 hours ago
$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
9 hours ago
2
2
$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
9 hours ago
$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
9 hours ago
$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
3 hours ago
$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
3 hours ago
add a comment |
draft saved
draft discarded
Thanks for contributing an answer to Cross Validated!
- Please be sure to answer the question. Provide details and share your research!
But avoid …
- Asking for help, clarification, or responding to other answers.
- Making statements based on opinion; back them up with references or personal experience.
Use MathJax to format equations. MathJax reference.
To learn more, see our tips on writing great answers.
draft saved
draft discarded
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstats.stackexchange.com%2fquestions%2f391072%2fwhy-this-simple-mixed-model-fail-to-converge%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
$begingroup$
I cannot replicate this problem, it works without errors for me. Did you update your packages?
$endgroup$
– user2974951
12 hours ago
$begingroup$
Seems I have everything updated. I will paste my session info.
$endgroup$
– Philippe Massicotte
12 hours ago
$begingroup$
I think it is from lmer in lmerTest. lmer from lme4 gives me the same results, but without the warning.
$endgroup$
– Philippe Massicotte
12 hours ago
1
$begingroup$
You have six data points in total, or is this a reduced example for posting?
$endgroup$
– Wayne
11 hours ago
1
$begingroup$
This is real data. So maybe there is a better analysis to use.
$endgroup$
– Philippe Massicotte
11 hours ago