How to correctly specify link function?

[Dunen]

I am new to Julia and MixedModels.jl. I am trying to fit a generalized mixed model where the link function is different to the typical for the error distribution family (to account for heteroscedasticity of variances), such as:

using DataFrames, MixedModels, StatsModels
gm1 = fit!(glmm(@formula(y ~ 1 + x + (1|group)), web_area, Gamma(IdentityLink)))

However, I get an Error: ERROR: MethodError: no method matching Gamma(::Type{IdentityLink})

I had a look in the MixedModels.jl code, and there it seems there needs to be another imput, as there is a comma infront of the link: glmm(f::Formula, fr::ModelFrame, d::Distribution[, l::GLM.Link])

...But I have no idea what that could be. How should I specify the link term so it is recognized? Any help is greatly appreciated!

[palday]

glmm is from an older version of MixedModels. You can try

fit(MixedModel, @formula(y ~ 1 + x + (1|group)), web_area, Gamma(), IdentityLink())

However, there are unfortunately known problems with gamma and inverse Gaussian models as well as Gaussian models with non-identity link (#609).

Beyond the technical issues -- are you using gamma with identity link to model reaction times (viz. Lo and Andrews, 2015)? I am generally skeptical of that advice and would tend to use a linear model, potentially with the response log or reciprocal transformed, depending on the particular data. (The reciprocal transformation is particularly nice because then you're modeling reaction speed, which is often the way we conceptually think of reaction time anyway.)

[Dunen]

Hi, thanks for the answer!

Unfortunally, I still get a different error: ERROR: UndefVarError: glmm not defined Stacktrace: [1] top-level scope

I am using the combination of Gamma and Identity link in order to get rid of heteroscedasticity in my model with linear relationship (not of reaction times, but changes of web size along an environmental gradient). I am actually unsure if there are any pitfalls in doing this (see this question I posted on stack exchange: https://stats.stackexchange.com/questions/589932/glmm-are-there-pitfalls-for-modeling-the-error-distribution-independently-of-t.

I tried understanding well your argument in the link you gave, but am afraid I am not deep enough in the topic to do so. Where can I find the list archives of R-SIG-mixed-models you were mentioning to understand the theoretical reasons to be skeptical better? Would a log-transformation of the response not change the relationship between the variables and thus not be correct if the only problem is heteroscedasticity? What would you recommend?

[palday]

Unfortunally, I still get a different error: ERROR: UndefVarError: glmm not defined Stacktrace: [1] top-level scope

Note that my answer did not use glmm and pointed out that it's from an older version of MixedModels, i.e. it's not present in the current release.

The link that I pointed to was about why gamma models in MixedModels.jl don't work, not about why not to use gamma in general. (If you fit a gamma model with MixedModels.jl, you will receive a warning to this effect.)

As to what choice of model family, link, etc. is good for your data -- it depends on your data and your inferential questions. If you can share your data, I might be able to help a little bit more, but I am not an expert in your domain (sounds like ecology/biology?). Heteroskedacity is a big topic, but based on your proposed solutions, are you more worried about skew?

[Dunen]

Hi,

It would be super nice if you could give me some input on this!
I have attached my data and a plot of the residuals.

The models I used in R notation are these:
lmer(web_area ~ rain + placement + canopy + understory + size + (1|altitude), data = web_areas)

glmer(web_area ~ rain + placement + canopy + understory + size + (1|altitude), data = web_areas, family = Gamma(link = "identity"))

The question is, if there is a within species trend of web_area with the specified independent variables (that is why I centered the independent variables).

Well, my main worry was the heteroscedasticity. Should I worry about the skew more? I think in the raw data there is some skew, but not along the regression line. Which family would you have picked to reduce heteroscedasticity?

Thanks a lot!

web_areas.csv

[palday]

A few background points

1. altitude is not in the provided data, but it is in the sample formulae.
2. it would be very unusual to have a continuous variable like altitude as a grouping variable (the bit behind |).
3. species would be a more typical grouping variable, so I tried that.
4. if you want to fit a gamma model, you can't do it with MixedModels.jl

With all that in mind, I looked at a few different standard transformations and log(web_area) worked quite well. sqrt(web_area) was also okay, after playing around with a few different possible random-effects structures. Both of these make sense as transformations because area grows as the square of radius (assuming a somewhat circular web, but a similar argument holds for size length in a more square web), so transformations that that account for that superlinear growth make sense. Yes, transforming the dependent variable changes the interpretation, but so does changing the model family (e.g., we often implicitly think of the mean, median and mode as being very similar, but they can be quite different for non-Gaussian distributions). You can see my tinkerings here.

NB: I am not an expert in ecology and I certainly don't know your data and research questions as well as you do. You are the person who has to justify a particular model structure as being appropriate for your desired inferences. I just took a quick look one evening after my day job.

[Dunen]

Hi,
thanks a lot for taking a look and taking the time to play around with my data!

Yes, you are correct, it is supossed to be species, not altitute! I don't know why I wrote altitute...

Thanks for your input, your argument on transformations that account for superlinear growth makes sense. I tried the log transformation before, and the results in terms of significance testing are basically the same. However, it seems to me that when log transformed, the linear relationship doesn't seem to fit perfectly (e.g. in the fitted vs response plot in your tinkers, more residuals are below the line for values <1.5, but >1.5 there are more above the line). With a typical transformation, I transform both, how big the spread is around the fitted line, and how the fitted line looks... That is why I went for just transforming the spread of the errors without changing the relation between independent and dependent variables themselfs.

That said, I will look more into the differences in interpretation that you mentioned, and maybe just use the normal log-transformation then.

Thanks and cheers!