The McCullagh (1989) distribution has density function
$$f(y;\theta,\nu) =
\frac{ \{ 1-y^2 \}^{\nu-\frac12}}
{ (1-2\theta y + \theta^2)^{\nu} \mbox{Beta}(\nu+\frac12, \frac12)}$$
where \(-1 < y < 1\) and \(-1 < \theta < 1\).
This distribution is equation (1) in that paper.
The parameter \(\nu\) satisfies \(\nu > -1/2\),
therefore the default is to use an log-offset link
with offset equal to 0.5, i.e.,
\(\eta_2=\log(\nu+0.5)\).
The mean is of \(Y\) is \(\nu \theta / (1+\nu)\),
and these are returned as the fitted values.
This distribution is related to the Leipnik distribution (see Johnson
et al. (1995)), is related to ultraspherical functions, and under
certain conditions, arises as exit distributions for Brownian motion.
Fisher scoring is implemented here and it uses a diagonal matrix so
the parameters are globally orthogonal in the Fisher information sense.
McCullagh (1989) also states that, to some extent, \(\theta\)
and \(\nu\) have the properties of a location parameter and a
precision parameter, respectively.