1
$\begingroup$

I have produced a mixed-effects model (lmer) that is analysing the difference in area between three protocols of growing mini-brains. I have takes into account various fixed effects including the type of Media (food) they have, the type of Cell Line they are from and the Timepoint, as well as random effects such as Individual, as each individual's Area is tracked over time in repeated imaging. I include Timepoint as a random slope and/or as an additional random effect if the model is able to generate without errors, as this is the idealistic model. Autocorrelation is not much of an issue with the random effect structure included. I could not include the full dataset due to the character limit, but I've included a header for reference.

    CellLine Media Individual Timepoint       Area
 1:      CL1    EB         IA         1 0.08902973
 2:      CL1    EB         IA         3 0.13397213
 3:      CL1    NE         IA         5 0.22353000
 4:      CL1    NE         IA         7 0.31261003
 5:      CL1    NE         IA         9 0.37988532
 6:      CL1    ND         IA        11 0.54168553
 7:      CL1    ND         IA        13 0.64331523
 8:      CL1    ND         IA        15 0.91611337
 9:      CL1    ND         IA        17 1.23737488
10:      CL1    ND         IA        19 1.55369448
11:      CL1    ND         IA        21 1.85185604
12:      CL2    EB         CA         1 0.06675972
13:      CL2    EB         CA         3 0.11718475
14:      CL2    NE         CA         5 0.19933605
15:      CL2    NE         CA         7 0.26540415
16:      CL2    NE         CA         9 0.33455309
17:      CL2    ND         CA        11 0.45397616
18:      CL2    ND         CA        13 0.76065087
19:      CL2    ND         CA        15 1.04349630
20:      CL2    ND         CA        17 1.33352196
21:      CL2    ND         CA        19 1.56763996
22:      CL2    ND         CA        21 1.84198481
23:      CL3    EB          D         1 0.07304713
24:      CL3    EB          D         3 0.11285901
25:      CL3    NE          D         5 0.19936120
26:      CL3    NE          D         7 0.23331321
27:      CL3    NE          D         9 0.24856647
28:      CL3    ND          D        11 0.31204416
29:      CL3    ND          D        13 0.38930386
30:      CL3    ND          D        15 0.54777174
31:      CL3    ND          D        17 0.79907952
32:      CL3    ND          D        19 1.14717570
33:      CL3    ND          D        21 1.59278960

Untransformed model:

library(lmer)
Model1 <- lmer(Area) ~ CellLine*Timepoint*Media + (1|Individual) + (1|Timepoint), REML = TRUE, data = dat)

The problem is with the output of estimates and confidence intervals that I require for reporting the model. If I transform my response variable (Area) via sqrt, after I back transform the fixed effect, not all of my fixed and random estimates fit within their respective confidence intervals, i.e. fixed effects estimate: CellLine = 9.423180e-06, CIs = 8.146938e-04, 0.0012028990, or random effects variance: Individual = 0.002341655, CIs = 1.562912e-04, 0.0003964329. This only seems to impact on this particular dataset, others appear fine for both fixed and random effects.

Sqrt Model:

sqrtModel <- lmer(sqrt(Area) ~ CellLine*Timepoint*Media + (1|Individual) + (1|Timepoint), REML = TRUE, data = dat)

CIs <- "^"(confint(sqrtModel, level = 0.95, method ="profile", oldNames=FALSE),2)
CIs
FixedEstimates <- "^"(sqrtModel@beta,2)
FixedEstimates
RandomVariance <- as.data.frame(VarCorr(sqrtModel))
RandomVariance

head(CIs)
                                2.5 %       97.5 %
sd_(Intercept)|Individual 1.562912e-04 0.0003964329
sd_(Intercept)|Timepoint  2.744643e-05 0.0002202416
sigma                     9.905215e-04 0.0012176109
(Intercept)               4.071123e-02 0.0730789521
CellLineCL2               2.047092e-03 0.0002154916
CellLineCL3               8.146938e-04 0.0012028990

head(FixedEstimates)
5.572006e-02 2.335659e-04 9.423180e-06 1.531326e-03 1.617956e-01 8.937379e-03

head(RandomVariance)
          grp        var1 var2        vcov      sdcor
1 Individual (Intercept) <NA> 0.002341655 0.04839065
2  Timepoint (Intercept) <NA> 0.001568108 0.03959934
3   Residual        <NA> <NA> 0.005674150 0.07532695

However, I have more than one of the same growth analysis (7 in fact) and those that have more extreme separation such as the one I've attached lend towards considerably better diagnostics if the response variable is transformed by log. When I go to use exactly the same throughput to generate my estimates/CIs I find the random estimates are significantly outside of CIs, no matter whether I back transform them or not; this occurs in all log transformed models from different datasets.

Log Model:

LogModel <- lmer((log(Area)) ~ CellLine*Timepoint*Media + (1|Individual) + (1|Timepoint), REML = TRUE, data = dat)

CIs <- exp(confint(LogModel, level = 0.95, method ="profile", oldNames=FALSE))
CIs
FixedEstimates <- exp(logModel@beta)
FixedEstimates
RandomVariance <- as.data.frame(VarCorr(logModel))
RandomVariance

head(CIs)
                              2.5 %    97.5 %
sd_(Intercept)|Individual 1.03965780 1.0600446
sd_(Intercept)|Timepoint  1.01721667 1.0452862
sigma                     1.07369807 1.0820281
(Intercept)               0.05331121 0.0649370
CellLineCL2               0.82065574 0.9437256
CellLineCL3               0.93693142 1.0857551

head(FixedEstimates)
0.05883771 0.88004128 1.00860125 1.28488963 1.41834589 2.16451530

head(RandomVariance)
         grp        var1 var2        vcov      sdcor
1 Individual (Intercept) <NA> 0.002341655 0.04839065
2  Timepoint (Intercept) <NA> 0.001568108 0.03959934
3   Residual        <NA> <NA> 0.005674150 0.07532695

For reference, to support model validity I plot the raw data against the model's predictions, followed by an ACF plot and finally tables dictating estimate/CI outputs for fixed and random effects. I really think it would be useful in this circumstances to have the residual effects variance so I can show the relevance of including random effects in reducing residual deviance and accommodating for temporal autocorrelation.

Any help would be greatly appreciated.

$\endgroup$
4
  • $\begingroup$ Not directly related to the problem you ask about, but is there some reason why you are including Timepoint as both a fixed and a random effect? That doesn't seem to make much sense and it might account for the NA values in some of your outputs. Choose one (I'd suggest changing Timepoint to a continuous fixed-effect predictor and modeling it flexibly with a regression spline, removing the random effect) and see if your problem still exists. $\endgroup$
    – EdM
    Commented Jul 19, 2022 at 14:23
  • $\begingroup$ For such time-course data, generalized least squares, also called a "growth curve model," might be a good alternate choice instead of a mixed model. See Chapter 7 of Frank Harrell's course notes or book for an explanation and worked-through examples. $\endgroup$
    – EdM
    Commented Jul 19, 2022 at 14:29
  • $\begingroup$ @EdM it was upon suggestion from a textbook to include anything that could affect autocorrelation to be included in the random effects, but I did think it odd. I'm not familiar with a regression spline, but have removed the Timepoint as a random, and used bs(Timepoint, knots = c(5,11)), with the knots representing change in Media, and it improves fitted.predicted fit but the fixed-effects matrix becomes rank deficient and the autocorrelation worsens. Not sure if I'm using it correctly? $\endgroup$ Commented Jul 19, 2022 at 15:42
  • $\begingroup$ Generalized least squares can directly handle the autocorrelation over time, with multiple choices of the form. See Chapter 7 of the Harrell references. Chapters 2 of his references explain regression splines; they are implemented in the rcs() function of his rms package. With about 10 timepoints you might try 3 or 4 knots with rcs(). You seem to have systematic changes of culture media over time (EB to NE to ND), so continuous Timepoint and Media might be linearly dependent and thus giving you the rank deficiency. $\endgroup$
    – EdM
    Commented Jul 19, 2022 at 16:36

1 Answer 1

2
$\begingroup$

There are a couple of important points here, only one of which specifically has to do with mixed models.

  1. Back-transforming coefficients that have been estimated on transformed data is not as easy as applying the inverse transformation to the coefficients, except in the case of log transformation.
  2. random-effects parameters don't have the same form or interpretation as fixed-effect parameters, so the process of back-transforming them would be different in any case.

Illustrating point 1 with a sqrt transform on a linear model: suppose $\sqrt{y} = \beta_0 + \beta_1 x$. If we define the "coefficient associated with $x$" as "the expected change in $y$ for a 1-unit change in $x$", then we should consider $(\beta_0 + \beta_1 (x+1))^2 - (\beta_0 + \beta_1 x)^2) = \beta_1^2(2x + 1)+2 \beta_0 \beta_1 = \beta_1(\beta_1 + 2x +2\beta_0)$. Just squaring $\beta_1$ isn't going to give us the right answer!

On the other hand, it's very common that people "back-transform" the coefficients of log-linear models. What's going on here? Suppose $\log(y) = \beta_0 + \beta_1 x$. Then the proportional change in $y$, going from $x$ to $x+1$, is $\exp(\beta_0 + \beta_1(x+1))/\exp(\beta_0 + \beta_1 x) = \exp(\beta_1)$. This is where the standard "exponentiate your coefficients from a log-transformed fit" procedure comes from.1

Point 2 is a little deeper. Random effects coefficients don't have units of "$\Delta y/\Delta x$", they have units (in the case of the scalar/intercept-level REs of your model) of "variation in $y$". If you did fit a model to a transformed variable (say $f(y)$) and wanted to back-transform RE standard deviations from the scale of $f(y)$ to the scale of $y$, you would actually need to multiply them by $df^{-1}(y)/dy$ or divide them by $df(y)/dy$. In the case where $f(y) = \log(y)$, this is equivalent to multiplying the standard deviations by some reference value of $y$ ($1/d(\log(y))/dy = 1/(1/y) = y$).

Finally, if we look at the RE standard deviation estimates all on the same scale, the estimates are (properly) contained within the 95% CIs. Since exponentiation is a monotonic function, this will also be true if we exponentiate everything [although as explained above, doing so wouldn't make much sense ...]. (In your example above, you compared the exponentiated CIs to the non-exponentiated standard deviations ...)

sd_ci <- confint(LOGModel,
                 parm = "theta_",
                 method ="profile",
                 oldNames = FALSE)
sd_ests <- as.data.frame(VarCorr(LOGModel))[,"sdcor"]
setNames(data.frame(sd_ests, sd_ci), c("est", "lwr", "upr"))
                                 est        lwr        upr
sd_(Intercept)|Individual 0.04839065 0.03889163 0.05831130
sd_(Intercept)|Timepoint  0.03959934 0.01707015 0.04429069
sigma                     0.07532695 0.07110883 0.07883718

1 Log-odds-transformed responses are often exponentiated to get answers in terms of odds ratios, but that's another can of worms.

$\endgroup$
4
  • $\begingroup$ Thank you for putting this together, it made an awful lot of sense. I've applied the refined code to fixed and random effects for all my models and they seem to work great now so thank you so much for that! I have a few follow-up questions just to double check my understanding if you don't mind. I understand your formulae for Point 1, but is that to say I should report any estimates/CIs derived as they are outputted, i.e. untransformed, for those derived from non-logged models? $\endgroup$ Commented Jul 19, 2022 at 17:25
  • $\begingroup$ Also, with your code I have not had any issue with fixed or random effects of either sqrt(y) or log(y) models, what was wrong with the code I was using, as they look fairly similar to me (minus the exponentiating). $\endgroup$ Commented Jul 19, 2022 at 17:25
  • $\begingroup$ It's certainly easiest to report the coefficients/CIs in their 'raw' form. Exactly how you present the model results is going to depend a lot on what kind of conclusions you want to draw from the model and who your audience is ... (by the way, performance::check_model() could be a useful shortcut for a lot of the diagnostic plots in the code you sent me) $\endgroup$
    – Ben Bolker
    Commented Jul 19, 2022 at 17:49
  • $\begingroup$ Thank you for making me aware, you're quite right, brilliant little shortcut and provided me a few extra diagnostics, collinearity etc. You've helped wrap up a 3 year long headache, I'm very grateful :) $\endgroup$ Commented Jul 19, 2022 at 18:18

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.