1

In my R code below, suppose I want to compare all unique 2 m objects using a similar R routine. For example, to compare m1 and m2, my routine is:

pchisq(2 * (logLik(m2) - logLik(m1)), df = abs(m1$df.residual - m2$df.residual), lower = F)

Question:

I was wondering how I could make a function to make all unique pairwise comparisons for all m objects using my routine?

Here is what I've tried without success:

## Suppose we have 4 `m` objects: `m1...m4` (6 unique pairwise comparisons possible)

m1 <- lm(hp ~ vs, mtcars)
m2 <- lm(hp ~ vs*wt, mtcars)
m3 <- lm(hp ~., mtcars)
m4 <- lm(hp ~ vs * gear * wt, mtcars)


 compare <- function(...){

   m <- list(...)
   L <- length(m) - 1

   lapply(1:L, function(i) pchisq(2 * (logLik(m[[i+1]]) - logLik(m[[i]])), df = abs(m[[i]]$df.residual - m[[i+1]]$df.residual), lower.tail = FALSE) )

  }
# Example of use:
compare(m1, m2, m3, m4)
2
  • 1
    Not what you asked for but possibly what you really need is: anova(m1, m2, m4, m3, test = "Chisq") Feb 5, 2019 at 22:02
  • 1
    The class of the models shown in the question is lm and anova supports that and there are many other methods as well. If you are using something not shown suggest you mention it. Feb 5, 2019 at 22:14

1 Answer 1

3

You can use combn to get all combinations you want to compare:

compare <- function(...){

  m <- list(...)
  n_mod <- length(m)
  names(m) <- sapply(substitute(list(...))[-1], deparse)
  combs <- t(combn(x = names(m), m = 2))

  comp_value <- apply(X = combs, MARGIN = 1, function(ind) pchisq(2 * (logLik(m[[ind[2]]]) - logLik(m[[ind[1]]])), df = abs(m[[ind[1]]]$df.residual - m[[ind[2]]]$df.residual), lower.tail = FALSE))
  df_out <- data.frame(combs, comp_value)
  names(df_out) <- c("mod_1", "mod_2", "comp_value")

  return(df_out)
}

So, to make it easier to read the result, you can return a data.frame with all the comparisons.

Then

compare(m1, m2, m3, m4)
  mod_1 mod_2   comp_value
1    m1    m2 2.391012e-02
2    m1    m3 7.253068e-08
3    m1    m4 1.248692e-06
4    m2    m3 2.735901e-07
5    m2    m4 4.256098e-06
6    m3    m4 1.000000e+00
4
  • Douglas, thank you once again. But I have realized that your function throws an error when I compare 2 fitted models from the glmmTMB package, like so: m1 <- glmmTMB(count~ mined + (1|site), zi=~mined, family=poisson, data=Salamanders) and m2 <- glmmTMB(count~spp + mined + (1|site), zi=~spp + mined, family=nbinom2, Salamanders). And now when I do compare(m1, m2), I get the following error: Error in -object$fit$objective : invalid argument to unary operator. The problem is with your third line (containing sapply)?
    – rnorouzian
    Feb 22, 2019 at 5:24
  • The problem is that the glmmTMB class does not have this entry: m1$df.residuals. Feb 23, 2019 at 22:03
  • Douglas, thank you. hummm, could you possibly think of a fix?
    – rnorouzian
    Feb 24, 2019 at 0:39
  • In this class you can get the degrees of freedom using glmmTMB:::df.residual.glmmTMB(mod). So you need to use this code instead of mod$df.residuals in the function. In addition, you can add a step in the function: if(class(mod) == "glmmTMB") df <- glmmTMB:::df.residual.glmmTMB(mod) otherwise use the current function. Feb 24, 2019 at 16:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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