Leipnik Regression Family Function

Estimates the two parameters of a (transformed) Leipnik distribution by maximum likelihood estimation.

leipnik(lmu = "logitlink", llambda = "logofflink(offset = 1)",
        imu = NULL, ilambda = NULL)

Arguments

lmu, llambda

Link function for the $\mu$ and $\lambda$ parameters. See Links for more choices.

imu, ilambda

Numeric. Optional initial values for $\mu$ and $\lambda$.

Details

The (transformed) Leipnik distribution has density function $$f(y;\mu,\lambda) = \frac{ \{ y(1-y) \}^{-\frac12}}{ \mbox{Beta}( \frac{\lambda+1}{2}, \frac12 )} \left[ 1 + \frac{(y-\mu)^2 }{y(1-y)} \right]^{ -\frac{\lambda}{2}}$$ where $0 < y < 1$ and $\lambda > -1$. The mean is $\mu$ (returned as the fitted values) and the variance is $1/\lambda$.

Jorgensen (1997) calls the above the transformed Leipnik distribution, and if $y = (x+1)/2$ and $\mu = (\theta+1)/2$, then the distribution of $X$ as a function of $x$ and $\theta$ is known as the the (untransformed) Leipnik distribution. Here, both $x$ and $\theta$ are in $(-1, 1)$.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm, rrvglm and vgam.

References

Jorgensen, B. (1997). The Theory of Dispersion Models. London: Chapman & Hall

Johnson, N. L. and Kotz, S. and Balakrishnan, N. (1995). Continuous Univariate Distributions, 2nd edition, Volume 2, New York: Wiley. (pages 612–617).

Author

T. W. Yee

Note

Convergence may be slow or fail. Until better initial value estimates are forthcoming try assigning the argument ilambda some numerical value if it fails to converge. Currently, Newton-Raphson is implemented, not Fisher scoring. Currently, this family function probably only really works for intercept-only models, i.e., y ~ 1 in the formula.

Examples

ldata <- data.frame(y = rnorm(2000, 0.5, 0.1))  # Improper data
fit <- vglm(y ~ 1, leipnik(ilambda = 1), ldata, trace = TRUE)
#> Iteration 1: loglikelihood = -41797.9896
#> Iteration 2: loglikelihood = -11775.9087
#> Iteration 3: loglikelihood = -1325.55842
#> Iteration 4: loglikelihood = 1984.60136
#> Iteration 5: loglikelihood = 2800.48874
#> Iteration 6: loglikelihood = 2907.19904
#> Iteration 7: loglikelihood = 2910.45711
#> Iteration 8: loglikelihood = 2910.46141
#> Iteration 9: loglikelihood = 2910.46141
head(fitted(fit))
#>           [,1]
#> [1,] 0.5020428
#> [2,] 0.5020428
#> [3,] 0.5020428
#> [4,] 0.5020428
#> [5,] 0.5020428
#> [6,] 0.5020428
with(ldata, mean(y))
#> [1] 0.50188
summary(fit)
#> 
#> Call:
#> vglm(formula = y ~ 1, family = leipnik(ilambda = 1), data = ldata, 
#>     trace = TRUE)
#> 
#> Coefficients: 
#>               Estimate Std. Error z value Pr(>|z|)    
#> (Intercept):1 0.008171   0.009743   0.839    0.402    
#> (Intercept):2 3.173675   0.030982 102.436   <2e-16 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Names of linear predictors: logitlink(mu), logofflink(lambda, offset = 1)
#> 
#> Log-likelihood: 2910.461 on 3998 degrees of freedom
#> 
#> Number of Fisher scoring iterations: 9 
#> 
#> No Hauck-Donner effect found in any of the estimates
#> 
coef(fit, matrix = TRUE)
#>             logitlink(mu) logofflink(lambda, offset = 1)
#> (Intercept)   0.008171442                       3.173675
Coef(fit)
#>         mu     lambda 
#>  0.5020428 22.8951413 

sum(weights(fit))  # Sum of the prior weights
#> [1] 2000
sum(weights(fit, type = "work"))  # Sum of the working weights
#> [1] 11576.51