Feature request: skew-normal family

[psmorris15]

The current application I am working on requires use of a skew-normal likelihood. The skew-normal family is useful for modeling unimodal responses that lie on the real line but for which the normal distribution is inappropriate due to the presence of skewness. I have added the distribution to my local glmmTMB version following the instructions given in the Hacking glmmTMB vignette.

Would there be any interest in working together to add the skew-normal family to the official glmmTMB release? While the brms R package has implemented the skew-normal, I am not aware of any R packages that allow for its use in a frequentist setting. Such an addition should attract users who need to model their data with a skew-normal likelihood but prefer to use frequentist methods or require a faster runtime than brms can provide.

[bbolker]

This seems like a good idea. The easiest way to proceed would probably be for you to make a fork of this repository, post your implementation publicly, then make a pull request.

Were there any surprises/challenges in implementing the SN, or did the vignette give you the information you needed?

[psmorris15]

I appreciate the quick reply. For the record, I've submitted this under PR #1029.

The vignette had the bulk of the information I needed. A few things:

• I was not sure how to implement the simulation step for the skew-normal. I didn't have an existing function like rnorm or rgamma to fall back on for this. This is potentially something that can be worked out.
• I was able to figure out how to get the skew parameter to display on print but not when the summary function is called.
• In my experience the implementation can be sensitive to the starting values to a larger extent than for the simpler distributions like normal or gamma. The default values often did not work (would get the error "negative log-likelihood is NaN at starting parameter values"), particularly for complex models, requiring me to make reasonable guesses.

[bbolker]

Poking around a little bit, this is the R code from fGarch:::.rsnorm (which generates a standardized skew-normal that then needs to be shifted and scaled). This "just" needs to be translated into C++ (need to be a little careful with the abs/sign stuff, which is problematic in an autodiff context, but I think should be OK here; we might need to implement with floating point rather than fancy AD types). See e.g. rtweedie in https://github.com/glmmTMB/glmmTMB/blob/master/glmmTMB/src/distrib.h

    weight = xi/(xi + 1/xi)
    z = runif(n, -weight, 1 - weight)
    Xi = xi^sign(z)
    Random = -abs(rnorm(n))/Xi * sign(z)
    m1 = 2/sqrt(2 * pi)
    mu = m1 * (xi - 1/xi)
    sigma = sqrt((1 - m1^2) * (xi^2 + 1/xi^2) + 2 * m1^2 - 1)
    Random = (Random - mu)/sigma
    Random

As we expand to more cases I definitely think we need to be more thoughtful about starting values ...