I am looking for an example to determine the coefficients of a GLM by maximizing the (log)likelihood, preferrably using R. The optimization (search algorithm) is not the problem; this can be done visually by inspecting plots. My key problem is how to calculate the likelihood values.
Consider some count data and a Poisson model.
For the simplest model, the calculation of the likelihood is staight foreward: the sum of the log densities of the poisson distribution and a plot of logLik versus lambda shows the maximum. This is the R solution for Y ~ 1:
n
I made a review of your code and it looks good to me,are you using R or Splus coding.it is strange it should work correctally i do not see any error if so it should be in the internal code of the loglike command
sorry
You may create a very simple program to construct a synthetic Poisson data set with coefficient values for predictors you define. Here I made two predictors as random uniform, but they can be binary, factor, or whatever. If you run this script you'll have a data set which is modeled using a Poisson family with an intercept of 2 and coefficients of 0.75 and -1.25 respectively for predictos x1 and x2. Its easy to create a model with offsets as well. There are other ways to do this by programming the Poisson component from first principles. I address this in detail in my book, Negative Binomial Regression, 2nd edition, 2011, Cambridge Univ Press, which discusses in detail all count models, including Bayesian count models.
At the bottom of the code is a way to calculate the log-likelihood following the glm. Unfortunately the statistic is not displayed in model output. You can use saved statistics to calulate it. Best, Joseph Hilbe
library(MASS)
nobs
Why don´t you use Eviews? An econometric program instead? It is easier than coding.
Thanks for your answers so far.
I do NOT look for a way to calculate GLMs in general with a computer software. I understood (at least I hope so) how the parameter of a Poisson-distributed variable is estimated by maximum likelihood (equally for a binomial or normally distributed variable; this doesn't matter). This was explored in the first part of the R script. Just to show: it works.
My problem is to use a similar stupid, "pedestrian way" to get the betas of a model where lambda (or pi or mu) is expressed as X*beta (the case when X is not just (1,1,1,...)). I know that there are software doing this smotthly and efficiently. But I want to see/explore/understand the WAY of the solution as used for the most simple model. I tried this in the second part of the R code and expected to get a result alt least quite similar to the result of the glm function, but this is not the case.
@Joseph: You showed that I should use glm. I know hoe to use glm. This does not answer the question I had.
@Carl: I am not sure if I understood your answer. What data is used to fit the model is not interesting here. It doesnt matter if the data suits the model or vice versa. The results (however sensible) should be similar for the "manual approach" as in the output from the glm. The coefficients calculated by glm are log-linked. You can see that exp(coef(m)) *equals* the given value of lambda (5) to generate the values (with rpois).
@Marco: Again, I'm not looking for a software. I'm looking for understanding.
Hi Jochem,
In R, indices and hence levels starts at 1, and not at 0 like in other programming languages, C for instance.
Hence, in your log.Lik function, you have nothing with y[g==0] and get the values for the first level with y[g==1].
Use y[g==1] then y[g==2] and it should work.
Best regards,
(Edit: Well, helps to have better figures, but still does not give the point on the top of the contours...)
In fact, you are not clear on what is beta I think.
In your density function, beta is the difference between the two parameters, so it should be 20 - 5 = 15.
When you are drawing the contours, with X=lambda and Y=beta, you consequently also uses this convention.
But when you place the glm point, you give to its X the sum( coef( m ) ) value, that is you have 20 and not 15, lambda + beta and not beta. Hence the shift.
If you replace the last line with
points(exp(coef(m)[1]),exp(sum(coef(m)[1:2])) - exp(coef(m)[1]),pch=15,col=3)
which gives to Y really the beta value as the difference between the two values, it is correctly placed at the peak of the likelyhood function.
Best regards,
- Badges
- Science topic
- Similar topics
- Mathematics
- Statistics