# R/covariance.R In dosresmeta: Multivariate Dose-Response Meta-Analysis

```#' Computes mean and standardized mean differences for continuous outcome with corresponding
#' co(variance) matrix
#'
#' @description This internal function computes mean and standardized mean of a continuous outcome with the corresponding
#' variances. It also reconstructs the covariance matrix from the available data.
#'
#' @param y a vector defining the mean outcome for each treatment level.
#' @param sd a vector defining the standard deviation of the outcome for each treatment level.
#' @param n a vector defining the number of subjects for each treatment level.
#' @param measure character string, indicating the measure to be calculated. Options are \code{md}
#' and \code{smd} for mean difference and standardized mean difference, respectively.
#' @param method character string indicating the method to be used. Options are \code{cohens}, \code{hedges}, and \code{glass}.
#' @param data an optional data frame (or object coercible by \code{\link{as.data.frame}} to a data frame) containing the variables in the previous arguments.
#'
#' @return A list containing the following
#' \tabular{ll}{
#' \code{y} \tab mean or standardized mean differences for each treatment level,
#' included the referent one (0 by calculation).\cr
#' \code{v} \tab variances corresponding to the mean or standardized mean differences for each
#' treatment level, included the referent one (0 by calculation)\cr
#' \code{S} \tab co(variance) matrix for the non-referent mean or standardized mean differences.\cr
#' }
#'
#' @details This is an internal function called by \code{\link{dosresmeta}} to reconstruct the (co)variance matrix of the
#' outcome variable. The function is expected to be extended and/or modified at every release of the package
#'
#'
#' @examples
#' data("ari")
#'
#' ## Obtaining standardized mean differences, variances, and (co)varinace
#' ## matrix for the first study (id = 1)
#' covar.smd(y, sd, n, measure = "smd", data = subset(ari, id == 1))
#'
#' ## Obtaining mean differences, variances, and (co)varinace matrices for the all the studies
#' cov.md <- by(ari, ari\$id, function(x) covar.smd(y, sd, n, "md", data = x))
#'
#' ## Extracting mean differences
#' unlist(lapply(cov.md, function(x) x\$y))
#' ## Extracting variances for the mean differences
#' unlist(lapply(cov.md, function(x) x\$v))
#' ## List of the (co)variance matrices for the mean differences
#' lapply(cov.md, function(x) x\$S)
#'
#' @author Alessio Crippa, \email{[email protected]@ki.se}
#'
#' @references
#'
#' Cooper, H., Hedges, L. V., & Valentine, J. C. (Eds.). (2009). The handbook of
#' research synthesis and meta-analysis. Russell Sage Foundation.
#'
#' @export covar.smd

covar.smd <- function(y, sd, n, measure = "md", method = "cohens", data){
if (missing(data))
data <- NULL
if (is.null(data)) {
data <- sys.frame(sys.parent())
}
else {
if (!is.data.frame(data)) {
data <- data.frame(data)
}
}
method <- match.arg(method, c("cohens", "hedges", "glass"))
mf <- match.call()
mf.y <- mf[[match("y", names(mf))]]
mf.sd <- mf[[match("sd", names(mf))]]
mf.n <- mf[[match("n", names(mf))]]
y <- eval(mf.y, data, enclos = sys.frame(sys.parent()))
sd <- eval(mf.sd, data, enclos = sys.frame(sys.parent()))
n <- eval(mf.n, data, enclos = sys.frame(sys.parent()))
y <- scale(y, y, F)
sdPooled <- sqrt(sum((n-1)*sd^2)/sum(n-1))
if (measure == "md"){
v <- sdPooled^2 * (n + n) / (n * n)
v <- 0
S <- matrix(sd^2/n, ncol = length(v[v != 0]),
nrow = length(v[v != 0]))
diag(S) <- v[v != 0]
}
if (measure == "smd"){
if (method == "cohens"){
y <- y/sdPooled
S <- 1/n + tcrossprod(y[-1])/(2*sum(n))
}
if (method == "hedges"){
y <- (y/sdPooled) * (1 - (3/((4 * sum(n)) - 9)))
S <- 1/n + tcrossprod(y[-1])/(2 * (sum(n) - 3.94))
}
else if (method == "glass"){
y <- y/sd
S <- 1/n + tcrossprod(y[-1])/(2 * (n - 1))
}
v <- c(0, 1/n[-1] + diag(S))
diag(S) <- v[-1]
}
list(y = y, v = v, S = S)
}

#' Computes the covariance matrix for a set of log relative risks
#'
#' @description Reconstructs the covariance matrix for a set of (reported) log relative risks, given the number of cases and
#' the number of total persons or person-years for each treatment (dose) level.
#'
#' @inheritParams hamling
#' @param covariance method to approximate the coviariance among set of reported log relative risks, "\code{gl}" for the method proposed by Greenland and Longnecker
#' (default), "\code{h}" for the method proposed by Hamling.
#'
#' @return The (co)variance matrix of the log relative risks.
#'
#' @details This is an internal function called by \code{\link{dosresmeta}} to reconstruct the (co)variance matrix of the (adjusted) log relative risks. The function
#' calls, depending on the choosen method, \code{\link{grl}} (default) or \code{\link{hamling}} to reconstruct the effective counts corresponding to the (adjusted) log
#' relative risks as well as their standard errors. From these it computes the covariance matrix; analytical formulas can be found in the referenced article.
#'
#' @author Alessio Crippa, \email{[email protected]@ki.se}
#'
#' @references
#' Orsini, N., Li, R., Wolk, A., Khudyakov, P., Spiegelman, D. (2012). Meta-analysis for linear and nonlinear dose-response relations:
#' examples, an evaluation of approximations, and software. American journal of epidemiology, 175(1), 66-73.
#'
#'
#' @examples
#' data("alcohol_cvd")
#'
#' ## Obtaining the (co)variance matrix of log RR for the first study (id = 1)
#' covar.logrr(y = logrr, v = I(se^2), cases = cases, n = n, type = type,
#'             data = subset(alcohol_cvd, id == 1))
#'
#' ## Obtaining the (co)variance matrices of log RRfor all study
#' by(alcohol_cvd, alcohol_cvd\$id, function(x)
#'    covar.logrr(y = logrr, v = I(se^2), cases = cases, n = n,
#'                type = type, data = x))
#'
#' ## Restructuring the previous results in a list of matrices
#' do.call("list", by(alcohol_cvd, alcohol_cvd\$id, function(x)
#'    covar.logrr(y = logrr, v = I(se^2), cases = cases, n = n, type = type,
#'                data = x)))
#'
#'@export covar.logrr

covar.logrr <- function(cases, n, y, v, type, data, covariance = "gl"){
if (missing(data))
data <- NULL
if (is.null(data)) {
data <- sys.frame(sys.parent())
}
else {
if (!is.data.frame(data)) {
data <- data.frame(data)
}
}
covariance <- match.arg(covariance, c("gl", "hamling"))
mf <- match.call(expand.dots = FALSE)
mf.cases <- mf[[match("cases", names(mf))]]
mf.n <- mf[[match("n", names(mf))]]
mf.y <- mf[[match("y", names(mf))]]
mf.v <- mf[[match("v", names(mf))]]
mf.type <- mf[[match("type", names(mf))]]
cases <- eval(mf.cases, data, enclos = sys.frame(sys.parent()))
n <- eval(mf.n, data, enclos = sys.frame(sys.parent()))
y <- eval(mf.y, data, enclos = sys.frame(sys.parent()))
v <- eval(mf.v, data, enclos = sys.frame(sys.parent()))
v[is.na(v)] <- 0
type <- eval(mf.type, data, enclos = sys.frame(sys.parent()))
if (is.null(type))
type <- as.vector(mf.type)
ps <- if (covariance == "gl")
grl(y, v, cases, n, type)
else
hamling(y, v, cases, n, type)
rcorr <- switch(as.character(type),
cc = {
s0 <- 1/ps[v==0, 1] + 1/(ps[v==0, 2] - ps[v==0, 1])
si <- s0 + 1/ps[v!=0, 1] + 1/(ps[v!=0, 2] - ps[v!=0, 1])
s0/tcrossprod(si)^0.5
},
ir = {
s0 <- 1/ps[v==0, 1]
si <- s0 + 1/ps[v!=0, 1]
s0/tcrossprod(si)^0.5
},
ci = {
s0 <- 1/ps[v==0, 1] - 1/ps[v==0, 2]
si <- s0 + 1/ps[v!=0, 1] - 1/ps[v!=0, 2]
s0/tcrossprod(si)^0.5
})
diag(rcorr) <- 1
tcrossprod(v[v != 0])^0.5 * rcorr
}

#' Approximating effective-counts as proposed by Hamling
#'
#' @description Reconstructs the set of pseudo-numbers (or "effective" numbers) of cases and non-cases consistent
#' with the input data (log relative risks). The method was first proposed in 2008 by Hamling.
#'
#' @param y a vector, defining the (reported) log relative risks.
#' @param v a vector, defining the variances of the reported log relative risks.
#' @param cases a vector, defining the number of cases for each exposure level.
#' @param n a vector, defining the total number of subjects for each exposure level. For incidence-rate data \code{n} indicates the amount of person-time within
#' each exposure level.
#' @param type a vector (or a character string), specifying the design of the study. Options are
#' \code{cc}, \code{ir}, and \code{ci}, for case-control, incidence-rate, and cumulative incidence data, respectively.
#' @param data an optional data frame (or object coercible by \code{\link{as.data.frame}} to a data frame) containing the variables in the previous arguments.
#'
#' @return A list containing the following
#' \tabular{ll}{
#' \code{y} \tab mean or standardized mean differences for each treatment level,
#' included the referent one (0 by calculation).\cr
#' \code{v} \tab variances corresponding to the mean or standardized mean differences for each
#' treatment level, included the referent one (0 by calculation)\cr
#' \code{S} \tab co(variance) matrix for the non-referent mean or standardized mean differences.\cr
#' }
#'
#' @details The function reconstructs the effective counts corresponding to the multivariable adjusted log relative risks as well as their standard errors.
#' A unique solution is guaranteed by keeping the ratio non-cases to cases and the fraction of unexposed subjects equal to the unadjusted data (Hamling).
#' See the referenced article for a complete description of the algorithm implementation.
#'
#' @examples
#' data("alcohol_cvd")
#'
#' ## Obtaining pseudo-counts for the first study (id = 1)
#' hamling(y = logrr, v = I(se^2), cases = cases, n = n, type = type,
#' data = subset(alcohol_cvd, id == 1))
#'
#' ## Obtaining pseudo-counts for all study
#' by(alcohol_cvd, alcohol_cvd\$id, function(x)
#' hamling(y = logrr, v = I(se^2), cases = cases, n = n, type = type, data = x))
#'
#' ## Restructuring the previous results in a matrix
#' do.call("rbind", by(alcohol_cvd, alcohol_cvd\$id, function(x)
#'    hamling(y = logrr, v = I(se^2), cases = cases, n = n, type = type,
#'       data = x)))
#'
#' @author Alessio Crippa, \email{[email protected]@ki.se}
#'
#'
#' @references
#' Hamling, J., Lee, P., Weitkunat, R., Ambuhl, M. (2008). Facilitating meta-analyses by deriving relative effect and precision estimates for alternative
#' comparisons from a set of estimates presented by exposure level or disease category. Statistics in medicine, 27(7), 954-970.
#'
#' Orsini, N., Li, R., Wolk, A., Khudyakov, P., Spiegelman, D. (2012). Meta-analysis for linear and nonlinear dose-response relations: examples, an evaluation
#' of approximations, and software. American journal of epidemiology, 175(1), 66-73.
#'
#' @export hamling

hamling <- function(y, v, cases, n, type, data){
if (missing(data))
data <- NULL
if (is.null(data)) {
data <- sys.frame(sys.parent())
}
else {
if (!is.data.frame(data)) {
data <- data.frame(data)
}
}
mf <- match.call(expand.dots = FALSE)
mf.y <- mf[[match("y", names(mf))]]
mf.v <- mf[[match("v", names(mf))]]
mf.cases <- mf[[match("cases", names(mf))]]
mf.n <- mf[[match("n", names(mf))]]
mf.type <- mf[[match("type", names(mf))]]
y <- eval(mf.y, data, enclos = sys.frame(sys.parent()))
v <- eval(mf.v, data, enclos = sys.frame(sys.parent()))
v[is.na(v)] <- 0
cases <- eval(mf.cases, data, enclos = sys.frame(sys.parent()))
n <- eval(mf.n, data, enclos = sys.frame(sys.parent()))
type <- eval(mf.type, data, enclos = sys.frame(sys.parent()))
if (is.null(type))
type <- as.vector(mf.type)
p0 <- if (as.character(type) == "cc")
(n - cases)[v==0]/sum(n - cases)
else
n[v==0]/sum(n)
z0 <- if (as.character(type) == "cc")
sum((n - cases)[v!=0])/sum(cases[v!=0])
else
sum(n[v!=0])/sum(cases[v!=0])
init <- c(cases[v==0], n[v==0])
opt <- optim(init, fun.h, v = v, y = y, type = type, p0 = p0, z0 = z0)
pscounts <- est.ps.h(opt\$par, v, y, type)
pscounts
}

#' @noRd
est.ps.h <- function(param, v, y, type){
A0 <- param
N0 <- param
ps <- switch(as.character(type),
cc = cbind(A <- (1+exp(y)*A0/(N0 - A0))/(v-1/A0-1/(N0 - A0)),
N = A + (1+(N0 - A0)/(A0*exp(y)))/(v-1/A0-1/(N0 - A0))),
ir = cbind(A = (1 - exp(y)*A0/N0)/(v - 1/A0 + 1/N0),
N = (N0/(A0*exp(y))-1)/(v-1/A0 + 1/N0)),
ci = cbind(A = 1/(v - 1/A0),
N = (N0/(A0*exp(y)))/(v-1/A0)))
ps[v==0, ] <- param
colnames(ps) <- c("A", "N")
return(ps)
}

#' @noRd
fun.h <- function(par, v, y, type, p0, z0){
ps <- est.ps.h(par, v, y, type)
if (as.character(type) == "cc"){
p1 <- (ps[v==0, 2] - ps[v==0, 1])/(sum(ps[, 2] - ps[, 1]))
z1 <- sum(ps[v!=0, 2] - ps[v!=0, 1])/sum(ps[v!=0, 1])
}
else {
p1 <- ps[v==0, 2]/sum(ps[, 2])
z1 <- sum(ps[v!=0, 2])/sum(ps[v!=0, 1])
}
((p1 - p0)/p0)^2+((z1 - z0)/z0)^2
}

#' Approximating effective-counts as proposed by Greenland & Longnecker
#'
#' @description Reconstructs the set of pseudo-numbers (or 'effective' numbers) of cases and non-cases consistent
#' with the input data (log relative risks). The method was first proposed in 1992 by Greenland and Longnecker.
#'
#' @inheritParams hamling
#' @param tol define the tolerance.
#'
#' @details The function reconstructs the effective counts corresponding to the multivariable adjusted log relative risks as well as their standard errors.
#' A unique solution is guaranteed by keeping the margins of the table of pseudo-counts equal to the margins of the crude or unadjusted data
#' (Greenland and Longnecker 1992). See the referenced article for a complete description of the algorithm implementation.
#'
#' @return The results are returned structured in a matrix
#' \tabular{ll}{
#' \code{A} \tab approximated number of effective cases. \cr
#' \code{N} \tab approximated total number of effective subjects. \cr
#' }
#'
#' @examples
#' data("alcohol_cvd")
#'
#' ## Obtaining pseudo-counts for the first study (id = 1)
#' grl(y = logrr, v = I(se^2), cases = cases, n = n, type = type,
#'    data = subset(alcohol_cvd, id == 1))
#'
#' ## Obtaining pseudo-counts for all study
#' by(alcohol_cvd, alcohol_cvd\$id, function(x)
#'    grl(y = logrr, v = I(se^2), cases = cases, n = n, type = type, data = x))
#'
#' ## Restructuring the previous results in a matrix
#' do.call("rbind", by(alcohol_cvd, alcohol_cvd\$id, function(x)
#'    grl(y = logrr, v = I(se^2), cases = cases, n = n, type = type, data = x)))
#'
#' @author Alessio Crippa, \email{[email protected]@ki.se}
#'
#'
#' @references
#' Greenland, S., Longnecker, M. P. (1992). Methods for trend estimation from summarized dose-response data, with applications to meta-analysis. American journal of epidemiology, 135(11), 1301-1309.
#'
#' Orsini, N., Li, R., Wolk, A., Khudyakov, P., Spiegelman, D. (2012). Meta-analysis for linear and nonlinear dose-response relations: examples, an evaluation of approximations, and software.
#' American journal of epidemiology, 175(1), 66-73.
#'
#' @export grl

grl <- function(y, v, cases, n, type, data, tol = 1e-05){
if (missing(data))
data <- NULL
if (is.null(data)) {
data <- sys.frame(sys.parent())
}
else {
if (!is.data.frame(data)) {
data <- data.frame(data)
}
}
mf <- match.call(expand.dots = FALSE)
mf.y <- mf[[match("y", names(mf))]]
mf.v <- mf[[match("v", names(mf))]]
mf.cases <- mf[[match("cases", names(mf))]]
mf.n <- mf[[match("n", names(mf))]]
mf.type <- mf[[match("type", names(mf))]]
y <- eval(mf.y, data, enclos = sys.frame(sys.parent()))
v <- eval(mf.v, data, enclos = sys.frame(sys.parent()))
v[is.na(v)] <- 0
cases <- eval(mf.cases, data, enclos = sys.frame(sys.parent()))
n <- eval(mf.n, data, enclos = sys.frame(sys.parent()))
type <- eval(mf.type, data, enclos = sys.frame(sys.parent()))
if (is.null(type))
type <- as.vector(mf.type)
Ax <- Axp <- cases
repeat{
A0 <- sum(cases) - sum(Ax[v!=0])
cx <- if (!as.character(type) == "ir")
1/Ax + 1/(n - Ax)
else
1/Ax
e <- if (!as.character(type) == "ir")
y[v!=0] + log(A0) + log(n[v!=0]-Ax[v!=0]) - log(Ax[v!=0]) -
log(n[v==0] - A0)
else
y[v!=0] + log(A0) + log(n[v!=0]) - log(Ax[v!=0]) - log(n[v==0])
H <- diag(cx[v!=0] + cx[v==0], nrow = sum(v!=0))
H[upper.tri(H)] <- H[lower.tri(H)] <- cx[v==0]
Axp[v==0] <- A0
Axp[v!=0] <- Ax[v!=0] + solve(H) %*% e
delta <- sum((Axp[v!=0] - Ax[v!=0])^2)
if (delta < tol)
break
Ax <- Axp
}
cbind(A = Axp, N = n)
}

#' @noRd
change_ref <- function(y, v, cases, n, type, data, ref = 1,
method = "hamling", expo = FALSE){
if (missing(data))
data <- NULL
if (is.null(data)) {
data <- sys.frame(sys.parent())
}
else {
if (!is.data.frame(data)) {
data <- data.frame(data)
}
}
mf <- match.call(expand.dots = FALSE)
mf.y <- mf[[match("y", names(mf))]]
mf.v <- mf[[match("v", names(mf))]]
mf.cases <- mf[[match("cases", names(mf))]]
mf.n <- mf[[match("n", names(mf))]]
mf.type <- mf[[match("type", names(mf))]]
y <- eval(mf.y, data, enclos = sys.frame(sys.parent()))
v <- eval(mf.v, data, enclos = sys.frame(sys.parent()))
v[is.na(v)] <- 0
cases <- eval(mf.cases, data, enclos = sys.frame(sys.parent()))
n <- eval(mf.n, data, enclos = sys.frame(sys.parent()))
type <- eval(mf.type, data, enclos = sys.frame(sys.parent()))
if (is.null(type))
type <- as.vector(mf.type)
ps <- if (method == "hamling"){
hamling(y, v, cases, n, type)
} else if (method == "gl"){
grl(y, v, cases, n, type)
}
rr <- if (type == "cc"){
(ps[, 1]*(ps[ref, 2] - ps[ref, 1]))/(ps[ref, 1]*(ps[, 2] - ps[, 1]))
} else {
(ps[, 1]/ps[, 2])/(ps[ref, 1]/ps[ref, 2])
}
var_logrr <- if (type == "cc"){
1/ps[, 1] + 1/(ps[, 2] - ps[, 1]) + 1/ps[ref, 1] + 1/(ps[ref, 2] - ps[ref, 1])
} else if (type == "ci") {
1/ps[, 1] - 1/ps[, 2] + 1/ps[ref, 1] - 1/ps[ref, 2]
} else if (type == "ir") {
1/ps[, 1] + 1/ps[, 2]
}
var_logrr[ref] <- 0
lb_rr <- exp(log(rr) - qnorm(.975)*sqrt(var_logrr))
ub_rr <- exp(log(rr) + qnorm(.975)*sqrt(var_logrr))
dat <- if (expo){
data.frame(ps, rr = rr, lb_rr = lb_rr, ub_rr = ub_rr)
} else {
data.frame(ps, logrr = log(rr), v = var_logrr,
loglbrr = log(lb_rr), logubrr = log(ub_rr))
}
colnames(dat) <- paste(colnames(dat), ref, sep = ".")
return(dat)
}
```

## Try the dosresmeta package in your browser

Any scripts or data that you put into this service are public.

dosresmeta documentation built on May 2, 2019, 6:30 a.m.