R help - Dec 2015 - package MASS - MLE of negative binomial distributions

Dear list,



I have a question about the exact estimate of the maximum likelihood for a
negative binomial fit. I'm trying to approach this in two different ways:
the first one is a fit using the glm.nb method, and the second one is a fit
using the fitdistr function for each condition separately, where I add up all
log likelihoods. These two methods do not yield the same values for the log
likelihood of the fit. They do yield the same log likelihood if all data are one
group (no summation), so I assume I'm doing something wrong when I sum up
log likelihoods. Am I not "allowed" to do this?


Example code:
library(MASS)
x<-c(601,619,637,609,594,499,494,507,477,450,400,367,428,359,400,276,260,262,304,342,216,189,152,231,200,104,85,85,85,112)
groups<-as.factor(c(rep("dist1",5),rep("dist2",5),rep("dist3",5),rep("dist4",5),rep("dist5",5),rep("dist6",5)))

glm.nb(x~groups)$twologlik

logliks<-NULL
for(group in levels(groups))
{
	NBfit<-fitdistr(x[groups==group],"Negative Binomial")
	logliks<-c(logliks,NBfit$loglik)
	rm(NBfit)
}

sum(logliks)*2


Many thanks!
Rik

glm.nb fits 6 mean parameters plus 1 theta. 6 x fitdistr fits two parameters
each.

-pd

On 10 Dec 2015, at 16:21 , Rik Verdonck <Rik.Verdonck at bio.kuleuven.be>
wrote:
> Dear list,
> 
> 
> 
> I have a question about the exact estimate of the maximum likelihood for a
negative binomial fit. I'm trying to approach this in two different ways:
the first one is a fit using the glm.nb method, and the second one is a fit
using the fitdistr function for each condition separately, where I add up all
log likelihoods. These two methods do not yield the same values for the log
likelihood of the fit. They do yield the same log likelihood if all data are one
group (no summation), so I assume I'm doing something wrong when I sum up
log likelihoods. Am I not "allowed" to do this?
> 
> 
> Example code:
> library(MASS)
>
x<-c(601,619,637,609,594,499,494,507,477,450,400,367,428,359,400,276,260,262,304,342,216,189,152,231,200,104,85,85,85,112)
>
groups<-as.factor(c(rep("dist1",5),rep("dist2",5),rep("dist3",5),rep("dist4",5),rep("dist5",5),rep("dist6",5)))
> 
> glm.nb(x~groups)$twologlik
> 
> logliks<-NULL
> for(group in levels(groups))
> {
> 	NBfit<-fitdistr(x[groups==group],"Negative Binomial")
> 	logliks<-c(logliks,NBfit$loglik)
> 	rm(NBfit)
> }
> 
> sum(logliks)*2
> 
> 
> Many thanks!
> Rik
> 
> 
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com

That makes sense indeed, but it also raises further questions to me if you
don't mind:
- Is there a good reason why the default assumption would be a constant
overdispersion?  Is it because you lose fewer degrees of freedom which would
give you more power in comparing models?
- Is there a way to tell glm.nb to fit the thetas separately per group?
- Is there a way to fix theta (i.e. feed it a value) and only make an estimate
of the expected value?

Many thanks!
Rik


________________________________________
Van: peter dalgaard [pdalgd at gmail.com]
Verzonden: donderdag 10 december 2015 16:53
Aan: Rik Verdonck
CC: r-help at r-project.org
Onderwerp: Re: [R] package MASS - MLE of negative binomial distributions

glm.nb fits 6 mean parameters plus 1 theta. 6 x fitdistr fits two parameters
each.

-pd

On 10 Dec 2015, at 16:21 , Rik Verdonck <Rik.Verdonck at bio.kuleuven.be>
wrote:
> Dear list,
>
>
>
> I have a question about the exact estimate of the maximum likelihood for a
negative binomial fit. I'm trying to approach this in two different ways:
the first one is a fit using the glm.nb method, and the second one is a fit
using the fitdistr function for each condition separately, where I add up all
log likelihoods. These two methods do not yield the same values for the log
likelihood of the fit. They do yield the same log likelihood if all data are one
group (no summation), so I assume I'm doing something wrong when I sum up
log likelihoods. Am I not "allowed" to do this?
>
>
> Example code:
> library(MASS)
>
x<-c(601,619,637,609,594,499,494,507,477,450,400,367,428,359,400,276,260,262,304,342,216,189,152,231,200,104,85,85,85,112)
>
groups<-as.factor(c(rep("dist1",5),rep("dist2",5),rep("dist3",5),rep("dist4",5),rep("dist5",5),rep("dist6",5)))
>
> glm.nb(x~groups)$twologlik
>
> logliks<-NULL
> for(group in levels(groups))
> {
>       NBfit<-fitdistr(x[groups==group],"Negative Binomial")
>       logliks<-c(logliks,NBfit$loglik)
>       rm(NBfit)
> }
>
> sum(logliks)*2
>
>
> Many thanks!
> Rik
>
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com