[R-sig-ME] [EXT] Re: Choice of distribution for random effects

[J.D. Haltigan]

Correct. The idea was to replicate the lmer *with the addition of a random
effects distribution* which is not possible in lmer (to see what numeric
differences in estimate effects might arise).

So, the fixed portion of the model, if I understand lmer correctly, should
be the same here. Whether I specify Beta or Gamma for the random effects
(in hglm2) yields the same error. I am wondering what might be the reason
for this.

At the risk of asking an elementary question here: When I add the random
effects specification for this model, I understand it as accounting for the
B/W cluster variance (union). In effect, by adding this random effects
component, I am saying that unions are *not* interchangeable and are drawn
from a universe of random draws of unions. I am borrowing a bit from the
parlance of Generalizability Theory, so apologies in advance if my language
is a bit cryptic.

JD

On Thu, Jul 21, 2022 at 6:11 PM Andrew Robinson <apro using unimelb.edu.au> wrote:

> As noted earlier, I’m not familiar with hglm, but a cursory glance
> suggests that those are not the same model.  You appear to have specified
> Beta distributed random effects in the call to hglm2.
>
> Cheers,
>
> Andrew
>
> --
> Andrew Robinson
> Chief Executive Officer, CEBRA and Professor of Biosecurity,
> School/s of BioSciences and Mathematics & Statistics
> University of Melbourne, VIC 3010 Australia
> Tel: (+61) 0403 138 955
> Email: apro using unimelb.edu.au
> Website: https://researchers.ms.unimelb.edu.au/~apro@unimelb/
>
> I acknowledge the Traditional Owners of the land I inhabit, and pay my
> respects to their Elders.
> On 22 Jul 2022, 7:45 AM +1000, J.D. Haltigan <jhaltiga using gmail.com>, wrote:
>
> *External email:* Please exercise caution
> ------------------------------
> Thanks for this continued deep dive as it is very instructive for me. I am
> currently also reading the Wood GAM book which is very helpful.
>
> As an aside, and perhaps relevant to the link discourse, in the current
> project I am working on, I fit an lmer with:
>
> ***lme4_1_B =
> lmer(posXsymp~treatment+proper_mask_base+prop_resp_ill_base_2 + pairID + (1
> | union), data = bdata.raw1)#lme4 package***
> which converges and provides results that are sensible (compared to models
> that have been run using fixed effects only).
>
> However, when I try to use HGLM to fit this model, it does not want to
> cooperate and I am wondering why. If I specify the HGLM as should be the
> case for linear models with family = Guassian, as:
>
> ***HGLM2_1_A = hglm2(posXsymp~
> treatment+proper_mask_base+prop_resp_ill_base_2 + pairID + (1 | union),
> family=gaussian(link=identity),
>                   rand.family=Beta(link=logit), maxit = 100, data =
> bdata.raw1)#HGLM package using hglm2***
>
> I get: Error in glm.fit(x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
> 1,  :
>   NA/NaN/Inf in 'x'
> In addition: Warning message:
> In hglm.default(X = X, y = Y, Z = Z, family = family, rand.family =
> rand.family,  :
>   Residuals numerically 0 are replaced by 1e-8
> >
>
> I am wondering if this is b/c "pairID" is a factor variable. I simply
> don't understand why lmer would fit this but HGLM returns that error.
>
> On Thu, Jul 21, 2022 at 4:18 PM Andrew Robinson <apro using unimelb.edu.au>
> wrote:
>
>> And, with my apologies for being fastidious, as you say the family
>> typically specifies the link function *in the software*, selecting the
>> (so-called) canonical link.
>>
>> That doesn’t make it the right link function for any given modeling
>> exercise. See, e.g., https://aip.scitation.org/doi/pdf/10.1063/1.5139815,
>> and particularly Figure 1.
>>
>> As Joe Hilbe and I wrote, it is not clear that there is any statistical
>> benefit in the canonization of a particular link function for each family.
>>
>> Cheers,
>>
>> Andrew
>> On 22 Jul 2022, 2:04 AM +1000, Ben Bolker <bbolker using gmail.com>, wrote:
>> I want to clarify one tiny thing.
>>
>> The family typically specifies the *link function* as well as the
>> conditional distribution, although this is often done by default (e.g.
>> poisson -> log, binomial -> logit/log-odds). The link function
>> specifies the scale on which the data are fitted, so although the random
>> effects are always Gaussian *on the scale of the linear predictor*
>> (i.e., the transformed scale), they vary on the scale of the data. If
>> the distribution of group-level intercepts is Gaussian on the log scale,
>> then it's log-Normal on the data (count) scale.
>>
>>
>> On 2022-07-20 11:36 p.m., Andrew Robinson wrote:
>> That is correct.
>>
>> Cheers,
>>
>> Andrew
>>
>> --
>> Andrew Robinson
>> Chief Executive Officer, CEBRA and Professor of Biosecurity,
>> School/s of BioSciences and Mathematics & Statistics
>> University of Melbourne, VIC 3010 Australia
>> Tel: (+61) 0403 138 955
>> Email: apro using unimelb.edu.au
>> Website: https://researchers.ms.unimelb.edu.au/~apro@unimelb/
>>
>> I acknowledge the Traditional Owners of the land I inhabit, and pay my
>> respects to their Elders.
>> On 21 Jul 2022, 1:04 PM +1000, J.D. Haltigan <jhaltiga using gmail.com>, wrote:
>> This cross-validated post seemed helpful and aligned with everything
>> previously mentioned here, and I wanted to share.
>>
>> https://stats.stackexchange.com/questions/190763/how-to-decide-which-glm-family-to-use
>>
>> Thanks for all of your insights as this is quite a deep dive for me.
>>
>> To be sure: in glmer, when one specifies the GLM family, that then applies
>> to both fixed & random effects, correct? In other words, there is no way
>> to
>> separately specify a distribution for the random effects in glmer? I know
>> I
>> asked this question to Ben on the lme4 Git, but I wanted to confirm again.
>>
>> JD
>>
>> On Wed, Jul 20, 2022 at 8:35 PM Andrew Robinson <apro using unimelb.edu.au>
>> wrote:
>>
>> Ben is correct that the traditional alignment of certain link functions
>> with certain members of the exponential family (the so-called 'canonical'
>> links) are merely for analytical convenience and have no statistical
>> justification.
>>
>> Canonical link functions exercise an unhealthy influence on statistical
>> practice, IMO. OTOH, there are circumstances in which they may afford
>> parameter estimates that are easier to interpret.
>>
>> I'm not familiar with hglm but I suspect you may need to use the log-link
>> function for the random effects in order to get something like the glmer
>> output.
>>
>> Cheers,
>>
>> Andrew
>>
>> --
>> Andrew Robinson
>> Chief Executive Officer, CEBRA and Professor of Biosecurity,
>> School/s of BioSciences and Mathematics & Statistics
>> University of Melbourne, VIC 3010 Australia
>> Tel: (+61) 0403 138 955
>> Email: apro using unimelb.edu.au
>> Website: https://researchers.ms.unimelb.edu.au/~apro@unimelb/
>>
>> I acknowledge the Traditional Owners of the land I inhabit, and pay my
>> respects to their Elders.
>> On 21 Jul 2022, 10:24 AM +1000, J.D. Haltigan <jhaltiga using gmail.com>,
>> wrote:
>>
>> Thanks, Ben. For example, here is a comparison of glmer & HGLM using
>> Guassian random effects for HGLM as suggested by Andrew. The magnitude of
>> the point estimates are about the same, but the nominal significance for
>> 'treatment' is only significant when considering the Gaussian random
>> effects from lmer.
>>
>> NB: In the project I am working on, the point is to show sensitivity of
>> results to different parameterizations of models. I understand the point
>> estimate magnitudes are about the same, but that is less relevant to this
>> exercise.
>>
>> *lme4_5_B = glmer(posXsymp~
>> treatment+proper_mask_base+prop_resp_ill_base_2
>> + pairID + (1 | union), family = "poisson", nAGQ=0, data =
>> bdata.raw3)#lme4
>> package using glmer*
>>
>> summary(lme4_5_B)
>>
>> Generalized linear mixed model fit by maximum likelihood (Adaptive
>> Gauss-Hermite Quadrature, nAGQ = 0) ['glmerMod']
>> Family: poisson ( log )
>> Formula: posXsymp ~ treatment + proper_mask_base + prop_resp_ill_base_2 +
>> pairID + (1 | union)
>> Data: bdata.raw3
>>
>> AIC BIC logLik deviance df.resid
>> 25054.5 27939.7 -12254.2 24508.5 287075
>>
>> Scaled residuals:
>> Min 1Q Median 3Q Max
>> -0.214597 -0.100237 -0.078562 -0.054499 40.646525
>>
>> Random effects:
>> Groups Name Variance Std.Dev.
>> union (Intercept) 0.007269442 0.08526102
>> Number of obs: 287348, groups: union, 538
>>
>> Fixed effects:
>> Estimate Std. Error z value
>> Pr(>|z|)
>> (Intercept) -5.987724526 0.581280589 -10.30092 <
>> 0.000000000000000222 ***
>> treatment -0.094657198 0.044546951 -2.12489
>> 0.03359612 *
>>
>>
>>
>> ***********************************************************************************************************************************************************************************************************************************
>>
>> *HGLM2_5_Test = hglm2(posXsymp~
>> treatment+proper_mask_base+prop_resp_ill_base_2 + pairID + (1 | union),
>> family =poisson (link = log), rand.family =gaussian(link=identity),
>> data = bdata.raw3)#HGLM package using hglm2*
>>
>> summary(HGLM2_5_Test)
>>
>> Call:
>> hglm2.formula(meanmodel = posXsymp ~ treatment + proper_mask_base +
>> prop_resp_ill_base_2 + pairID + (1 | union), data = bdata.raw3,
>> family = poisson(link = log), rand.family = gaussian(link = identity))
>>
>> ----------
>> MEAN MODEL
>> ----------
>>
>> Summary of the fixed effects estimates:
>>
>> Estimate Std. Error t-value Pr(>|t|)
>> (Intercept) -6.02638 1.09414 -5.508 0.0000000364 ***
>> treatment -0.02977 0.13624 -0.218 0.8271
>>
>> On Wed, Jul 20, 2022 at 10:42 AM Ben Bolker <bbolker using gmail.com> wrote:
>>
>> I believe that there are some mathematically natural pairings (i.e. a
>> Gamma random effect + a Poisson response, a Beta-distributed random
>> effect + a binomial response), but I don't know if there's much
>> theoretical justification other than analytical convenience.
>>
>> (If you have a categorical response then the natural random effect
>> would be Dirichlet, i.e. a Dirichlet-multinomial marginal distribution,
>> but do you really want to go down that rabbit hole?)
>>
>> I strongly second Andrew's recommendation to check that the
>> difference is not a difference between packages/estimation procedures.
>>
>> On 2022-07-20 5:26 a.m., Andrew Robinson wrote:
>>
>> You should really also verify that the hglm is doing what you expect by
>>
>> fitting with Gaussian random effects.
>>
>>
>> The random effects in glmer are Gaussian. The effects enter the model
>>
>> via the linear predictor, which is then translated to the mean function of
>> the chosen distribution of the response variable via the link function (NB
>> this is a conceptual explanation, not an algorithmic one).
>>
>>
>> Cheers,
>>
>> Andrew
>> On 20 Jul 2022, 6:40 PM +1000, J.D. Haltigan <jhaltiga using gmail.com>,
>>
>> wrote:
>>
>> External email: Please exercise caution
>>
>> ________________________________
>> Thanks, I will inspect the BLUPS.
>>
>> Re: choice of distribution, what I meant was, for example, if my fixed
>>
>> effects family is estimated using a Poisson model, does that inform the
>> choice for random effects as well? (i.e., why would one invoke a gaussian
>> distribution for random effects if the response variable is, say,
>> categorical, or a count?)
>>
>>
>> On Wed, Jul 20, 2022 at 3:40 AM Andrew Robinson <apro using unimelb.edu.au
>>
>> <mailto:apro using unimelb.edu.au>> wrote:
>>
>> You can use a qq-plot of the BLUPS to guide that decision.
>>
>> The difference might also be due to something else, though. Have you
>>
>> tried the call to hglm2 with a gaussian random effects family? Does it
>> give the same output as the glmer?
>>
>>
>> For: does the choice of distribution for the fixed effects portion of
>>
>> the model inform the choice for the random effects? I’m not sure what you
>> mean - do you mean the exponential family for the response variable?
>>
>>
>>
>> Cheers,
>>
>> Andrew
>>
>> --
>> Andrew Robinson
>> Chief Executive Officer, CEBRA and Professor of Biosecurity,
>> School/s of BioSciences and Mathematics & Statistics
>> University of Melbourne, VIC 3010 Australia
>> Tel: (+61) 0403 138 955
>> Email: apro using unimelb.edu.au<mailto:apro using unimelb.edu.au>
>> Website: https://researchers.ms.unimelb.edu.au/~apro@unimelb/
>>
>> I acknowledge the Traditional Owners of the land I inhabit, and pay my
>>
>> respects to their Elders.
>>
>> On 20 Jul 2022, 2:27 PM +1000, J.D. Haltigan <jhaltiga using gmail.com<mailto:
>>
>> jhaltiga using gmail.com>>, wrote:
>>
>> Hi:
>>
>> Is there clear best practice or guidance when it comes to choosing the
>> distribution of random effects where multiple choices exist (e.g.,
>> gaussian, gamma, etc.)? I ask in the context of extending some analyses
>> from an RCT in which the outcome is symptomatic seropositivity (so a
>> count). The random effects I am modeling are village cluster [union]
>>
>> (it's
>>
>> a cluster randomized trial). I get different results (significance-wise)
>> depending on whether I choose a normal or gamma distribution for the
>>
>> random
>>
>> effects.
>>
>> The basic model (proportional outcome) is:
>>
>> lme4_5_B = glmer(posXsymp~
>>
>> treatment+proper_mask_base+prop_resp_ill_base_2
>>
>> + pairID + (1 | union), family = "poisson", nAGQ=0, data =
>>
>> bdata.raw3)#lme4
>>
>> package using glmer
>>
>>
>> HGLM2_5_A = hglm2(posXsymp~
>>
>> treatment+proper_mask_base+prop_resp_ill_base_2
>>
>> + pairID + (1 | union), family =poisson (link = log), rand.family
>> =Gamma(link=log),
>> data = bdata.raw3)#HGLM package using hglm2
>>
>> Does, for example, the choice of distribution for the fixed effects
>>
>> portion
>>
>> of the model inform the choice for the random effects?
>>
>> Thank you for any insights.
>>
>> Best regards,
>> J.D.
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models using r-project.org<mailto:R-sig-mixed-models using r-project.org>
>>
>> mailing list
>>
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models<
>>
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>
>>
>>
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models using r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>
>> --
>> Dr. Benjamin Bolker
>> Professor, Mathematics & Statistics and Biology, McMaster University
>> Director, School of Computational Science and Engineering
>> (Acting) Graduate chair, Mathematics & Statistics
>>
>> _______________________________________________
>> R-sig-mixed-models using r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models using r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models using r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models using r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>> --
>> Dr. Benjamin Bolker
>> Professor, Mathematics & Statistics and Biology, McMaster University
>> Director, School of Computational Science and Engineering
>> (Acting) Graduate chair, Mathematics & Statistics
>>
>> _______________________________________________
>> R-sig-mixed-models using r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models using r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>

	[[alternative HTML version deleted]]