R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Estimate the Shape Parameter of the Gamma Distribution in a...

gamma.shape

R Documentation

Estimate the Shape Parameter of the Gamma Distribution in a GLM Fit

Description

Find the maximum likelihood estimate of the shape parameter of the gamma distribution after fitting a Gamma generalized linear model.

Usage

## S3 method for class 'glm'
gamma.shape(object, it.lim = 10,
            eps.max = .Machine$double.eps^0.25, verbose = FALSE, ...)

Arguments

`object`	Fitted model object from a `Gamma` family or `quasi` family with `variance = "mu^2"`.
`it.lim`	Upper limit on the number of iterations.
`eps.max`	Maximum discrepancy between approximations for the iteration process to continue.
`verbose`	If `TRUE`, causes successive iterations to be printed out. The initial estimate is taken from the deviance.
`...`	further arguments passed to or from other methods.

Details

A glm fit for a Gamma family correctly calculates the maximum likelihood estimate of the mean parameters but provides only a crude estimate of the dispersion parameter. This function takes the results of the glm fit and solves the maximum likelihood equation for the reciprocal of the dispersion parameter, which is usually called the shape (or exponent) parameter.

Value

List of two components

`alpha`	the maximum likelihood estimate
`SE`	the approximate standard error, the square-root of the reciprocal of the observed information.

References

Venables, W. N. and Ripley, B. D. (2002) Modern Applied Statistics with S. Fourth edition. Springer.

Examples

clotting <- data.frame(
    u = c(5,10,15,20,30,40,60,80,100),
    lot1 = c(118,58,42,35,27,25,21,19,18),
    lot2 = c(69,35,26,21,18,16,13,12,12))
clot1 <- glm(lot1 ~ log(u), data = clotting, family = Gamma)
gamma.shape(clot1)

gm <- glm(Days + 0.1 ~ Age*Eth*Sex*Lrn,
          quasi(link=log, variance="mu^2"), quine,
          start = c(3, rep(0,31)))
gamma.shape(gm, verbose = TRUE)
summary(gm, dispersion = gamma.dispersion(gm))  # better summary

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(MASS)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/MASS/gamma.shape.glm.Rd_%03d_medium.png", width=480, height=480)
> ### Name: gamma.shape
> ### Title: Estimate the Shape Parameter of the Gamma Distribution in a GLM
> ###   Fit
> ### Aliases: gamma.shape gamma.shape.glm print.gamma.shape
> ### Keywords: models
> 
> ### ** Examples
> 
> clotting <- data.frame(
+     u = c(5,10,15,20,30,40,60,80,100),
+     lot1 = c(118,58,42,35,27,25,21,19,18),
+     lot2 = c(69,35,26,21,18,16,13,12,12))
> clot1 <- glm(lot1 ~ log(u), data = clotting, family = Gamma)
> gamma.shape(clot1)
               
Alpha: 538.1315
SE:    253.5991
> 
> gm <- glm(Days + 0.1 ~ Age*Eth*Sex*Lrn,
+           quasi(link=log, variance="mu^2"), quine,
+           start = c(3, rep(0,31)))
> gamma.shape(gm, verbose = TRUE)
Initial estimate: 1.060344
Iter. 1 Alpha: 1.238408
Iter. 2 Alpha: 1.276997
Iter. 3 Alpha: 1.278343
Iter. 4 Alpha: 1.278345
                
Alpha: 1.2783449
SE:    0.1345175
> summary(gm, dispersion = gamma.dispersion(gm))  # better summary

Call:
glm(formula = Days + 0.1 ~ Age * Eth * Sex * Lrn, family = quasi(link = log, 
    variance = "mu^2"), data = quine, start = c(3, rep(0, 31)))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.0385  -0.7164  -0.1532   0.3863   1.3087  

Coefficients: (4 not defined because of singularities)
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)            3.06105    0.44223   6.922 4.46e-12 ***
AgeF1                 -0.61870    0.59331  -1.043 0.297041    
AgeF2                 -2.31911    0.98885  -2.345 0.019014 *  
AgeF3                 -0.37623    0.53149  -0.708 0.479020    
EthN                  -0.13789    0.62540  -0.220 0.825496    
SexM                  -0.48844    0.59331  -0.823 0.410369    
LrnSL                 -1.92965    0.98885  -1.951 0.051009 .  
AgeF1:EthN             0.10249    0.82338   0.124 0.900942    
AgeF2:EthN            -0.50874    1.39845  -0.364 0.716017    
AgeF3:EthN             0.06314    0.74584   0.085 0.932534    
AgeF1:SexM             0.40695    0.94847   0.429 0.667884    
AgeF2:SexM             3.06173    1.11626   2.743 0.006091 ** 
AgeF3:SexM             1.10841    0.74208   1.494 0.135267    
EthN:SexM             -0.74217    0.82338  -0.901 0.367394    
AgeF1:LrnSL            2.60967    1.10114   2.370 0.017789 *  
AgeF2:LrnSL            4.78434    1.36304   3.510 0.000448 ***
AgeF3:LrnSL                 NA         NA      NA       NA    
EthN:LrnSL             2.22936    1.39845   1.594 0.110899    
SexM:LrnSL             1.56531    1.18112   1.325 0.185077    
AgeF1:EthN:SexM       -0.30235    1.32176  -0.229 0.819065    
AgeF2:EthN:SexM        0.29742    1.57035   0.189 0.849780    
AgeF3:EthN:SexM        0.82215    1.03277   0.796 0.425995    
AgeF1:EthN:LrnSL      -3.50803    1.54655  -2.268 0.023311 *  
AgeF2:EthN:LrnSL      -3.33529    1.92481  -1.733 0.083133 .  
AgeF3:EthN:LrnSL            NA         NA      NA       NA    
AgeF1:SexM:LrnSL      -2.39791    1.51050  -1.587 0.112400    
AgeF2:SexM:LrnSL      -4.12161    1.60698  -2.565 0.010323 *  
AgeF3:SexM:LrnSL            NA         NA      NA       NA    
EthN:SexM:LrnSL       -0.15305    1.66253  -0.092 0.926653    
AgeF1:EthN:SexM:LrnSL  2.13480    2.08685   1.023 0.306317    
AgeF2:EthN:SexM:LrnSL  2.11886    2.27882   0.930 0.352473    
AgeF3:EthN:SexM:LrnSL       NA         NA      NA       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasi family taken to be 0.7822615)

    Null deviance: 190.40  on 145  degrees of freedom
Residual deviance: 128.36  on 118  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 7

> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>