incorrect ns from metafor objects

[acoppock]

Thanks as always for this fabulous package. I'm encountering an issue with ns from metafor objects, can you point me in the right direction?

library(modelsummary)
library(metafor)

dat <-
  data.frame(
    estimate = rnorm(100),
    std.error = 0.4,
    X = rbinom(100, 1, 0.5)
  )

fit_1 <- rma(estimate ~ 1, sei = std.error, data = dat)
fit_2 <- rma(estimate ~ X, sei = std.error, data = dat)

# gives the wrong ns
modelsummary(models = list(fit_1, fit_2))

# but these are correct
glance(fit_1)$nobs
glance(fit_2)$nobs

[vincentarelbundock]

Thanks for the report. The number of observations reported is the one extracted by the nobs() method supplied by the metafor package itself.

I paste the code of the metafor-provided function below.

Do you (or perhaps @bwiernik ) know why metafor would return these figures as nobs?

library(metafor)
dat <- data.frame( 
  estimate = rnorm(100), std.error = 0.4, X = rbinom(100, 1, 0.5) 
)

fit_1 <- rma(estimate ~ 1, sei = std.error, data = dat) 
fit_2 <- rma(estimate ~ X, sei = std.error, data = dat)

stats::nobs(fit_1)
    [1] 99

stats::nobs(fit_2)

    [1] 98

Function:

> metafor:::nobs.rma
function (object, ...) 
{
    mstyle <- .get.mstyle()
    .chkclass(class(object), must = "rma")
    n.obs <- object$k.eff - ifelse(object$method == "REML", 1, 
        0) * object$p.eff
    return(n.obs)
}
<bytecode: 0x000001b79fca9c80>
<environment: namespace:metafor>

[bwiernik]

It's returning the effective n that's used for computing log-likelihood/REML criteria and related stats (eg, AIC, BIC). It seems a bit of an odd choice to me, not sure what @wviechtb's reasoning was there.

You can get the raw number of rows of data with

model$k.all

[wviechtb]

It's been a while since I added this method (in version 1.7-0, which was released in 2013). I am vaguely recalling that this had something to do with package interoperability, where some package was expecting nobs() to return the effective n (where, for $p$ fixed effects, this is $n-p$ for REML estimation). It may also have been motivated by what help(nobs) says, namely that this is "principally intended to be used in computing BIC", where one would also use $n-p$ for the sample size for REML estimation.

Note that model$k.all is not what you want, since this is the number of studies before subsetting. It should be model$k.

In any case, I will consider changing what nobs() returns.

[bwiernik]

Thanks Wolfgang!

[vincentarelbundock]

Thanks a lot for your input @bwiernik. And thanks @wviechtb for looking into this. I really appreciate both of your help!

@acoppock, if/when this is changed upstream in metafor, the change should automatically be reflected in modelsummary tables, regardless of the version you use. I really think this is best fixed upstream, rather than as a hackish workaround here, so I’ll close the issue for now.

In the meantime, you can easily define a glance_custom.rma() method to get the info you want in the tables:

library(metafor)
library(modelsummary)

dat <- data.frame( 
  estimate = rnorm(100), std.error = 0.4, X = rbinom(100, 1, 0.5) 
)

mod <- list(
  rma(estimate ~ 1, sei = std.error, data = dat),
  rma(estimate ~ X, sei = std.error, data = dat)
)

glance_custom.rma <- \(x) data.frame("nobs" = x$k)

modelsummary(mod, gof_omit = "i2")

+-----------+---------+---------+
|           | (1)     | (2)     |
+===========+=========+=========+
| overall   | 0.168   |         |
+-----------+---------+---------+
|           | (0.109) |         |
+-----------+---------+---------+
| intercept |         | 0.266   |
+-----------+---------+---------+
|           |         | (0.147) |
+-----------+---------+---------+
| X         |         | -0.218  |
+-----------+---------+---------+
|           |         | (0.218) |
+-----------+---------+---------+
| Num.Obs.  | 100     | 100     |
+-----------+---------+---------+
| AIC       | 303.4   | 304.4   |
+-----------+---------+---------+
| BIC       | 308.6   | 312.2   |
+-----------+---------+---------+ 

[acoppock]

Great, I'll use that custom glance method, TY!

[wviechtb]

Just to let you know, I have adjusted nobs() in metafor to return the number of studies. This is now in the development version (https://wviechtb.github.io/metafor/#installation) which will go to CRAN eventually.

[vincentarelbundock]

Thanks a lot, this is awesome!