[R] simultaneously maximizing two independent log likelihood functions using mle2

[Adam Zeilinger]

Dear Dr. Bolker and R list,

Thanks so much for the help!  Apologies for taking so long to write 
back, I've been distracted with other work.  Between your response and 
the suggested section in the EMD book, I was able to figure it out!

I have another question.  I ran three optimization runs with my 
likelihood function: one run in which I simultaneously maximized the 
function with both years of my data, in the way described in Section 
6.6.1 in the EMD book (lets call this mle2 run NLL.both), one run in 
which I maximize the function with only first year (NLL1), and a final 
run in which I maximize the function with only the second year of data 
(NLL2).

The sum of the log-likelihood values of the second two models should 
equal the log-likelihood value of the first model, correct?  In other 
words, NLL1 + NLL2 = NLL.both?  I do not get equivalent values.  NLL1 + 
NLL2 = -339.07, whereas NLL.both = -481.73.  However, the parameter 
estimates are equivalent.  Do you have any thoughts on what is going on 
here?  Should I be concerned about unequal log-likelihood values though 
my parameter estimates are equivalent?

I tried the same thing with the DBH example from Section 6.6.1 (from 
data set FirDBHec).  The log-likelihood values were equivalent when I 
use the following code:

data(FirDBHFec)
X = na.omit(FirDBHFec[, c("TOTCONES", "DBH", "WAVE_NON")])
X$TOTCONES = round(X$TOTCONES)

# Both wave groups simultaneously maximized
nbNLL.both = function(a.w, b.w, a.n, b.n, k) {
   wcode = as.numeric(X$WAVE_NON)
   a = c(a.n, a.w)[wcode]
   b = c(b.n, b.w)[wcode]
   predcones = a * DBH^b
   -sum(dnbinom(TOTCONES, mu = predcones, size = k,
   log = TRUE))
}

start.both = list(a.n = 0.28, a.w = 0.4, b.n = 2.36, b.w = 2.15, k = 
1.51) # Starting values from EMD book Section 6.6.1
nbfit.both = mle2(nbNLL.both, start = start.ab2, data = X,
   parameters = list(a ~ WAVE_NON, b ~ WAVE_NON))

# Separate MLE for each wave group
w.group = X[X$WAVE_NON == "w",]
nw.group = X[X$WAVE_NON == "n",]

nbNLL.bywave = function(a, b, k) {
   predcones = a * DBH^b
   -sum(dnbinom(TOTCONES, mu = predcones, size = k,
   log = TRUE))
}

start.bywave = list(a = 0.30, b = 2.32, k = 1.50)
nbfit.wave = mle2(nbNLL.bywave, start = start.bywave, data = w.group)
nbfit.nonwave = mle2(nbNLL.bywave, start = start.bywave data = nw.group)

I am using a similar approach when maximizing my function as I've 
described with the DBH example.  Any help would be much appreciated.

Adam Zeilinger



On 10/17/2011 5:47 AM, Ben Bolker wrote:
> Adam Zeilinger<zeil0006<at>  umn.edu>  writes:
>
>> I have a log likelihood function that I was able to optimize using
>> mle2.  I have two years of the data used to fit the function and I would
>> like to fit both years simultaneously to test if the model parameter
>> estimates differ between years, using likelihood ratio tests and AIC.
>> Can anyone give advice on how to do this?
>>
>> My likelihood functions are long so I'll use the tadpole predation
>> example from Ben Bolker's book, Ecological Data and Models in R (p.
>> 268-270).
>>
>> library(emdbook)
>> data(ReedfrogFuncresp)
>> attach(ReedfrogFuncresp)
>> # Holling Type II Equation
>> holling2.pred = function(N0, a, h, P, T) {
>>     a * N0 * P * T/(1 + a * h * N0)
>> }
>> # Negative log likelihood function
>> NLL.holling2 = function(a, h, P = 1, T = 1) {
>>     -sum(dbinom(Killed, prob = a * T * P/(1 + a * h * Initial),
>>     size = Initial, log = TRUE))
>> }
>> # MLE statement
>> FFR.holling2 = mle2(NLL.holling2, start = list(a = 0.012,
>>     h = 0.84), data = list(T = 14, P = 3))
>>
>> I have my negative log likelihood function setup similarly to the above
>> example.  Again, my goal is to simultaneously estimate parameters from
>> the same function for two years, such that I can test if the parameters
>> from the two years are different.  Perhaps an important difference from
>> the above example is that I am using a multinomial distribution (dmnom)
>> because my data are trinomially distributed.
>    The most convenient way to do this would be to use the
> 'parameters' argument to mle2.  For example,
>
> FFR.holling2.byyear<- mle2(NLL.holling2, start = list(a = 0.012,
>      h = 0.84), data = c(list(T = 14, P = 3), ReedfrogFuncresp),
>      parameters=list(a~year, h~year))
>
> -- or you could fit a model where only one or the other of the
> parameters varied between years.
>     R should automatically construct a set of parameters (in
> this case 'a for year 1', 'delta-a between year 1 and year 2',
> 'h for year 1', 'delta-h between year 1 and year 2'); it
> sets the starting values of the 'delta' variables to zero
> by default.  (Things get a little bit tricky if you also
> want to set lower and upper bounds.)  [See section 6.6.1 of
> the book for more details.]
>
>     Note that I have included ReedfrogFuncresp in the data list.
> I now regret using attach() in the book examples.
>
>    Ben Bolker
>
> ______________________________________________
> R-help at r-project.org  mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guidehttp://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Adam Zeilinger
Post Doctoral Scholar
Department of Entomology
University of California Riverside
www.linkedin.com/in/adamzeilinger