[R-sig-ME] Comparing Gaussian and bêta regression

D. Rizopoulos d@rizopoulo@ @ending from er@@mu@mc@nl
Fri Sep 21 16:10:46 CEST 2018

In the case of only random intercepts, there are formulas that you can use to obtain the marginalized coefficients. I.e.,

\beta^M = \beta^SS / sqrt(1 + \kappa * \sigma_b^2),

where \beta^SS are the subject-specific coefficients, \kappa is a constant that depends on the link function, and \sigma_b^2 the variance of the random intercepts.

For the case with more random effects, Heagerty and colleagues have worked on marginalized models, but these are rather complicated.

However, recently Hedeker et al. (2017) came up with a nice a simple idea to obtain the marginalized coefficients. This is implemented in the function marginal_coefs() of my GLMMadaptive package; check https://drizopoulos.github.io/GLMMadaptive/articles/Methods_MixMod.html#marginalized-coefficients

Best,
Dimitris

From: Ben Bolker <bbolker using gmail.com<mailto:bbolker using gmail.com>>
Date: Friday, 21 Sep 2018, 4:02 PM
To: Emmanuel Curis <emmanuel.curis using parisdescartes.fr<mailto:emmanuel.curis using parisdescartes.fr>>, D. Rizopoulos <d.rizopoulos using erasmusmc.nl<mailto:d.rizopoulos using erasmusmc.nl>>
Cc: r-sig-mixed-models using r-project.org <r-sig-mixed-models using r-project.org<mailto:r-sig-mixed-models using r-project.org>>
Subject: Re: [R-sig-ME] Comparing Gaussian and bêta regression

Does anyone know offhand if there's R code (ideally a package)
floating around that implements these marginalization calculations for
mixed model estimates, by delta method or quadrature or simulation or
... ?  (The mixed-model ecosystem is getting pretty big and messy ...)
Do SAS/Stata/whatever have straightforward ways to do this that we could
copy?

cheers
Ben Bolker

On 2018-09-21 09:20 AM, Emmanuel Curis wrote:
> Thanks. It's more clear with the example on Cross Validated, but I'll
> have to think about it carefully.
>
> As I understand it now, the idea is as follow [I do not try to deal
> with the scale-change exactly, it's just to illustrate has I
> understand it]:
>
>  - let assume the fixed effect coefficient for sex is one, for
>    instance, and is is the only covariate for simplicity, except the
>    random effect for clinicians
>
>     => « score is increased by one in males »
>
>  - let assume clinician 1 has random effect +0.1 and clinician 2 -0.1
>
>  - interpretation should be done assuming all other variables
>    constants, including random effects, here clinicians
>
>     => score is increased by one in males for the same clinician
>
>  - if now I compare between two clinicians, and assuming no
>    interaction between sex and clinician, the difference between
>    "score male, clinician 1" and "score female, clinician 2" is
>      (+1 + 0.1) - (0 - 0.1) = 1.2
>
>  - now the question is what happens when averaging on all clinicians
>    (on the random effect)
>
>  - in the « usual » linear model, the difference would be
>      E( 1 + U1 - (0 + U2) ) = 1 + E(U1 - U2) = 1
>
>    => the « conditionnal » effet is the same than when averaging other
>    all clinicians
>      (but comparing between two specific clinicians still will give a
>       coefficient slightly different than +1)
>
>  - in any model with a non-linear link function, in the original
>    scale, the difference would be (with f the reciprocal of the link
>    function, and assuming « small » random effects to use a second
>    order Taylor expansion)
>
>     E( f(1 + U1) - f(0 - U2) ) = E( f(1 + U1) ) - E( f(0 - U2) )
>
>       ~ f(1) + E(f'(1) * U1) + E(f"(1) * U1²/2) -
>         f(0) - E(f'(0) * U2) - E(f"(0) * U2²/2)
>
>       = f(1) - f(0) + f"(1) E(U1²)/2 - f"(0) E(U2²)/2
>         -----------
>            delta
>
> and delta is the fixed effect interpretation « conditionnal » on the
> random effect, which is now different than the difference in the
> population averaged over all clinicians.
>
> Is this correct? [I guess Taylor expansion is unneeded here, I just
> used it to see more clearly the "bias"]
>
> Best regards,
>
> On Fri, Sep 21, 2018 at 07:54:45AM +0000, D. Rizopoulos wrote:
> « The estimates of the coefficients are asymptotically unbiased, and indeed this stems from the properties of maximum likelihood. But what I'm talking about is the *interpretation* of these coefficients. The use of a nonlinear link function complicates things in that regard.
> «
> « You may check more on this in this discussion on Cross Validated: https://stats.stackexchange.com/questions/365907/interpretation-of-fixed-effects-from-mixed-effect-logistic-regression
> «
> « Best,
> « Dimitris
> «
> « -----Original Message-----
> « From: Emmanuel Curis <emmanuel.curis using parisdescartes.fr>
> « Sent: Friday, September 21, 2018 9:34 AM
> « To: D. Rizopoulos <d.rizopoulos using erasmusmc.nl>
> « Cc: Ben Bolker <bbolker using gmail.com>; r-sig-mixed-models using r-project.org
> « Subject: Re: [R-sig-ME] Comparing Gaussian and bêta regression
> «
> « Hi,
> «
> « Thanks for the information.  I must say that presently, I do not understand what this means exactly, so I'm really curious about it.
> «
> « I was naively thinking that fixed effects coefficients obtained were estimates of the « true » coefficients in the model, and that because of the maximum likelihood asymptotic properties, were asymptotically unbiased.  Does this mean it is false? Or that they estimate something else than these « true » coefficients? Or may be that even the notion of « true coefficient » is not so clear?
> «
> « Would it be possible to have a detailed, simple example, either directly or through a reference?
> «
> « Best regards,
> «
> « On Wed, Sep 19, 2018 at 04:08:28PM +0000, D. Rizopoulos wrote:
> « « An additional issue with the interpretation from the Beta mixed model would be that because of the nonlinear link function, the fixed effects coefficients will have an interpretation conditional on the random effects. This same issue you also for example have in the mixed effects logistic and Poisson regression models. Most often this is not the interpretation you want, i.e., a marginal/population interpretation.
> « «
> « « Best,
> « « Dimitris
> « «
> « «
> « « -----Original Message-----
> « « From: R-sig-mixed-models <r-sig-mixed-models-bounces using r-project.org> On Behalf Of Emmanuel Curis « Sent: Wednesday, September 19, 2018 5:02 PM « To: Ben Bolker <bbolker using gmail.com> « Cc: r-sig-mixed-models using r-project.org « Subject: Re: [R-sig-ME] Comparing Gaussian and bêta regression « « Thank you very much for the hints.
> « «
> « « The « not very satisfactory » is from a "theoretical" point of view:
> « « I'm not very comfortable with modeling with a Gaussian a value constrained between 0 and 10, with the extremes obtained not so rarely.  From a practical point of view, it does not seem to produce unexpected results.  Of course, there are some effects that are borderline significant, that also makes the question uprise: what is the part of true signal and basically inadequate model in these effects? Still finding them with a more sounded model would make them a little bit more "trustable"...
> « «
> « « For the ordinal outcome: I have wrongly selected my example value, it induced in error, sorry ; the step is 0.1 and not 0.5 ; in practice, « 46 different values were observed. Integer values and, to a less extent, half-integer values are clearly over-represented, I guess because of inconscient rounding during scoring.  I don't know how to handle this in a model, however, but that's another problem, and may be there is no need for that.  But, for the ordinal aspect, I fear that would make too much parameters in the model...
> « «
> « « Just thinking... Would it be imaginable to make inferences on the beta-distribution model, since it seems to much better describe the data, but use the linear model on the raw scale just to have point-estimates of the changes in an easiest-to-interpret way?
> « « [despite it is problematic close to the boundaries...] « « Is the Gelmann & Hill book you're thinking about this one: ?
> « «
> « « Data Analysis Using Regression and Multilevel/Hierarchical Models  Cambridge University Press « ISBN-10: 052168689X « « On Wed, Sep 19, 2018 at 09:49:51AM -0400, Ben Bolker wrote:
> « « «
> « « «
> « « « On 2018-09-19 03:30 AM, Emmanuel Curis wrote:
> « « « > Hello,
> « « « >
> « « « > I'm doing my first try on bêta regression, with mixed effects model, « > and was wondering if my reasonning is correct...
> « « « >
> « « « > The context is a clinical study where the outcome is a score variable, « > with continuous values between 0 and 10 (both included) and, in « > practice, values with only one decimal digit (eg. 1.5) There is « > about 400 patients. Random effect is the clinician who does the « > examination and afterthat collects the score that evaluates its « > intervention.
> « « « >
> « « « > As a quick-and-dirty analysis, I did a linear mixed effect model on « > the raw data, with lmer. Residuals and random effects are not so bad, « > and results consistent & easy to interpret, but assuming a Gaussian « > distribution is not very satisfactory.
> « « «
> « « « Can you expand on why "not very satisfactory"?  Do you get unrealistic « predictions etc.?
> « « «
> « « «   This sounds like it could also be treated as an ordinal response (with
> « « « 21 values {0, 0.5, 1, ... 9.5, 10}).
> « « « >
> « « « > Hence, I tried a bêta regression on the data after the transformation « > (y/10 * (n-1) + 0.5) / n, and used glmmTMB for that. And of course I « > wondered if the fit was better.
> « « « >
> « « « > 1) Is it right that ln-likelihood of the model on the raw data
> « « « >    (Gaussian) and on the transformed data (bêta) cannot be compared,
> « « « >    because they involve probability densities and not probabilities,
> « « « >    hence depend on the data scale ?
> « « «
> « « «   You can compare log-likelihoods (actually technically they're
> « « « log-likelihood *densities*, which is where the problem comes from) if « you account for the scaling.  In this case since you're doing a linear « transformation the scaling should be pretty easy.
> « « « >
> « « « > 2) Is it right that the lmer model done on the raw data and the same
> « « « >    one done on the transformed data are conceptually the same, since
> « « « >    the transformation is linear — so that the ln-likelihood it gives
> « « « >    is « the same » expressed in the two different scales? (of course,
> « « « >    coefficients and so on will be different because of the scale
> « « « >    change)
> « « «
> « « «    Should be. (You could do a simple test of this ...)
> « « « >
> « « « > 3) And so, is it correct to compare the ln-likelihood (using logLik)
> « « « >    or the AIC given by glmmTMB with the bêta model and by lmer on
> « « « >    transformed data to compare the two models (raw data Gaussian vs
> « « « >    bêta)?
> « « «
> « « «   I would think so.
> « « « >
> « « « >    If so, the bêta model seems better than the Gaussian one. But now
> « « « >    comes the interpretation problem, other than « are coefficients
> « « « >    significantly different from 0? ».
> « « « >
> « « « > 4) Since the default link is the logit for the mean, interpretation is
> « « « >    not quite clear for me.  For the Gaussian model on raw data,
> « « « >    interpretation is clear, for instance « men score 1 point lower
> « « « >    than women in average ».  But how can the coefficients of the
> « « « >    bêta-model be back-converted in a similar fashion ?
> « « «
> « « «    You probably need to go read stuff about interpretation of
> « « « logit/log-odds  parameters: Gelman and Hill's book is good.
> « « «
> « « « Quick rules of thumb:
> « « «
> « « « * for β∆x small, as for log (proportional) « * for intermediate values, linear change in probability with « slope ≈ β/4 « * for large values, as for log ( 1 − x ) « >
> « « « >    Would it be easier to use a log link and expression changes in the
> « « « >    scale as percent changes on the mean?
> « « «
> « « «   This will work fine for low score values, but will run into trouble at
> « « « the upper end of the score range.
> « « «
> « « « >
> « « « > Thanks in advance,
> « « « >
> « « «
> « « « _______________________________________________
> « « « R-sig-mixed-models using r-project.org mailing list « https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> « «
> « « --
> « «                                 Emmanuel CURIS
> « «                                 emmanuel.curis using parisdescartes.fr
> « «
> « « Page WWW: http://emmanuel.curis.online.fr/index.html
> « «
> « « _______________________________________________
> « « R-sig-mixed-models using r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> «
> « --
> «                                 Emmanuel CURIS
> «                                 emmanuel.curis using parisdescartes.fr
> «
> « Page WWW: http://emmanuel.curis.online.fr/index.html
>

[[alternative HTML version deleted]]