Dear All,
I am trying to find out how R meta packages has been calculating ratio of means from the Metacont package. It will be helpful to know what formula R is using. If anyone knows your help will be very much appreciated.
Sincerely,
Joe
Are you referencing the metacont() function from, the package meta?
https://www.rdocumentation.org/packages/meta/versions/4.9-9/topics/metacont
Having trouble finding an R package "metacont" via Google...
Thanks Nate for your response
I wanted to find out the exact mathematical formula it uses. We are trying to see whether it is standardised using the Control rather than Mean(intervention)/Mean Control. The documentation is still not clear about the mathematics that was used in R.
Regards,
Joe
If I understand and you want to look at the source code of metacont() you can get there two ways from an R console:
-
getAnywhere(metacont)this is a useful strategy for getting to any function's source code (eg even if it is a hidden function that you could not call directly) -
meta::metacontyou would use three colons:::if the function was hidden (not-exported)
> getAnywhere(metacont)
A single object matching ‘metacont’ was found
It was found in the following places
namespace:meta
with value
function (n.e, mean.e, sd.e, n.c, mean.c, sd.c, studlab, data = NULL,
subset = NULL, exclude = NULL, sm = gs("smcont"), pooledvar = gs("pooledvar"),
method.smd = gs("method.smd"), sd.glass = gs("sd.glass"),
exact.smd = gs("exact.smd"), level = gs("level"), level.comb = gs("level.comb"),
comb.fixed = gs("comb.fixed"), comb.random = gs("comb.random"),
overall = comb.fixed | comb.random, overall.hetstat = comb.fixed |
comb.random, hakn = gs("hakn"), method.tau = gs("method.tau"),
method.tau.ci = if (method.tau == "DL") "J" else "QP", tau.preset = NULL,
TE.tau = NULL, tau.common = gs("tau.common"), prediction = gs("prediction"),
level.predict = gs("level.predict"), method.bias = gs("method.bias"),
backtransf = gs("backtransf"), title = gs("title"), complab = gs("complab"),
outclab = "", label.e = gs("label.e"), label.c = gs("label.c"),
label.left = gs("label.left"), label.right = gs("label.right"),
byvar, bylab, print.byvar = gs("print.byvar"), byseparator = gs("byseparator"),
keepdata = gs("keepdata"), warn = gs("warn"), control = NULL)
{
chknull(sm)
chklevel(level)
chklevel(level.comb)
chklogical(comb.fixed)
chklogical(comb.random)
chklogical(overall)
chklogical(overall.hetstat)
chklogical(hakn)
method.tau <- setchar(method.tau, .settings$meth4tau)
method.tau.ci <- setchar(method.tau.ci, .settings$meth4tau.ci)
chklogical(tau.common)
chklogical(prediction)
chklevel(level.predict)
method.bias <- setchar(method.bias, c("rank", "linreg", "mm",
"count", "score", "peters"))
chklogical(keepdata)
fun <- "metacont"
sm <- setchar(sm, .settings$sm4cont)
chklogical(pooledvar)
method.smd <- setchar(method.smd, c("Hedges", "Cohen", "Glass"))
sd.glass <- setchar(sd.glass, c("control", "experimental"))
chklogical(warn)
nulldata <- is.null(data)
if (nulldata)
data <- sys.frame(sys.parent())
mf <- match.call()
n.e <- eval(mf[[match("n.e", names(mf))]], data, enclos = sys.frame(sys.parent()))
chknull(n.e)
k.All <- length(n.e)
mean.e <- eval(mf[[match("mean.e", names(mf))]], data, enclos = sys.frame(sys.parent()))
chknull(mean.e)
sd.e <- eval(mf[[match("sd.e", names(mf))]], data, enclos = sys.frame(sys.parent()))
chknull(sd.e)
n.c <- eval(mf[[match("n.c", names(mf))]], data, enclos = sys.frame(sys.parent()))
chknull(n.c)
mean.c <- eval(mf[[match("mean.c", names(mf))]], data, enclos = sys.frame(sys.parent()))
chknull(mean.c)
sd.c <- eval(mf[[match("sd.c", names(mf))]], data, enclos = sys.frame(sys.parent()))
chknull(sd.c)
studlab <- eval(mf[[match("studlab", names(mf))]], data,
enclos = sys.frame(sys.parent()))
studlab <- setstudlab(studlab, k.All)
byvar <- eval(mf[[match("byvar", names(mf))]], data, enclos = sys.frame(sys.parent()))
by <- !is.null(byvar)
subset <- eval(mf[[match("subset", names(mf))]], data, enclos = sys.frame(sys.parent()))
missing.subset <- is.null(subset)
exclude <- eval(mf[[match("exclude", names(mf))]], data,
enclos = sys.frame(sys.parent()))
missing.exclude <- is.null(exclude)
chklength(mean.e, k.All, fun)
chklength(sd.e, k.All, fun)
chklength(n.c, k.All, fun)
chklength(mean.c, k.All, fun)
chklength(sd.c, k.All, fun)
chklength(studlab, k.All, fun)
if (by)
chklength(byvar, k.All, fun)
if (!by & tau.common) {
warning("Value for argument 'tau.common' set to FALSE as argument 'byvar' is missing.")
tau.common <- FALSE
}
if (by & !tau.common & !is.null(tau.preset)) {
warning("Argument 'tau.common' set to TRUE as argument tau.preset is not NULL.")
tau.common <- TRUE
}
if (!missing.subset)
if ((is.logical(subset) & (sum(subset) > k.All)) || (length(subset) >
k.All))
stop("Length of subset is larger than number of studies.")
if (!missing.exclude) {
if ((is.logical(exclude) & (sum(exclude) > k.All)) ||
(length(exclude) > k.All))
stop("Length of argument 'exclude' is larger than number of studies.")
exclude2 <- rep(FALSE, k.All)
exclude2[exclude] <- TRUE
exclude <- exclude2
}
else exclude <- rep(FALSE, k.All)
if (by) {
chkmiss(byvar)
byvar.name <- byvarname(mf[[match("byvar", names(mf))]])
bylab <- if (!missing(bylab) && !is.null(bylab))
bylab
else byvar.name
}
if (keepdata) {
if (nulldata)
data <- data.frame(.n.e = n.e)
else data$.n.e <- n.e
data$.mean.e <- mean.e
data$.sd.e <- sd.e
data$.n.c <- n.c
data$.mean.c <- mean.c
data$.sd.c <- sd.c
data$.studlab <- studlab
if (by)
data$.byvar <- byvar
if (!missing.subset) {
if (length(subset) == dim(data)[1])
data$.subset <- subset
else {
data$.subset <- FALSE
data$.subset[subset] <- TRUE
}
}
if (!missing.exclude)
data$.exclude <- exclude
}
if (!missing.subset) {
n.e <- n.e[subset]
mean.e <- mean.e[subset]
sd.e <- sd.e[subset]
n.c <- n.c[subset]
mean.c <- mean.c[subset]
sd.c <- sd.c[subset]
studlab <- studlab[subset]
exclude <- exclude[subset]
if (by)
byvar <- byvar[subset]
}
k.all <- length(n.e)
if (k.all == 0)
stop("No studies to combine in meta-analysis.")
if (k.all == 1) {
comb.fixed <- FALSE
comb.random <- FALSE
prediction <- FALSE
overall <- FALSE
overall.hetstat <- FALSE
}
chknumeric(n.e)
chknumeric(mean.e)
chknumeric(sd.e)
chknumeric(n.c)
chknumeric(mean.c)
chknumeric(sd.c)
n.e <- int2num(n.e)
mean.e <- int2num(mean.e)
sd.e <- int2num(sd.e)
n.c <- int2num(n.c)
mean.c <- int2num(mean.c)
sd.c <- int2num(sd.c)
npn.n <- npn(n.e) | npn(n.c)
N <- n.e + n.c
if (sm == "MD" | sm == "ROM")
var.pooled <- ((n.e - 1) * sd.e^2 + (n.c - 1) * sd.c^2)/(N -
2)
if (any(npn.n) & warn)
warning("Note, studies with non-positive values for n.e and / or n.c get no weight in meta-analysis.")
if (sm == "MD") {
TE <- ifelse(npn.n, NA, mean.e - mean.c)
if (pooledvar)
seTE <- ifelse(npn.n, NA, sqrt(var.pooled * (1/n.e +
1/n.c)))
else seTE <- ifelse(npn.n, NA, sqrt(sd.e^2/n.e + sd.c^2/n.c))
seTE[is.na(TE)] <- NA
}
else if (sm == "SMD") {
J <- function(x) gamma(x/2)/(sqrt(x/2) * gamma((x - 1)/2))
K <- function(x) 1 - (x - 2)/(x * J(x)^2)
if (method.smd %in% c("Hedges", "Cohen"))
S.within <- sqrt(((n.e - 1) * sd.e^2 + (n.c - 1) *
sd.c^2)/(N - 2))
else S.within <- if (sd.glass == "control")
sd.c
else sd.e
smd <- ifelse(npn.n, NA, (mean.e - mean.c)/S.within)
if (method.smd == "Cohen") {
TE <- smd
if (exact.smd) {
J <- function(x) gamma(x/2)/(sqrt(x/2) * gamma((x -
1)/2))
K <- function(x) 1 - (x - 2)/(x * J(x)^2)
seTE <- ifelse(npn.n, NA, sqrt(N/(n.e * n.c) +
(J(N - 2) * smd)^2 * K(N - 2)))
}
else seTE <- ifelse(npn.n, NA, sqrt(N/(n.e * n.c) +
TE^2/(2 * N)))
}
else if (method.smd == "Hedges") {
if (exact.smd) {
J <- function(x) gamma(x/2)/(sqrt(x/2) * gamma((x -
1)/2))
K <- function(x) 1 - (x - 2)/(x * J(x)^2)
}
else {
J <- function(x) 1 - 3/(4 * x - 1)
K <- function(x) 1/(2 * (x - 1.94))
}
TE <- J(N - 2) * smd
seTE <- ifelse(npn.n, NA, sqrt(N/(n.e * n.c) + TE^2 *
K(N - 2)))
}
else if (method.smd == "Glass") {
n.g <- if (sd.glass == "control")
n.c
else n.e
TE <- smd
seTE <- ifelse(npn.n, NA, sqrt(N/(n.e * n.c) + TE^2/(2 *
n.g - 1)))
}
seTE[is.na(TE)] <- NA
}
else if (sm == "ROM") {
npn.mean <- npn(mean.e) | npn(mean.c)
if (any(npn.mean) & warn)
warning("Note, studies with negative or zero means get no weight in meta-analysis.")
TE <- ifelse(npn.n | npn.mean, NA, log(mean.e/mean.c))
if (pooledvar)
seTE <- ifelse(npn.n, NA, sqrt(var.pooled * (1/(n.e *
mean.e^2) + 1/(n.c * mean.c^2))))
else seTE <- ifelse(npn.n | npn.mean, NA, sqrt(sd.e^2/(n.e *
mean.e^2) + sd.c^2/(n.c * mean.c^2)))
seTE[is.na(TE)] <- NA
}
sel <- sd.e <= 0 | sd.c <= 0
if (any(sel, na.rm = TRUE) & warn)
warning("Note, studies with non-positive values for sd.e or sd.c get no weight in meta-analysis.")
seTE[sel] <- NA
if (sm == "SMD")
TE[sel] <- NA
m <- metagen(TE, seTE, studlab, exclude = if (missing.exclude)
NULL
else exclude, sm = sm, level = level, level.comb = level.comb,
comb.fixed = comb.fixed, comb.random = comb.random, overall = overall,
overall.hetstat = overall.hetstat, hakn = hakn, method.tau = method.tau,
method.tau.ci = method.tau.ci, tau.preset = tau.preset,
TE.tau = TE.tau, tau.common = FALSE, prediction = prediction,
level.predict = level.predict, method.bias = method.bias,
backtransf = backtransf, title = title, complab = complab,
outclab = outclab, label.e = label.e, label.c = label.c,
label.left = label.left, label.right = label.right, keepdata = FALSE,
warn = warn, control = control)
if (by & tau.common) {
hcc <- hetcalc(TE, seTE, method.tau, "", TE.tau, level.comb,
byvar, control)
}
res <- list(n.e = n.e, mean.e = mean.e, sd.e = sd.e, n.c = n.c,
mean.c = mean.c, sd.c = sd.c, pooledvar = pooledvar,
method.smd = method.smd, sd.glass = sd.glass, exact.smd = exact.smd)
m$n.e <- NULL
m$n.c <- NULL
res <- c(res, m)
res$call <- match.call()
if (keepdata) {
res$data <- data
if (!missing.subset)
res$subset <- subset
}
class(res) <- c(fun, "meta")
if (by) {
res$byvar <- byvar
res$bylab <- bylab
res$print.byvar <- print.byvar
res$byseparator <- byseparator
res$tau.common <- tau.common
if (!tau.common) {
res <- c(res, subgroup(res))
res$tau2.resid <- res$lower.tau2.resid <- res$upper.tau2.resid <- NA
res$tau.resid <- res$lower.tau.resid <- res$upper.tau.resid <- NA
}
else if (!is.null(tau.preset)) {
res <- c(res, subgroup(res, tau.preset))
res$tau2.resid <- res$lower.tau2.resid <- res$upper.tau2.resid <- NA
res$tau.resid <- res$lower.tau.resid <- res$upper.tau.resid <- NA
}
else {
res <- c(res, subgroup(res, hcc$tau))
res$Q.w.random <- hcc$Q
res$df.Q.w.random <- hcc$df.Q
res$tau2.resid <- hcc$tau2
res$lower.tau2.resid <- hcc$lower.tau2
res$upper.tau2.resid <- hcc$upper.tau2
res$tau.resid <- hcc$tau
res$lower.tau.resid <- hcc$lower.tau
res$upper.tau.resid <- hcc$upper.tau
res$sign.lower.tau.resid <- hcc$sign.lower.tau
res$sign.upper.tau.resid <- hcc$sign.upper.tau
}
if (!tau.common || method.tau == "DL") {
ci.H.resid <- calcH(res$Q.w.fixed, res$df.Q.w, level.comb)
res$H.resid <- ci.H.resid$TE
res$lower.H.resid <- ci.H.resid$lower
res$upper.H.resid <- ci.H.resid$upper
}
else {
res$H.resid <- hcc$H.resid
res$lower.H.resid <- hcc$lower.H.resid
res$upper.H.resid <- hcc$upper.H.resid
}
if (!tau.common || method.tau == "DL") {
ci.I2.resid <- isquared(res$Q.w.fixed, res$df.Q.w,
level.comb)
res$I2.resid <- ci.I2.resid$TE
res$lower.I2.resid <- ci.I2.resid$lower
res$upper.I2.resid <- ci.I2.resid$upper
}
else {
res$I2.resid <- hcc$I2.resid
res$lower.I2.resid <- hcc$lower.I2.resid
res$upper.I2.resid <- hcc$upper.I2.resid
}
res$event.e.w <- NULL
res$event.c.w <- NULL
res$event.w <- NULL
res$n.w <- NULL
res$time.e.w <- NULL
res$time.c.w <- NULL
}
class(res) <- c(fun, "meta")
res
}
<bytecode: 0x7feff144e3d0>
<environment: namespace:meta>
Thank you again Nate
@Nate
Is there a way to Amend the internal metacont functions and code, so that i can incorporate the standardisation I need?
Regards
Joe
I think there are three reasonable options forward, ranked by my preference.
- See if you can build a new variable
scaled = FALSE(for the default to keep with current behavior) intometaconf()that does what you want. If that succeeds I would open an issue/pull-request on the repository to see if the maintainer wants to add your scaling functionality. (I bet you aren't the only person that would benefit from scaling) - Build a function to be used right before calling
metaconf()that handles the variable scaling for you, similar to?scale() - Build a new function by copying the code of the
metacont()into a new onemetacont_scaled()and then in the new one you make the scaled changes. This is very similar to #1 and it would be your function so you'd be the one making updates to it as the mainmetaconf()changes.
This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.