[R] glm() function not finding the maximum

[Richard Nixon]

For those of you who have been following my random ranting about the
glm(,family=gamma) function, for completeness here are my final
(hopefully!) remarks.

Thanks to Brett Presnell for explaining how the glm function estimates its
parameters.

The dispersion parameter in a gamma distribution is the same thing as
dispersion as defined by McCullagh and Nelder (Generalized linear models).
It is estimated by a moment estimator.

If one requires a maximum likelihood estimator, then one can use the
gamma.shape() function from the MASS library. shape=1/dispersion.

These two methods will generally yield similar results, but for highly
skewed data, like that given in my example, they can differ markedly.

Code showing this is included below.

Hope this is of some use to someone, and has iterated close enough to how
this function actually works.

Richard.


> gamma(data)

$loglik.mom
[1] 875.0035 #-2log likelihood using a moment estimator for shape

$loglik.mle
[1] 793.3913 #-2log likelihood using a mle estimator for shape

$par
[1]   0.09622557   0.51814501 415.99465517
#par[1] = moment estmator for shape
#par[2] = mle estimator for shape

data_c(51.47, 210.19, 49.55, 61.93, 60.61, 744.57, 338.59, 133.93,
191.57, 111.43, 432.83, 185.23, 155.61, 84.72, 120.2, 15.33,
77.05, 115.77, 25.23, 657.94, 108.39, 61.08, 142.42, 87.86, 272.87,
213.78, 65.23, 102.45, 58.16, 176.58, 76.58, 434.12, 362.35,
102.53, 103.6, 25.23, 97.19, 88.52, 118.55, 151.9, 2.7, 156.41,
21.79, 272.27, 23.16, 32.07, 6325.23, 92.37, 8340.04, 51.08,
55.59, 94.08, 69.98, 554.13, 104.88, 170.15, 945.1, 143.52)

#Fits data to a gamma distribution using glm()
gamma_function(data){
require(MASS, quietly=T)
model_glm(data~1, family=Gamma)
shape.mom_1/summary(model)$dispersion #this is a moment estimator
shape.mle_gamma.shape(m)$alpha        #this is a mle estimator
mu_mean(data)

dev.res.mom_-2*log(dgamma(data,shape=shape.mom,scale=mu/shape.mom))
loglik.mom_sum(dev.res.mom)  #actually -2 * log like
dev.res.mle_-2*log(dgamma(data,shape=shape.mle,scale=mu/shape.mle))
loglik.mle_sum(dev.res.mle)  #actually -2 * log like

list(loglik.mom=loglik.mom,loglik.mle=loglik.mle,par=c(shape.mom,shape.mle,mu))
}

--
Dr. Richard Nixon
MRC Biostatistics Unit, Institute of Public Health,
Robinson Way, Cambridge, CB2 2SR
http://www.mrc-bsu.cam.ac.uk/personal/richard
Tel: +44 (0)1223 330382, Fax: +44 (0)1223 330388

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._