3
$\begingroup$

Let's say I have an outcome that is exponentially distributed, so $p(y|\lambda_i) = \lambda_i e^{-\lambda_i y}$. However, we also know that $\log{\lambda_i} = \beta_0+\beta_1x_i$. If we want to estimate the MLE's for $\beta_0$ and $\beta_1$, what's the best way to do it given data?

I was thinking of using the following code:

summary(glm(y~x,family=Gamma(link="log")))

Not 100% sure what the various arguments are going though.

$\endgroup$
  • $\begingroup$ your idea won't restrict the shape parameter of the Gamma distribution to 1, which would be necessary for your plan ... bbmle::mle2 will do this. You could try to derive and solve the score equations but I don't know how far you'd get ... $\endgroup$ – Ben Bolker Oct 16 '16 at 0:34
  • $\begingroup$ How would I implement that into the code? Also I found the score equations but I don't think there is a closed form solutions of the estimates $\endgroup$ – Lzydude Oct 16 '16 at 0:35
  • $\begingroup$ You don't estimate maximum likelihood estimators ... unless you are constructing an estimate-of-an-estimate, you calculate MLEs. The MLEs themselves are usually estimators of parameters, or functions of them $\endgroup$ – Glen_b -Reinstate Monica Oct 16 '16 at 0:39
  • $\begingroup$ It's not clear what you mean in your final sentence. Did you mean 'doing' rather than 'going'? Are you asking for an explanation of the arguments of the function? Beyond the help, probably the best thing to do is read a basic treatment of GLMs. $\endgroup$ – Glen_b -Reinstate Monica Oct 16 '16 at 0:46
  • $\begingroup$ @Lzydude, "how do I ... in R?" is really a StackOverflow question rather than a CrossValidated question. bbmle::mle2(y~dexp(exp(lograte)), parameters=list(lograte~x), start=list(logmu=0), data=...) should do it, I think ... $\endgroup$ – Ben Bolker Oct 16 '16 at 0:46
7
$\begingroup$

This has been answered on the R help list by Adelchi Azzalini: the important point is that the dispersion parameter (which is what distinguishes an exponential distribution from the more general Gamma distribution) does not affect the parameter estimates in a generalized linear model, only the standard errors of the parameters/confidence intervals/p-values etc.; in R an estimate of the dispersion parameter is automatically reported, but as Azzalini comments, summary.glm allows the user to specify the dispersion parameter. So, as stated by Azzalini,

The Gamma family is parametrised in glm() by two parameters: mean and dispersion; the "dispersion" regulates the shape.

So [one] must fit a GLM with the Gamma family, and then produce a "summary" with dispersion parameter set equal to 1, since this value corresponds to the exponential distribution in the Gamma family.

In practice:

fit <- glm(formula =...,  family = Gamma(link="log"))
summary(fit,dispersion=1)   

[Azzalini had family=Gamma, i.e. using the default inverse link; I changed it to specify the log link as in your question.]

$\endgroup$
1
$\begingroup$

Using a GLM call as you suggest there is the easiest correct approach, but to actually make the Gamma into an exponential you can specify the dispersion to be 1. It will not change the fitted mean / coefficients, but it impacts the standard errors.

[i.e. Your suggested call of summary(glm(y~x,family=Gamma(link="log"))) should give you what you want, but if you're interested in significance of coefficients and so on under the exponential assumption, you'd add ,dispersion=1 before the final parenthesis. If you want to fit Gamma GLMs more generally, there are several useful helper functions in the package MASS that comes with R but is not loaded by default.]

(Another alternative might be to consider using parametric survival models which also offer ways to fit exponential, Weibull and various other models.)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.