Looping all pairwise comparisons from a list in R

Question

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)

Not what you asked for but possibly what you really need is: `anova(m1, m2, m4, m3, test = "Chisq")` — G. Grothendieck, Feb 05 '19 at 22:02
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. — G. Grothendieck, Feb 05 '19 at 22:14

score 3 · Accepted Answer · answered Feb 05 '19 at 22:15

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

answered Feb 05 '19 at 22:15

Douglas Mesquita

1,011
7
12

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 '19 at 05:24
The problem is that the `glmmTMB` class does not have this entry: `m1$df.residuals`. – Douglas Mesquita Feb 23 '19 at 22:03
Douglas, thank you. hummm, could you possibly think of a fix? – rnorouzian Feb 24 '19 at 00: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. – Douglas Mesquita Feb 24 '19 at 16:16

Looping all pairwise comparisons from a list in R

1 Answers1