R help - Oct 2012 - model selection with spg and AIC (or, convert list to fitted model object)

Dear R Help,

I have two nested negative log-likelihood functions that I am optimizing 
with the spg function [BB package].  I would like to perform model 
selection on these two objective functions using AIC (and possibly 
anova() too).  However, the spg() function returns a list and I need a 
fitted model object for AIC(), ICtab() [bbmle package], or anova().

How can I perform AIC-based model selection on two spg-optimized 
objective functions?  Alternatively, how can I convert the list returned 
by spg into a fitted model object that can be run in AIC, ICtab, or anova?

Here is an example:

##########################################################################

library(BB)
library(bbmle)

# define multinomial distribution
dmnom2 <- function(x,prob,log=FALSE) {
   r <- lgamma(sum(x) + 1) + sum(x * log(prob) - lgamma(x + 1))
   if (log) r else exp(r)
}


# data frame y
y <- structure(list(tv = 1:20, n1 = c(17L, 18L, 16L, 18L, 16L, 16L,
                18L, 18L, 20L, 16L, 20L, 16L, 17L, 17L, 18L, 19L, 18L, 
16L, 17L,
                17L), n2 = c(2L, 2L, 3L, 2L, 4L, 4L, 2L, 2L, 0L, 4L, 0L, 
4L,
                3L, 3L, 2L, 1L, 2L, 4L, 3L, 3L), n3 = c(1L, 0L, 1L, 0L, 
0L, 0L,
                0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L)), .Names = c("tv",
                "n1", "n2", "n3"), row.names =
c(NA, -20L), class =
"data.frame")

# Negative log likelihood functions
NLL1 <- function(par, y){
   p1 <- par[1]
   p2 <- par[2]
   mu1 <- par[3]
   mu2 <- par[4]
   n1 <- y$n1
   n2 <- y$n2
   n3 <- y$n3
   t <- y$tv
   P1 <- (p1*((-1 + exp(sqrt((mu1 + mu2 + p1 + p2)^2 -
     4*(mu2*p1 + mu1*(mu2 + p2)))*t))*((-mu2)*(mu2 - p1 + p2) +
     mu1*(mu2 + 2*p2)) - mu2*sqrt((mu1 + mu2 + p1 + p2)^2 -
     4*(mu2*p1 + mu1*(mu2 + p2))) -
     exp(sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))*t)*
     mu2*sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2))) +
     2*exp((1/2)*(mu1 + mu2 + p1 + p2 + sqrt((mu1 + mu2 + p1 + p2)^2 -
     4*(mu2*p1 + mu1*(mu2 + p2))))*t)*mu2*
     sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))))/
     exp((1/2)*(mu1 + mu2 + p1 + p2 + sqrt((mu1 + mu2 + p1 + p2)^2 -
     4*(mu2*p1 + mu1*(mu2 + p2))))*t)/(2*(mu2*p1 + mu1*(mu2 + p2))*
     sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))) 

   P2 <- (p2*((-1 + exp(sqrt((mu1 + mu2 + p1 + p2)^2 -
     4*(mu2*p1 + mu1*(mu2 + p2)))*t))*(-mu1^2 + 2*mu2*p1 +
     mu1*(mu2 - p1 + p2)) - mu1*sqrt((mu1 + mu2 + p1 + p2)^2 -
     4*(mu2*p1 + mu1*(mu2 + p2))) -
     exp(sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))*t)*
     mu1*sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2))) +
     2*exp((1/2)*(mu1 + mu2 + p1 + p2 + sqrt((mu1 + mu2 + p1 + p2)^2 -
     4*(mu2*p1 + mu1*(mu2 + p2))))*t)*mu1*
     sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))))/
     exp((1/2)*(mu1 + mu2 + p1 + p2 + sqrt((mu1 + mu2 + p1 + p2)^2 -
     4*(mu2*p1 + mu1*(mu2 + p2))))*t)/(2*(mu2*p1 + mu1*(mu2 + p2))*
     sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2))))
   P3 <- 1 - P1 - P2
   p.all <- c(P1, P2, P3)
   -sum(dmnom2(c(n1, n2, n3), prob = p.all, log = TRUE))
}

NLL2 <- function(par, y) {
   p1 <- par[1]
   p2 <- p1
   mu1 <- par[2]
   mu2 <- mu1
   n1 <- y$n1
   n2 <- y$n2
   n3 <- y$n3
   t <- y$t
   P1 <- (p1*((-1 + exp(sqrt((mu1 + mu2 + p1 + p2)^2 -
     4*(mu2*p1 + mu1*(mu2 + p2)))*t))*((-mu2)*(mu2 - p1 + p2) +
     mu1*(mu2 + 2*p2)) - mu2*sqrt((mu1 + mu2 + p1 + p2)^2 -
     4*(mu2*p1 + mu1*(mu2 + p2))) -
     exp(sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))*t)*
     mu2*sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2))) +
     2*exp((1/2)*(mu1 + mu2 + p1 + p2 + sqrt((mu1 + mu2 + p1 + p2)^2 -
     4*(mu2*p1 + mu1*(mu2 + p2))))*t)*mu2*
     sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))))/
     exp((1/2)*(mu1 + mu2 + p1 + p2 + sqrt((mu1 + mu2 + p1 + p2)^2 -
     4*(mu2*p1 + mu1*(mu2 + p2))))*t)/(2*(mu2*p1 + mu1*(mu2 + p2))*
     sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))) 

   P2 <- (p2*((-1 + exp(sqrt((mu1 + mu2 + p1 + p2)^2 -
     4*(mu2*p1 + mu1*(mu2 + p2)))*t))*(-mu1^2 + 2*mu2*p1 +
     mu1*(mu2 - p1 + p2)) - mu1*sqrt((mu1 + mu2 + p1 + p2)^2 -
     4*(mu2*p1 + mu1*(mu2 + p2))) -
     exp(sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))*t)*
     mu1*sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2))) +
     2*exp((1/2)*(mu1 + mu2 + p1 + p2 + sqrt((mu1 + mu2 + p1 + p2)^2 -
     4*(mu2*p1 + mu1*(mu2 + p2))))*t)*mu1*
     sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))))/
     exp((1/2)*(mu1 + mu2 + p1 + p2 + sqrt((mu1 + mu2 + p1 + p2)^2 -
     4*(mu2*p1 + mu1*(mu2 + p2))))*t)/(2*(mu2*p1 + mu1*(mu2 + p2))*
     sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2))))
   P3 <- 1 - P1 - P2
   p.all <- c(P1, P2, P3)
   -sum(dmnom2(c(n1, n2, n3), prob = p.all, log = TRUE))
}

# Optimization with spg
spg1 <- spg(par = c(2, 0.3, 0.001, 0.001), fn = NLL1, y = y,
             lower = c(0.0001, 0.0001, 0.0001, 0.0001),
             control = list(maxit = 5000, trace = FALSE))

spg2 <- spg(par = c(2, 0.001), fn = NLL2, y = y,
             lower = c(0.0001, 0.0001),
             control = list(maxit = 5000, trace = FALSE))

# Model selection with ICtab; currently gives an error
ICtab(spg1, spg2)

##########################################################################

Any suggestions would be greatly appreciated.

Thanks in advance,
Adam Zeilinger

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

Adam,

See the attached R code that solves your problem and beyond.  One important
issue is that you are enforcing constraints only indirectly.  You need to make
sure that P1, P2, and P3 (which are functions of original parameters and time)
are all between 0 and 1.  It is not enough to impose constraints on original
parameters p1, p2, mu1 and mu2.

I also show you how to do a likelihood ratio test with the results from spg. 
You can also do the same for nlminb or Rvmmin.

Finally, I also show you how to use optimx to compare different algorithms. This
shows you that in addition to spg, you also get very good results (as seen from
objective function values and KKT conditions) with a number of other algorithms
(e.g., Rvmmin, nlminb), many of which are much faster than spg.

This example illustrates the utility of optimx().  I am always surprised as to
why more R users doing optimization are not using optimx.  This is a very
powerful for benchmarking unconstrained and box-constrained optimization
problems.  It deserves to be used widely, in my biased, but correct, opinion.

Ravi

Ravi Varadhan, Ph.D.
Assistant Professor
The Center on Aging and Health
Division of Geriatric Medicine & Gerontology
Johns Hopkins University
rvaradhan at jhmi.edu<mailto:rvaradhan at jhmi.edu>
410-502-2619

Adam Zeilinger <zeil0006 <at> umn.edu> writes:
> 
> Dear R Help,
> 
> I have two nested negative log-likelihood functions that I am optimizing 
> with the spg function [BB package].  I would like to perform model 
> selection on these two objective functions using AIC (and possibly 
> anova() too).  However, the spg() function returns a list and I need a 
> fitted model object for AIC(), ICtab() [bbmle package], or anova().
> 
> How can I perform AIC-based model selection on two spg-optimized 
> objective functions?  Alternatively, how can I convert the list returned 
> by spg into a fitted model object that can be run in AIC, ICtab, or anova?
> 
  I believe you can use spg within mle2 by specifying
optimizer="optimx", method="spg" ... that should give you a
fitted
mle2 object that you can work with.

  Alternately, it's just not that hard to compute the AICs yourself ...
or a likelihood ratio test ...

  Ben Bolker