Nothing
nma.ruecker <- function(TE, seTE,
treat1, treat2,
treat1.pos, treat2.pos,
narms, studlab,
sm = "",
level, level.ma,
seTE.orig, tau.direct = 0, sep.trts = ":",
method.tau = "DL",
func.inverse) {
w.pooled <- 1 / seTE^2
m <- length(TE) # Number of pairwise comparisons (edges)
n <- length(unique(c(treat1, treat2))) # Number of treatments (vertices)
W <- diag(w.pooled, # Weighted degree diagonal matrix
nrow = length(w.pooled)) #
df1 <- 2 * sum(1 / narms) # Sum of degrees of freedom per study
##
## B is the edge-vertex incidence matrix (m x n)
##
B <- createB(treat1.pos, treat2.pos, ncol = n)
##
## B.full is the full edge-vertex incidence matrix (m x n)
##
B.full <- createB(ncol = n)
##
## M is the unweighted Laplacian, D its diagonal,
## A is the adjacency matrix
##
M <- t(B) %*% B # unweighted Laplacian matrix
D <- diag(diag(M)) # diagonal matrix
A <- D - M # adjacency matrix (n x n)
##
## L is the weighted Laplacian (Kirchhoff) matrix (n x n)
## Lplus is its Moore-Penrose pseudoinverse
##
L <- t(B) %*% W %*% B
Lplus <- do.call(func.inverse, list(X = L))
Lplus[is.zero(Lplus)] <- 0
##
## R resistance distance (variance) matrix (n x n)
##
R <- matrix(0, nrow = n, ncol = n)
for (i in 1:n) {
for (j in 1:n) {
R[i, j] <- Lplus[i, i] + Lplus[j, j] - 2 * Lplus[i, j]
}
}
##
## V is the vector of effective variances
##
V <- vector(length = m, mode = "numeric")
for (i in 1:m) {
V[i] <- R[treat1.pos[i], treat2.pos[i]]
}
##
## G is the matrix B %*% Lplus %*% t(B)
## H is the projection matrix (also called "hat matrix")
##
## Interpretation:
## (i) diag(G) = V The effective variances
## (ii) diag(H) = V %*% W = V * w The leverages
## (iii) sum(diag(H)) = n - 1 Rank of projection
## (iv) mean(diag(H)) = (n - 1) / m Mean leverage = average efficiency
##
G <- B %*% Lplus %*% t(B)
H <- G %*% W
##
## Variance-covariance matrix for all comparisons
##
Cov <- B.full %*% Lplus %*% t(B.full)
##
## Resulting effects and variances at numbered edges
##
v <- as.vector(H %*% TE)
ci.v <- ci(v, sqrt(V), level = level)
##
## Resulting effects, all edges, as a n x n matrix:
##
all <- matrix(NA, nrow = n, ncol = n)
##
for (i in 1:m) {
all[treat1.pos[i], treat2.pos[i]] <- v[i]
}
##
for (i in 1:n) {
for (j in 1:n) {
for (k in 1:n) {
if (!is.na(all[i, k]) & !is.na(all[j, k])) {
all[i, j] <- all[i, k] - all[j, k]
all[j, i] <- all[j, k] - all[i, k]
}
if (!is.na(all[i, j]) & !is.na(all[k, j])) {
all[i, k] <- all[i, j] - all[k, j]
all[k, i] <- all[k, j] - all[i, j]
}
if (!is.na(all[i, k]) & !is.na(all[i, j])) {
all[j, k] <- all[i, k] - all[i, j]
all[k, j] <- all[i, j] - all[i, k]
}
}
}
}
##
## Test of total heterogeneity / inconsistency:
##
Q <- as.vector(t(TE - v) %*% W %*% (TE - v))
df <- df1 - (n - 1)
if (df == 0)
pval.Q <- NA
else
pval.Q <- pchisq(Q, df, lower.tail = FALSE)
##
## Heterogeneity variance
##
I <- diag(m)
E <- matrix(0, nrow = m, ncol = m)
for (i in 1:m)
for (j in 1:m)
E[i, j] <- as.numeric(studlab[i] == studlab[j])
##
if (df == 0) {
tau2 <- NA
tau <- NA
I2 <- lower.I2 <- upper.I2 <- NA
}
else {
tau2 <-
max(0, (Q - df) / sum(diag((I - H) %*% (B %*% t(B) * E / 2) %*% W)))
tau <- sqrt(tau2)
ci.I2 <- isquared(Q, df, level.ma)
I2 <- ci.I2$TE
lower.I2 <- ci.I2$lower
upper.I2 <- ci.I2$upper
}
##
## Decomposition of total Q into parts from pairwise meta-analyses
## and residual inconsistency
##
Q.matrix <- matrix(0, nrow = n, ncol = n)
n.pairwise <- 0
data <- data.frame(treat1.pos, treat2.pos,
TE, w.pooled,
stringsAsFactors = FALSE)
##
for (i in 1:(n - 1)) {
for (j in (i + 1):n) {
if (A[i, j] > 1) {
n.pairwise <- n.pairwise + 1
sub <- data[(data$treat1.pos == i & data$treat2.pos == j) |
(data$treat1.pos == j & data$treat2.pos == i), ]
l <- nrow(sub)
Q.matrix[i, j] <-
suppressWarnings(metagen(TE = sub$TE, seTE = 1 / sqrt(sub$w.pooled),
method.tau = "DL", method.tau.ci = "")$Q)
}
}
}
q <- t1 <- t2 <- dfs <- vector(length = n.pairwise, mode = "numeric")
p <- 0
for (i in 1:(n - 1)) {
for (j in (i + 1):n) {
if (A[i, j] > 0) {
p <- p + 1
t1[p] <- i
t2[p] <- j
q[p] <- Q.matrix[i, j]
dfs[p] <- A[i, j] - 1
}
}
}
##
Q.heterogeneity <- sum(Q.matrix)
Q.inconsistency <- Q - Q.heterogeneity
##
names.treat <- sort(unique(c(treat1, treat2)))
##
Q.decomp <- data.frame(treat1 = names.treat[t1],
treat2 = names.treat[t2],
Q = q,
df = dfs,
pval.Q = pchisq(q, dfs, lower.tail = FALSE),
stringsAsFactors = FALSE)
##
Q.decomp$pval.Q[Q.decomp$df == 0] <- NA
TE.pooled <- all
seTE.pooled <- sqrt(R)
##
ci.pooled <- meta::ci(all, sqrt(R), level = level.ma)
##
lower.pooled <- ci.pooled$lower
upper.pooled <- ci.pooled$upper
statistic.pooled <- ci.pooled$statistic
pval.pooled <- ci.pooled$p
##
rownames(TE.pooled) <- colnames(TE.pooled) <- names.treat
rownames(seTE.pooled) <- colnames(seTE.pooled) <- names.treat
rownames(lower.pooled) <- colnames(lower.pooled) <- names.treat
rownames(upper.pooled) <- colnames(upper.pooled) <- names.treat
rownames(statistic.pooled) <- colnames(statistic.pooled) <- names.treat
rownames(pval.pooled) <- colnames(pval.pooled) <- names.treat
A.matrix <- A
L.matrix <- L
Lplus.matrix <- Lplus
##
rownames(A.matrix) <- colnames(A.matrix) <- names.treat
rownames(L.matrix) <- colnames(L.matrix) <- names.treat
rownames(Lplus.matrix) <- colnames(Lplus.matrix) <- names.treat
rownames(Q.matrix) <- colnames(Q.matrix) <- names.treat
G.matrix <- G
H.matrix <- H
##
rownames(G.matrix) <- colnames(G.matrix) <- studlab
rownames(H.matrix) <- colnames(H.matrix) <- studlab
names.Cov <- rep("", nrow(B.full))
##
k <- 0
for (i in 1:(n - 1)) {
for (j in (i + 1):n) {
k <- k + 1
names.Cov[k] <- paste(names.treat[i], names.treat[j], sep = sep.trts)
}
}
##
rownames(Cov) <- colnames(Cov) <- names.Cov
## Calculate direct treatment estimates
##
Bo <- B
Bo[Bo == 1] <- -2
rownames(B) <- studlab
colnames(B) <- names.treat
##
B.matrix <- B[do.call(order, data.frame(Bo)), , drop = FALSE]
##
TE.direct <- seTE.direct <-
Q.direct <- tau2.direct <- I2.direct <-
matrix(NA, ncol = length(names.treat), nrow = length(names.treat))
##
rownames(TE.direct) <- colnames(TE.direct) <-
rownames(seTE.direct) <- colnames(seTE.direct) <-
rownames(Q.direct) <- colnames(Q.direct) <-
rownames(tau2.direct) <- colnames(tau2.direct) <-
rownames(I2.direct) <- colnames(I2.direct) <-
names.treat
##
C.matrix <- B.matrix[!duplicated(B.matrix), , drop = FALSE]
##
for (i in 1:dim(C.matrix)[1]) {
sel1 <- C.matrix[i, ] == 1
sel2 <- C.matrix[i, ] == -1
sel.treat1 <- colnames(C.matrix)[sel1]
sel.treat2 <- colnames(C.matrix)[sel2]
##
selstud <- treat1 == sel.treat1 & treat2 == sel.treat2
##
m.i.tau.preset <-
suppressWarnings(metagen(TE, seTE.orig, subset = selstud,
tau.preset = tau.direct,
method.tau = "DL", method.tau.ci = ""))
##
if (is.na(tau.direct) | tau.direct == 0) {
TE.i <- m.i.tau.preset$TE.common
seTE.i <- m.i.tau.preset$seTE.common
}
else {
TE.i <- m.i.tau.preset$TE.random
seTE.i <- m.i.tau.preset$seTE.random
}
##
TE.direct[sel.treat1, sel.treat2] <- TE.i
seTE.direct[sel.treat1, sel.treat2] <- seTE.i
##
TE.direct[sel.treat2, sel.treat1] <- -TE.i
seTE.direct[sel.treat2, sel.treat1] <- seTE.i
##
m.i <-
suppressWarnings(metagen(TE, seTE.orig, subset = selstud,
method.tau = method.tau, method.tau.ci = ""))
##
Q.direct[sel.treat1, sel.treat2] <-
Q.direct[sel.treat2, sel.treat1] <- m.i$Q
tau2.direct[sel.treat1, sel.treat2] <-
tau2.direct[sel.treat2, sel.treat1] <- m.i$tau2
I2.direct[sel.treat1, sel.treat2] <-
I2.direct[sel.treat2, sel.treat1] <- m.i$I2
}
##
ci.direct <- meta::ci(TE.direct, seTE.direct, level = level.ma)
##
lower.direct <- ci.direct$lower
upper.direct <- ci.direct$upper
statistic.direct <- ci.direct$statistic
pval.direct <- ci.direct$p
##
## Contribution of individual studies to Q
##
Q.pooled <- w.pooled * (TE - v)^2
res <- list(studlab = studlab,
treat1 = treat1, treat2 = treat2,
TE = TE, seTE = seTE,
seTE.orig = seTE.orig,
TE.nma = ci.v$TE,
seTE.nma = ci.v$seTE,
lower.nma = ci.v$lower,
upper.nma = ci.v$upper,
statistic.nma = ci.v$statistic,
pval.nma = ci.v$p,
leverage = diag(H),
w.pooled = w.pooled,
Q.pooled = Q.pooled,
treat1.pos = treat1.pos,
treat2.pos = treat2.pos,
##
TE.pooled = TE.pooled,
seTE.pooled = seTE.pooled,
lower.pooled = lower.pooled,
upper.pooled = upper.pooled,
statistic.pooled = statistic.pooled,
pval.pooled = pval.pooled,
##
k = length(unique(studlab)),
m = length(TE),
n = dim(TE.pooled)[[1]],
Q = Q,
df = df,
pval.Q = pval.Q,
I2 = I2,
lower.I2 = lower.I2,
upper.I2 = upper.I2,
tau = tau,
Q.heterogeneity = Q.heterogeneity,
Q.inconsistency = Q.inconsistency,
##
sm = sm,
level = level,
level.ma = level.ma,
##
A.matrix = A.matrix,
B.matrix = B,
L.matrix = L.matrix,
Lplus.matrix = Lplus.matrix,
Q.matrix = Q.matrix,
G.matrix = G.matrix,
H.matrix = H.matrix,
##
Cov = Cov,
##
TE.direct = TE.direct,
seTE.direct = seTE.direct,
lower.direct = lower.direct,
upper.direct = upper.direct,
statistic.direct = statistic.direct,
pval.direct = pval.direct,
##
Q.direct = Q.direct,
tau2.direct = tau2.direct,
I2.direct = I2.direct,
##
Q.decomp = Q.decomp
)
res$version <- packageDescription("netmeta")$Version
res
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.