a first 'beta' version

man/RobergeEtAl2017.Rd

@@ -12,7 +12,7 @@
     \bold{year}            \tab \code{numeric}   \tab publication year \cr
     \bold{N}               \tab \code{numeric}   \tab number of patients \cr
     \bold{onset}           \tab \code{factor}    \tab treatment onset (up to 16 weeks' gestation or later) \cr
-    \bold{dose}            \tab \code{numeric}   \tab dose (mg) \cr
+    \bold{dose}            \tab \code{numeric}   \tab dose (mg/day) \cr
     \bold{control}         \tab \code{factor}    \tab type of control group \cr
     \bold{asp.PE.events}   \tab \code{numeric}   \tab number of PE events in aspirin group \cr
     \bold{asp.PE.total}    \tab \code{numeric}   \tab number of PE cases in aspirin group \cr
@@ -25,12 +25,13 @@
   }
 }
 \details{A systematic literature review was performed in order to
-  investigate effects of aspirin administered during pregnancy. Of
-  particular interest were the occurrence of \emph{preeclampsia (PE)}
-  and \emph{fetal growth restriction (FGR)}. A total of 45 relevant
-  studies were found, out of which 40 reported on PE, and 30 reported on
-  FGR. Besides event rates, the mode of administration (treatment onset
-  (early vs. late) and dose (in mg)) was recorded.
+  summarize the evidence on effects of aspirin administered during
+  pregnancy. Of particular interest were occurrences of
+  \emph{preeclampsia (PE)} and \emph{fetal growth restriction (FGR)}. A
+  total of 45 relevant randomized controlled trials (RCTs) were found,
+  out of which 40 reported on PE, and 30 reported on FGR. Besides event
+  rates, the mode of administration (treatment onset (early vs. late)
+  and dose (in mg)) was also recorded for each study.
 }
 \source{S. Roberge, K. Nicolaides, S. Demers, J. Hyett, N. Chaillet, E. Bujold.
   \href{https://doi.org/10.1016/j.ajog.2016.09.076}{The role of aspirin
@@ -58,10 +59,11 @@ es.pe  <- escalc(measure="OR",
                  subset=complete.cases(RobergeEtAl2017[,7:10]))
 # show forest plot:
 forestplot(es.pe, title="preeclampsia (PE)")
-# show "bubble plot":
+# show "bubble plot" (bubble sizes are
+# inversely proportional to standard errors):
 plot(es.pe$dose, es.pe$yi, cex=1/sqrt(es.pe$vi),
      col=c("blue","red")[as.numeric(es.pe$onset)],
-     xlab="dose", ylab="log-OR (PE)", main="Roberge et al. (2017)")
+     xlab="dose (mg)", ylab="log-OR (PE)", main="Roberge et al. (2017)")
 legend("topright", col=c("blue","red"), c("early onset", "late onset"), pch=1)
 
 # fetal growth restriction (FGR):
@@ -75,7 +77,7 @@ forestplot(es.fgr, title="fetal growth restriction (FGR)")
 # show "bubble plot":
 plot(es.fgr$dose, es.fgr$yi, cex=1/sqrt(es.fgr$vi),
      col=c("blue","red")[as.numeric(es.fgr$onset)],
-     xlab="dose", ylab="log-OR (FGR)", main="Roberge et al. (2017)")
+     xlab="dose (mg)", ylab="log-OR (FGR)", main="Roberge et al. (2017)")
 legend("topright", col=c("blue","red"), c("early onset", "late onset"), pch=1)
 
 \dontrun{
@@ -97,6 +99,28 @@ print(X02)
 # perform regression:
 bmr02 <- bmr(es.pe, X=X02)
 bmr02$summary
+
+# derive predictions from the model;
+# specify corresponding "regressor matrices":
+newx.early <- cbind(1, 0, seq(50, 150, by=5), 0)
+newx.late  <- cbind(0, 1, 0, seq(50, 150, by=5))
+# (note: columns correspond to "beta" parameters)
+
+# compute predicted medians and 95 percent intervals: 
+pred.early <- cbind("median"=bmr02$qpred(0.5, x=newx.early),
+                    bmr02$pred.interval(x=newx.early))
+pred.late <- cbind("median"=bmr02$qpred(0.5, x=newx.late),
+                    bmr02$pred.interval(x=newx.late))
+
+# draw "bubble plot": 
+plot(es.pe$dose, es.pe$yi, cex=1/sqrt(es.pe$vi),
+     col=c("blue","red")[as.numeric(es.pe$onset)],
+     xlab="dose (mg)", ylab="log-OR (PE)", main="Roberge et al. (2017)")
+legend("topright", col=c("blue","red"), c("early onset", "late onset"), pch=1)
+# add predictions to bubble plot:
+matlines(newx.early[,3], pred.early, col="blue", lty=c(1,2,2))
+matlines(newx.late[,4], pred.late, col="red", lty=c(1,2,2))
+
 }
 }
 \keyword{datasets}

man/bmr.Rd

@@ -36,7 +36,7 @@
 }
 \arguments{
   \item{y}{
-    vector of estimates \emph{or} an \code{\link[metafor]{escalc}} object.
+    vector of estimates, \emph{or} an \code{\link[metafor]{escalc}} object.
 }
   \item{sigma}{
     vector of standard errors associated with \code{y}.
@@ -62,9 +62,9 @@
       \item \code{"shrinkage"}: the \sQuote{uniform shrinkage} prior
       \item \code{"I2"}: a uniform prior on the \sQuote{relative heterogeneity} \eqn{I^2}
     }
-    The default is \code{"uniform"} (which should be used with caution;
-    see remarks below). The above priors are described in some more
-    detail below.
+    The default is \code{"uniform"} (which should be used with
+    caution). The above priors are described in some more detail in the
+    \code{\link{bayesmeta}()} help.
 }
   \item{beta.prior.mean, beta.prior.sd, beta.prior.cov}{
     the mean and standard deviations, or covariance of the normal prior
@@ -104,7 +104,7 @@
   \eqn{i}-th estimate, that is associated with standard error
   \eqn{\sigma_i > 0}{sigma[i] > 0}, where \eqn{i=1,...,k}. In addition to
   estimates and standard errors for the \eqn{i}th observation,
-  covariabes \eqn{x_{i,j}}{x[i,j]} with \eqn{j=1,...,d} are available
+  a set of covariables \eqn{x_{i,j}}{x[i,j]} with \eqn{j=1,...,d} are available
   for each estimate \eqn{y_i}{y[i]}.
 
   The model includes \eqn{d+1} unknown parameters,
@@ -140,17 +140,19 @@
 
   Meta-regression allows the consideration of (study-level) covariables
   in a meta-analysis. Quite often, these may also be indicator variables
-  (\sQuote{zero/one} variables) idenfying subgroups of studies.
+  (\sQuote{zero/one} variables) simply identifying subgroups of studies.
 
 \subsection{Connection to the simple random-effects model}{
   The meta-regression model is a generalisation of the \sQuote{simple}
-  random-effects model that is implemented e.g. in the
+  random-effects model that is implemented in the
   \code{\link{bayesmeta}()} function. Meta-regression reduces to the
   estimation of a single \dQuote{intercept} term when the regressor
   matrix (\eqn{X}) has a single column (\eqn{d=1}) and consists only of
-  ones. The single regression coefficient \eqn{\beta_1}{beta[1]} then is
-  equivalent to the \eqn{\mu} parameter from the simple random effects
-  model (see also the \sQuote{Examples} section below).
+  ones (which is also the default setting in case the \sQuote{\code{X}}
+  argument is left unspecified). The single regression coefficient
+  \eqn{\beta_1}{beta[1]} then is equivalent to the \eqn{\mu} parameter
+  from the simple random effects model (see also the \sQuote{Examples}
+  section below).
   }
 
 \subsection{Prior specification}{
@@ -165,16 +167,16 @@
 
   For sensible prior choices for the heterogeneity parameter \eqn{\tau},
   see also Roever (2020), Roever \emph{et al.} (2021) and the
-  \code{\link{bayesmeta}()} help.
+  \sQuote{\code{\link{bayesmeta}()}} function's help.
   }
   
 \subsection{Computation}{
   The \code{bmr()} function utilizes the same computational method
   as the \code{\link{bayesmeta}()} function to derive the posterior
-  distribution, namely, the \acronym{DIRECT} algorithm (Roever and Friede,
-  2017). The model as outlined above may be adapted by specifying
-  appropriate design matrices \eqn{X} and prior distributions for the
-  effect parameter vector \eqn{\beta} and the heterogeneity \eqn{\tau}.
+  distribution, namely, the \acronym{DIRECT} algorithm. Numerical
+  accuracy of the computations is determined by the \sQuote{\code{delta}}
+  and \sQuote{\code{epsilon}} arguments (Roever and Friede,
+  2017).
   }
 }
 \value{
@@ -186,16 +188,17 @@
   \item{k}{number of data points (length of \code{y}, or rows of \eqn{X}).}
   \item{d}{number of coefficients (columns of \code{X}).}
   \item{labels}{vector of labels corresponding to \code{y} and \code{sigma}.}
-  \item{variables}{variable names for the \eqn{\beta} coefficients.}
+  \item{variables}{variable names for the \eqn{\beta} coefficients
+  (determined by the column names of the supplied \code{X} argument).}
   \item{tau.prior}{the prior probability density function for \eqn{\tau}.}
   \item{tau.prior.proper}{a \code{logical} flag indicating whether the
     heterogeneity prior appears to be proper (which is judged based on
     an attempted numerical integration of the density function).} 
-  \item{beta.prior}{a \code{list} containing the prior mean and
-  covariance matrix for the coefficient vector \eqn{\beta}.}
+  \item{beta.prior}{a \code{list} containing the prior mean vector and
+  covariance matrix for the coefficients \eqn{\beta}.}
   \item{beta.prior.proper}{a \code{logical} vector (of length \eqn{d})
-  indicating whether the corresponding \eqn{beta} coefficient's prior is
-  proper (finite prior mean and variance specified).} 
+  indicating whether the corresponding \eqn{\beta} coefficient's prior is
+  proper (i.e., finite prior mean and variance were specified).} 
   \item{dprior}{a \code{function(tau, beta, which.beta, log=FALSE)} of
   \eqn{\tau} and/or \eqn{\beta} parameters, returning either the joint or
     marginal prior probability density, depending on which parameter(s)
@@ -223,10 +226,10 @@
   corresponding to the supplied \eqn{y_i}{y[i]} data values
   (\eqn{i=1,\ldots,k}{i=1,...,k}). May be identified using the
   \sQuote{\code{which}} argument via its index (\eqn{i}) or a character
-  sting giving the corresponding study label.}
+  string giving the corresponding study label.}
   \item{cond.moment}{a \code{function(tau)}
     returning conditional posterior moments (mean and covariance) of
-    the \eqn{\beta} as a function of \eqn{\tau}.}
+    \eqn{\beta} as a function of \eqn{\tau}.}
   \item{summary}{a \code{matrix} listing some summary statistics, namely
     marginal posterior mode, median, mean, standard deviation
     and a (shortest) 95\% credible intervals,
@@ -286,10 +289,15 @@
   Christian Roever \email{christian.roever@med.uni-goettingen.de}
 }
 \seealso{
-  \code{\link{bayesmeta}}, \code{\link[metafor]{escalc}}.
+  \code{\link{bayesmeta}}, \code{\link[metafor]{escalc}},
+  \code{\link[stats]{model.matrix}}, \code{\link{CrinsEtAl2014}},
+  \code{\link{RobergeEtAl2017}}.
 }
 \examples{
 \dontrun{
+######################################################################
+# (1)  A simple example with two groups of studies
+#
 # load data:
 data("CrinsEtAl2014")
 # compute effect measures (log-OR):
@@ -304,20 +312,118 @@ crins.es[,c("publication", "randomized", "exp.AR.events", "exp.total",
 # specify regressor matrix:
 X <- cbind("intercept"=1, "nrand"=(crins.es$randomized=="no"))
 print(X)
-# NB: regressor matrix specifes an overall intercept
+# NB: regressor matrix specifies an overall intercept
 #     and an offset term for non-randomized studies.
 
 # perform regression:
 bmr01 <- bmr(y=crins.es$yi, sigma=sqrt(crins.es$vi),
              labels=crins.es$publication, X=X)
 
-# alternatively, may simply supply the "escalc" object:
+# alternatively, may simply supply the "escalc" object
+# (yields identical results):
 bmr01 <- bmr(crins.es, X=X)
 
+# show results:
 bmr01
 bmr01$summary
 plot(bmr01)
 pairs(bmr01)
+
+# NOTE: there are many ways to set up the regressor matrix "X".
+# See the above specification and check out the following alternatives:
+X <- cbind("intercept"=1.0, "offset"=0.5*c(-1,1,1,1,-1,1))
+X <- cbind("rand"=c(1,0,0,0,1,0), "nrand"=c(0,1,1,1,0,1))
+
+
+######################################################################
+# (2)  A simple example reproducing a "bayesmeta" analysis:
+
+data("CrinsEtAl2014")
+crins.es <- escalc(measure="OR",
+                   ai=exp.AR.events,  n1i=exp.total,
+                   ci=cont.AR.events, n2i=cont.total,
+                   slab=publication, data=CrinsEtAl2014)
+
+# a "simple" meta-analysis:
+bma02 <- bayesmeta(crins.es,
+                   tau.prior=function(t){dhalfnormal(t, scale=0.5)},
+                   mu.prior.mean=0, mu.prior.sd=4)
+
+# the equivalent "intercept-only" meta-regression:
+bmr02 <- bmr(crins.es,
+             tau.prior=function(t){dhalfnormal(t, scale=0.5)},
+             beta.prior.mean=0, beta.prior.sd=4)
+# the corresponding (default) regressor matrix:
+bmr02$X
+
+# compare computation time and numbers of bins used internally:
+cbind("seconds" = c("bayesmeta" = unname(bma02$init.time),
+                    "bmr"       = unname(bmr02$init.time)),
+      "bins"    = c("bayesmeta" = nrow(bma02$support),
+                    "bmr"       = nrow(bmr02$support$tau)))
+
+# compare heterogeneity estimates:
+rbind("bayesmeta"=bma02$summary[,1], "bmr"=bmr02$summary[,1])
+
+# compare effect estimates:
+rbind("bayesmeta"=bma02$summary[,2], "bmr"=bmr02$summary[,2])
+
+
+######################################################################
+# (3)  An example with binary as well as continuous covariables:
+#
+# load data:
+?RobergeEtAl2017
+data("RobergeEtAl2017")
+str(RobergeEtAl2017)
+head(RobergeEtAl2017)
+
+# compute effect sizes (log odds ratios) from count data:
+es.pe  <- escalc(measure="OR",
+                 ai=asp.PE.events,  n1i=asp.PE.total,
+                 ci=cont.PE.events, n2i=cont.PE.total,
+                 slab=study, data=RobergeEtAl2017,
+                 subset=complete.cases(RobergeEtAl2017[,7:10]))
+
+# show "bubble plot" (bubble sizes are
+# inversely proportional to standard errors):
+plot(es.pe$dose, es.pe$yi, cex=1/sqrt(es.pe$vi),
+     col=c("blue","red")[as.numeric(es.pe$onset)],
+     xlab="dose (mg)", ylab="log-OR (PE)", main="Roberge et al. (2017)")
+legend("topright", col=c("blue","red"), c("early onset", "late onset"), pch=1)
+
+# set up regressor matrix:
+# (individual intercepts and slopes for two subgroups):
+X <- model.matrix(~ -1 + onset + onset:dose, data=es.pe)
+colnames(X) <- c("intEarly", "intLate", "slopeEarly", "slopeLate")
+# check out regressor matrix (and compare to original data):
+print(X)
+
+# perform regression:
+bmr03 <- bmr(es.pe, X=X)
+bmr03$summary
+
+# derive predictions from the model;
+# specify corresponding "regressor matrices":
+newx.early <- cbind(1, 0, seq(50, 150, by=5), 0)
+newx.late  <- cbind(0, 1, 0, seq(50, 150, by=5))
+# (note: columns correspond to "beta" parameters)
+
+# compute predicted medians and 95 percent intervals: 
+pred.early <- cbind("median"=bmr03$qpred(0.5, x=newx.early),
+                    bmr03$pred.interval(x=newx.early))
+pred.late <- cbind("median"=bmr03$qpred(0.5, x=newx.late),
+                    bmr03$pred.interval(x=newx.late))
+
+# draw "bubble plot": 
+plot(es.pe$dose, es.pe$yi, cex=1/sqrt(es.pe$vi),
+     col=c("blue","red")[as.numeric(es.pe$onset)],
+     xlab="dose (mg)", ylab="log-OR (PE)", main="Roberge et al. (2017)")
+legend("topright", col=c("blue","red"), c("early onset", "late onset"), pch=1)
+# add predictions to bubble plot:
+matlines(newx.early[,3], pred.early, col="blue", lty=c(1,2,2))
+matlines(newx.late[,4], pred.late, col="red", lty=c(1,2,2))
+
 }
 }
 \keyword{ models }