Hello all, I was wondering if someone can enlighten me as to the difference between the logLik in R vis-a-vis Stata for a GLM model with the gamma family. Stata calculates the loglikelihood of the model as (in R notation) some equivalent function of -1/scale * sum(Y/mu+log(mu)+(scale-1)*log(Y)+log(scale)+scale*lgamma(1/scale)) where scale (or dispersion) = 1, Y = the response variable, and mu is the fitted values given by the fitted model. R seems to use a very similar approach, but scale is set equal to the calculated dispersion for the gamma model. However, when I calculate the logLik by hand this way the answer differs slightly (about .5) from the logLik(glm.m1). I haven't been able to figure out why looking at the help. If anyone has any ideas, the insight would be much appreciated. Cheers, Skipper

jseabold wrote:> > Hello all, > > I was wondering if someone can enlighten me as to the difference > between the logLik in R vis-a-vis Stata for a GLM model with the gamma > family. > > Stata calculates the loglikelihood of the model as (in R notation) > some equivalent function of > > -1/scale * > sum(Y/mu+log(mu)+(scale-1)*log(Y)+log(scale)+scale*lgamma(1/scale)) > > where scale (or dispersion) = 1, Y = the response variable, and mu is > the fitted values given by the fitted model. > > R seems to use a very similar approach, but scale is set equal to the > calculated dispersion for the gamma model. However, when I calculate > the logLik by hand this way the answer differs slightly (about .5) > from the logLik(glm.m1). > > I haven't been able to figure out why looking at the help. If anyone > has any ideas, the insight would be much appreciated. > >This is not a full answer, but the beginning of an answer about how to find out the answer. stats:::logLik.glm shows that R computes the log-likelihood (a bit oddly) by extracting the AIC component of a model, dividing by 2, and subtracting it from the number of parameters (since AIC = -2 (logLik) + 2 p). Gamma()$aic in turn gives the formula for the AIC of a gamma GLM: function (y, n, mu, wt, dev) { n <- sum(wt) disp <- dev/n -2 * sum(dgamma(y, 1/disp, scale = mu * disp, log = TRUE) * wt) + 2 } and ?dgamma gives the The Gamma distribution with parameters 'shape' = a and 'scale' = s has density f(x)= 1/(s^a Gamma(a)) x^(a-1) e^-(x/s) By the way, I find it deeply confusing that Stata refers to what R calls the "shape" parameter (a, the thing that appears inside the Gamma or lgamma() function) as "shape" <http://www.stata.com/help.cgi?glm>. Oh well. good luck Ben Bolker -- View this message in context: http://www.nabble.com/GLM-Gamma-Family-logLik-formula--tp24504030p24508346.html Sent from the R help mailing list archive at Nabble.com.